text:
module MPS

import ..MathOptFormat

import MathOptInterface
const MOI = MathOptInterface

MOI.Utilities.@model(InnerModel,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (MOI.SOS1, MOI.SOS2),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    ()
)

struct Options
    warn::Bool
end

get_options(m::InnerModel) = get(m.ext, :MPS_OPTIONS, Options(false))

"""
    Model(; kwargs...)

Create an empty instance of MathOptFormat.MPS.Model.

Keyword arguments are:

 - `warn::Bool=false`: print a warning when variables or constraints are renamed.
"""
function Model(;
        warn::Bool = false
)
    model = InnerModel{Float64}()
    model.ext[:MPS_OPTIONS] = Options(warn)
    return model
end

function Base.show(io::IO, ::InnerModel)
    print(io, "A Mathematical Programming System (MPS) model")
    return
end

# ==============================================================================
#
#   MOI.write_to_file
#
# ==============================================================================

function MOI.write_to_file(model::InnerModel, io::IO)
    options = get_options(model)
    MathOptFormat.create_unique_names(
        model, warn = options.warn, replacements = [' ' => '_']
    )
    write_model_name(io, model)
    write_rows(io, model)
    discovered_columns = write_columns(io, model)
    write_rhs(io, model)
    write_ranges(io, model)
    write_bounds(io, model, discovered_columns)
    write_sos(io, model)
    println(io, "ENDATA")
    return
end

# ==============================================================================
#   Model name
# ==============================================================================

function write_model_name(io::IO, model::InnerModel)
    model_name = MOI.get(model, MOI.Name())
    println(io, "NAME          ", model_name)
    return
end

# ==============================================================================
#   ROWS
# ==============================================================================

const LINEAR_CONSTRAINTS = (
    (MOI.LessThan{Float64}, 'L'),
    (MOI.GreaterThan{Float64}, 'G'),
    (MOI.EqualTo{Float64}, 'E'),
    (MOI.Interval{Float64}, 'L')  # See the note in the RANGES section.
)

function write_rows(io::IO, model::InnerModel)
    println(io, "ROWS\n N  OBJ")
    for (set_type, sense_char) in LINEAR_CONSTRAINTS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{
                                        MOI.ScalarAffineFunction{Float64},
                                        set_type}())
            row_name = MOI.get(model, MOI.ConstraintName(), index)
            if row_name == ""
                error("Row name is empty: $(index).")
            end
            println(io, " ", sense_char, "  ", row_name)
        end
    end
    return
end

# ==============================================================================
#   COLUMNS
# ==============================================================================

function list_of_integer_variables(model::InnerModel)
    integer_variables = Set{String}()
    for set_type in (MOI.ZeroOne, MOI.Integer)
        for index in MOI.get(model, MOI.ListOfConstraintIndices{
                MOI.SingleVariable, set_type}())
            v_index = MOI.get(model, MOI.ConstraintFunction(), index)
            v_name = MOI.get(model, MOI.VariableName(), v_index.variable)
            push!(integer_variables, v_name)
        end
    end
    return integer_variables
end

function add_coefficient(coefficients, variable_name, row_name, coefficient)
    if haskey(coefficients, variable_name)
        push!(coefficients[variable_name], (row_name, coefficient))
    else
        coefficients[variable_name] = [(row_name, coefficient)]
    end
    return
end

function extract_terms(
        model::InnerModel, coefficients, row_name::String,
        func::MOI.ScalarAffineFunction, discovered_columns::Set{String})
    for term in func.terms
        variable_name = MOI.get(model, MOI.VariableName(), term.variable_index)
        add_coefficient(coefficients, variable_name, row_name, term.coefficient)
        push!(discovered_columns, variable_name)
    end
    return
end

function extract_terms(
        model::InnerModel, coefficients, row_name::String, func::MOI.SingleVariable,
        discovered_columns::Set{String})
    variable_name = MOI.get(model, MOI.VariableName(), func.variable)
    add_coefficient(coefficients, variable_name, row_name, 1.0)
    push!(discovered_columns, variable_name)
    return
end

function write_columns(io::IO, model::InnerModel)
    # Many MPS readers (e.g., CPLEX and GAMS) will error if a variable (column)
    # appears in the BOUNDS section but did not appear in the COLUMNS section.
    # This is likely because such variables are meaningless - they don't appear
    # in the objective or constraints and so can be trivially removed.
    # To avoid generating MPS files that crash existing readers, we cache the
    # names of all variables seen in COLUMNS into `discovered_columns`, and then
    # pass this set to `write_bounds` so that it can act appropriately.
    discovered_columns = Set{String}()
    println(io, "COLUMNS")
    coefficients = Dict{String, Vector{Tuple{String, Float64}}}()
    for (set_type, sense_char) in LINEAR_CONSTRAINTS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{
                                        MOI.ScalarAffineFunction{Float64},
                                        set_type}())
            row_name = MOI.get(model, MOI.ConstraintName(), index)
            func = MOI.get(model, MOI.ConstraintFunction(), index)
            extract_terms(
                model, coefficients, row_name, func, discovered_columns)
        end
    end
    obj_func_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_func_type}())
    extract_terms(model, coefficients, "OBJ", obj_func, discovered_columns)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        # MPS doesn't support maximization so we flip the sign on the objective
        # coefficients.
        for (v_name, terms) in coefficients
            for (idx, (row_name, coef)) in enumerate(terms)
                if row_name == "OBJ"
                    terms[idx] = (row_name, -coef)
                end
            end
        end
    end
    integer_variables = list_of_integer_variables(model)
    for (variable, terms) in coefficients
        if variable in integer_variables
            println(io, "    MARKER    'MARKER'                 'INTORG'")
        end
        for (constraint, coefficient) in terms
            print(io, "     ", rpad(variable, 8), " ", rpad(constraint, 8), " ")
            Base.Grisu.print_shortest(io, coefficient)
            println(io)
        end
        if variable in integer_variables
            println(io, "    MARKER    'MARKER'                 'INTEND'")
        end
    end
    return discovered_columns
end

# ==============================================================================
#   RHS
# ==============================================================================

value(set::MOI.LessThan) = set.upper
value(set::MOI.GreaterThan) = set.lower
value(set::MOI.EqualTo) = set.value
value(set::MOI.Interval) = set.upper  # See the note in the RANGES section.

function write_rhs(io::IO, model::InnerModel)
    println(io, "RHS")
    for (set_type, sense_char) in LINEAR_CONSTRAINTS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{
                                        MOI.ScalarAffineFunction{Float64},
                                        set_type}())
            row_name = MOI.get(model, MOI.ConstraintName(), index)
            print(io, "    rhs       ", rpad(row_name, 8), "  ")
            set = MOI.get(model, MOI.ConstraintSet(), index)
            Base.Grisu.print_shortest(io, value(set))
            println(io)
        end
    end
    return
end

# ==============================================================================
#  RANGES
#
# Here is how RANGE information is encoded.
#
#     Row type | Range value |  lower bound  |  upper bound
#     ------------------------------------------------------
#         G    |     +/-     |     rhs       | rhs + |range|
#         L    |     +/-     | rhs - |range| |     rhs
#         E    |      +      |     rhs       | rhs + range
#         E    |      -      | rhs + range   |     rhs
#
# We elect to write out ScalarAffineFunction-in-Interval constraints in terms of
# LessThan (L) constraints with a range shift. The RHS term is set to the upper
# bound, and the RANGE term to upper - lower.
# ==============================================================================

function write_ranges(io::IO, model::InnerModel)
    println(io, "RANGES")
    for index in MOI.get(model, MOI.ListOfConstraintIndices{
                                    MOI.ScalarAffineFunction{Float64},
                                    MOI.Interval{Float64}}())
        row_name = MOI.get(model, MOI.ConstraintName(), index)
        print(io, "    rhs       ", rpad(row_name, 8), "  ")
        set = MOI.get(model, MOI.ConstraintSet(), index)::MOI.Interval{Float64}
        Base.Grisu.print_shortest(io, set.upper - set.lower)
        println(io)
    end
    return
end

# ==============================================================================
#   BOUNDS
# Variables default to [0, ∞).
#     FR    free variable      -∞ < x < ∞
#     FX    fixed variable     x == b
#     MI    lower bound -inf   -∞ < x
#     UP    upper bound             x <= b
#     PL    upper bound +inf        x < ∞
#     LO    lower bound        b <= x
#     BV    binary variable    x = 0 or 1
#
#  Not yet implemented:
#     LI    integer variable   b <= x
#     UI    integer variable        x <= b
#     SC    semi-cont variable x = 0 or l <= x <= b
#           l is the lower bound on the variable. If none set then defaults to 1
# ==============================================================================

function write_single_bound(io::IO, var_name::String, lower, upper)
    name = rpad(var_name, 8)
    if lower == upper
        print(io, " FX bounds    ", name, " ")
        Base.Grisu.print_shortest(io, lower)
        println(io)
    elseif lower == -Inf && upper == Inf
        # Skip this for now, we deal with it at the end of write_bounds.
    else
        if lower == -Inf
            println(io, " MI bounds    ", name)
        else
            print(io, " LO bounds    ", name, " ")
            Base.Grisu.print_shortest(io, lower)
            println(io)
        end
        if upper == Inf
            println(io, " PL bounds    ", name)
        else
            print(io, " UP bounds    ", name, " ")
            Base.Grisu.print_shortest(io, upper)
            println(io)
        end
    end
    return
end

function update_bounds(current::Tuple{Float64, Float64}, set::MOI.GreaterThan)
    return (max(current[1], set.lower), current[2])
end

function update_bounds(current::Tuple{Float64, Float64}, set::MOI.LessThan)
    return (current[1], min(current[2], set.upper))
end

function update_bounds(::Tuple{Float64, Float64}, set::MOI.Interval)
    return (set.lower, set.upper)
end

function update_bounds(::Tuple{Float64, Float64}, set::MOI.EqualTo)
    return (set.value, set.value)
end

function write_bounds(io::IO, model::InnerModel, discovered_columns::Set{String})
    println(io, "BOUNDS")
    free_variables = Set(MOI.get(model, MOI.ListOfVariableIndices()))
    bounds = Dict{MOI.VariableIndex, Tuple{Float64, Float64}}()
    for (set_type, sense_char) in LINEAR_CONSTRAINTS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{
                MOI.SingleVariable, set_type}())
            func = MOI.get(model, MOI.ConstraintFunction(), index)
            variable_index = func.variable::MOI.VariableIndex
            if !haskey(bounds, variable_index)
                bounds[variable_index] = (-Inf, Inf)
            end
            set = MOI.get(model, MOI.ConstraintSet(), index)::set_type
            bounds[variable_index] = update_bounds(bounds[variable_index], set)
        end
    end
    for (index, (lower, upper)) in bounds
        var_name = MOI.get(model, MOI.VariableName(), index)
        if var_name in discovered_columns
            write_single_bound(io, var_name, lower, upper)
        else
            @warn("Variable $var_name is mentioned in BOUNDS, but is not " *
                  "mentioned in the COLUMNS section. We are ignoring it.")
        end
        pop!(free_variables, index)
    end
    for index in MOI.get(model, MOI.ListOfConstraintIndices{
            MOI.SingleVariable, MOI.ZeroOne}())
        func = MOI.get(model, MOI.ConstraintFunction(), index)
        variable_index = func.variable::MOI.VariableIndex
        var_name = MOI.get(model, MOI.VariableName(), variable_index)
        if var_name in discovered_columns
            println(io, " BV bounds    ", var_name)
        else
            @warn("Variable $var_name is mentioned in BOUNDS, but is not " *
                  "mentioned in the COLUMNS section. We are ignoring it.")
        end
        if variable_index in free_variables
            # We can remove the variable because it has a bound, but first check
            # that it is still there because some variables might have two
            # bounds and so might have already been removed.
            pop!(free_variables, variable_index)
        end
    end
    for variable_index in free_variables
        var_name = MOI.get(model, MOI.VariableName(), variable_index)
        if var_name in discovered_columns
            println(io, " FR bounds    ", var_name)
        else
            @warn("Variable $var_name is mentioned in BOUNDS, but is not " *
                  "mentioned in the COLUMNS section. We are ignoring it.")
        end
    end
    return
end

# ==============================================================================
#   SOS
# ==============================================================================

function write_sos_constraint(io::IO, model::InnerModel, index)
    func = MOI.get(model, MOI.ConstraintFunction(), index)
    set = MOI.get(model, MOI.ConstraintSet(), index)
    for (variable, weight) in zip(func.variables, set.weights)
        var_name = MOI.get(model, MOI.VariableName(), variable)
        print(io, "    ", rpad(var_name, 8), "  ")
        Base.Grisu.print_shortest(io, weight)
        println(io)
    end
end

function write_sos(io::IO, model::InnerModel)
    sos1_indices = MOI.get(model, MOI.ListOfConstraintIndices{
        MOI.VectorOfVariables, MOI.SOS1{Float64}}())
    sos2_indices = MOI.get(model, MOI.ListOfConstraintIndices{
        MOI.VectorOfVariables, MOI.SOS2{Float64}}())
    if length(sos1_indices) + length(sos2_indices) > 0
        println(io, "SOS")
        idx = 1
        for (sos_type, indices) in enumerate([sos1_indices, sos2_indices])
            for index in indices
                println(io, " S", sos_type, " SOS", idx)
                write_sos_constraint(io, model, index)
                idx += 1
            end
        end
    end
    return
end

# ==============================================================================
#
#   MOI.read_from_file
#
# Here is a template for an MPS file, reproduced from
# http://lpsolve.sourceforge.net/5.5/mps-format.htm.
#
# Field:    1           2          3         4         5         6
# Columns:  2-3        5-12      15-22     25-36     40-47     50-61
#           NAME   problem name
#           ROWS
#            type     name
#           COLUMNS
#                     column      row       value     row      value
#                     name        name                name
#           RHS
#                     rhs         row       value     row      value
#                     name        name                name
#           RANGES
#                     range       row       value     row      value
#                     name        name                name
#           BOUNDS
#            type     bound       column    value
#                     name        name
#           SOS
#            type     CaseName    SOSName   SOSpriority
#                     CaseName    VarName1  VarWeight1
#                     CaseName    VarName2  VarWeight2
#                     CaseName    VarNameN  VarWeightN
#           ENDATA
# ==============================================================================

mutable struct TempColumn
    lower::Float64
    upper::Float64
    is_int::Bool
    TempColumn() = new(0.0, Inf, false)
end

mutable struct TempRow
    lower::Float64
    upper::Float64
    sense::String
    terms::Dict{String, Float64}
    TempRow() = new(-Inf, Inf, "E", Dict{String, Float64}())
end

mutable struct TempMPSModel
    name::String
    obj_name::String
    columns::Dict{String, TempColumn}
    rows::Dict{String, TempRow}
    intorg_flag::Bool  # A flag used to parse COLUMNS section.
    TempMPSModel() = new("", "", Dict{String, TempColumn}(),
        Dict{String, TempRow}(), false)
end

const HEADERS = ("ROWS", "COLUMNS", "RHS", "RANGES", "BOUNDS", "SOS", "ENDATA")

function MOI.read_from_file(model::InnerModel, io::IO)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end
    data = TempMPSModel()
    header = "NAME"
    multi_objectives = String[]
    while !eof(io) && header != "ENDATA"
        line = strip(readline(io))
        if line == "" || startswith(line, "*")
            # Skip blank lines and comments.
            continue
        end
        if uppercase(string(line)) in HEADERS
            header = uppercase(string(line))
            continue
        end
        # TODO: split into hard fields based on column indices.
        items = String.(split(line, " ", keepempty = false))
        if header == "NAME"
            # A special case. This only happens at the start.
            parse_name_line(data, items)
        elseif header == "ROWS"
            multi_obj = parse_rows_line(data, items)
            multi_obj !== nothing && push!(multi_objectives, multi_obj)
        elseif header == "COLUMNS"
            parse_columns_line(data, items, multi_objectives)
        elseif header == "RHS"
            parse_rhs_line(data, items)
        elseif header == "RANGES"
            parse_ranges_line(data, items)
        elseif header == "BOUNDS"
            parse_bounds_line(data, items)
        end
    end
    copy_to(model, data)
    return
end

function bounds_to_set(lower, upper)
    if -Inf < lower < upper < Inf
        return MOI.Interval(lower, upper)
    elseif -Inf < lower && upper == Inf
        return MOI.GreaterThan(lower)
    elseif -Inf == lower && upper < Inf
        return MOI.LessThan(upper)
    elseif lower == upper
        return MOI.EqualTo(upper)
    end
    return nothing
end

function copy_to(model::InnerModel, temp::TempMPSModel)
    MOI.set(model, MOI.Name(), temp.name)
    variable_map = Dict{String, MOI.VariableIndex}()
    # Add variables.
    for (name, column) in temp.columns
        x = MOI.add_variable(model)
        variable_map[name] = x
        MOI.set(model, MOI.VariableName(), x, name)
        set = bounds_to_set(column.lower, column.upper)
        if set === nothing && column.is_int
            # TODO: some solvers may interpret this as binary.
            MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
        elseif set !== nothing && column.is_int
            if set == MOI.Interval(0.0, 1.0)
                MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
            else
                MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
                MOI.add_constraint(model, MOI.SingleVariable(x), set)
            end
        elseif set !== nothing
            MOI.add_constraint(model, MOI.SingleVariable(x), set)
        end
    end
    # Add linear constraints.
    for (c_name, row) in temp.rows
        if c_name == temp.obj_name
            # Set objective.
            MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
            obj_func = if length(row.terms) == 1 &&
                    first(row.terms).second == 1.0
                MOI.SingleVariable(variable_map[first(row.terms).first])
            else
                MOI.ScalarAffineFunction([
                    MOI.ScalarAffineTerm(coef, variable_map[v_name])
                        for (v_name, coef) in row.terms],
                0.0)
            end
            MOI.set(model, MOI.ObjectiveFunction{typeof(obj_func)}(), obj_func)
        else
            constraint_function = MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(coef, variable_map[v_name])
                    for (v_name, coef) in row.terms],
                0.0)
            set = bounds_to_set(row.lower, row.upper)
            if set !== nothing
                c = MOI.add_constraint(model, constraint_function, set)
                MOI.set(model, MOI.ConstraintName(), c, c_name)
            else
                error("Expected a non-empty set for $(c_name). Got row=$(row)")
            end
        end
    end
    return
end

# ==============================================================================
#   NAME
# ==============================================================================

function parse_name_line(data::TempMPSModel, items::Vector{String})
    if !(1 <= length(items) <= 2) || uppercase(items[1]) != "NAME"
        error("Malformed NAME line: $(join(items, " "))")
    end
    if length(items) == 2
        data.name = items[2]
    elseif length(items) == 1
        data.name = ""
    end
    return
end

# ==============================================================================
#   ROWS
# ==============================================================================

function parse_rows_line(data::TempMPSModel, items::Vector{String})
    if length(items) != 2
        error("Malformed ROWS line: $(join(items, " "))")
    end
    sense, name = items
    if haskey(data.rows, name)
        error("Duplicate row encountered: $(line).")
    elseif sense != "N" && sense != "L" && sense != "G" && sense != "E"
        error("Invalid row sense: $(join(items, " "))")
    end
    row = TempRow()
    row.sense = sense
    if sense == "N"
        if data.obj_name != ""
            # Detected a duplicate objective. Skip it.
            return name
        end
        data.obj_name = name
    end
    # Add some default bounds for the constraints.
    if sense == "G"
        row.lower = 0.0
    elseif sense == "L"
        row.upper = 0.0
    elseif sense == "E"
        row.lower = 0.0
        row.upper = 0.0
    end
    data.rows[name] = row
    return
end

# ==============================================================================
#   COLUMNS
# ==============================================================================

function parse_single_coefficient(data, row_name, column_name, value)
    row = get(data.rows, row_name, nothing)
    if row === nothing
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    if haskey(row.terms, column_name)
        row.terms[column_name] += parse(Float64, value)
    else
        row.terms[column_name] = parse(Float64, value)
    end
    return
end

function parse_columns_line(data::TempMPSModel, items::Vector{String}, multi_objectives::Vector{String})
    if length(items) == 3
        # [column name] [row name] [value]
        column_name, row_name, value = items
        if uppercase(row_name) == "'MARKER'" && uppercase(value) == "'INTORG'"
            data.intorg_flag = true
            return
        elseif uppercase(row_name) == "'MARKER'" &&
                uppercase(value) == "'INTEND'"
            data.intorg_flag = false
            return
        elseif row_name in multi_objectives
            return
        end
        if !haskey(data.columns, column_name)
            data.columns[column_name] = TempColumn()
        end
        parse_single_coefficient(data, row_name, column_name, value)
        if data.columns[column_name].is_int && !data.intorg_flag
            error("Variable $(column_name) appeared in COLUMNS outside an" *
                  " `INT` marker after already being declared as integer.")
        end
        data.columns[column_name].is_int = data.intorg_flag
    elseif length(items) == 5
        # [column name] [row name] [value] [row name 2] [value 2]
        column_name, row_name_1, value_1, row_name_2, value_2 = items
        if !haskey(data.columns, column_name)
            data.columns[column_name] = TempColumn()
        end
        if row_name_1 in multi_objectives || row_name_2 in multi_objectives
            return
        end
        parse_single_coefficient(data, row_name_1, column_name, value_1)
        parse_single_coefficient(data, row_name_2, column_name, value_2)
        if data.columns[column_name].is_int && !data.intorg_flag
            error("Variable $(column_name) appeared in COLUMNS outside an" *
                  " `INT` marker after already being declared as integer.")
        end
        data.columns[column_name].is_int = data.intorg_flag
    else
        error("Malformed COLUMNS line: $(join(items, " "))")
    end
    return
end

# ==============================================================================
#   RHS
# ==============================================================================

function parse_single_rhs(data, row_name, value)
    if !haskey(data.rows, row_name)
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    value = parse(Float64, value)
    sense = data.rows[row_name].sense
    if sense == "E"
        data.rows[row_name].upper = value
        data.rows[row_name].lower = value
    elseif sense == "G"
        data.rows[row_name].lower = value
    elseif sense == "L"
        data.rows[row_name].upper = value
    elseif sense == "N"
        error("Cannot have RHS for objective: $(join(items, " "))")
    end
    return
end

# TODO: handle multiple RHS vectors.
function parse_rhs_line(data::TempMPSModel, items::Vector{String})
    if length(items) == 3
        # [rhs name] [row name] [value]
        rhs_name, row_name, value = items
        parse_single_rhs(data, row_name, value)
    elseif length(items) == 5
        # [rhs name] [row name 1] [value 1] [row name 2] [value 2]
        rhs_name, row_name_1, value_1, row_name_2, value_2 = items
        parse_single_rhs(data, row_name_1, value_1)
        parse_single_rhs(data, row_name_2, value_2)
    else
        error("Malformed RHS line: $(join(items, " "))")
    end
    return
end

# ==============================================================================
#  RANGES
#
# Here is how RANGE information is encoded. (We repeat this comment because it
# is so non-trivial.)
#
#     Row type | Range value |  lower bound  |  upper bound
#     ------------------------------------------------------
#         G    |     +/-     |     rhs       | rhs + |range|
#         L    |     +/-     | rhs - |range| |     rhs
#         E    |      +      |     rhs       | rhs + range
#         E    |      -      | rhs + range   |     rhs
# ==============================================================================

function parse_single_range(data, row_name, value)
    if !haskey(data.rows, row_name)
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    value = parse(Float64, value)
    row = data.rows[row_name]
    if row.sense == "G"
        row.upper = row.lower + abs(value)
    elseif row.sense == "L"
        row.lower = row.upper - abs(value)
    elseif row.sense == "E"
        if value > 0.0
            row.upper = row.upper + value
        else
            row.lower = row.lower + value
        end
    end
    return
end

# TODO: handle multiple RANGES vectors.
function parse_ranges_line(data::TempMPSModel, items::Vector{String})
    if length(items) == 3
        # [rhs name] [row name] [value]
        rhs_name, row_name, value = items
        parse_single_range(data, row_name, value)
    elseif length(items) == 5
        # [rhs name] [row name] [value] [row name 2] [value 2]
        rhs_name, row_name, value, row_name_2, value_2 = items
        parse_single_range(data, row_name, value)
        parse_single_range(data, row_name_2, value_2)
    else
        error("Malformed RANGES line: $(join(items, " "))")
    end
    return
end

# ==============================================================================
#   BOUNDS
# ==============================================================================

function parse_bounds_line(data::TempMPSModel, items::Vector{String})
    if length(items) == 3
        bound_type, bound_name, column_name = items
        if !haskey(data.columns, column_name)
            error("Column name $(column_name) not found.")
        end
        column = data.columns[column_name]
        if bound_type == "PL"
            column.upper = Inf
        elseif bound_type == "MI"
            column.lower = -Inf
        elseif bound_type == "FR"
            column.lower = -Inf
            column.upper = Inf
        elseif bound_type == "BV"
            column.lower = 0.0
            column.upper = 1.0
            column.is_int = true
        else
            error("Invalid bound type $(bound_type): $(join(items, " "))")
        end
    elseif length(items) == 4
        bound_type, bound_name, column_name, value = items
        if !haskey(data.columns, column_name)
            error("Column name $(column_name) not found.")
        end
        column = data.columns[column_name]
        value = parse(Float64, value)
        if bound_type == "FX"
            column.lower = column.upper = value
        elseif bound_type == "UP"
            column.upper = value
        elseif bound_type == "LO"
            column.lower = value
        elseif bound_type == "LI"
            column.lower = value
            column.is_int = true
        elseif bound_type == "UI"
            column.upper = value
            column.is_int = true
        elseif bound_type == "FR"
            # So even though FR bounds should be of the form:
            #  FR BOUND1    VARNAME
            # there are cases in MIPLIB2017 (e.g., leo1 and leo2) like so:
            #  FR BOUND1    C0000001       .000000
            # In these situations, just ignore the value.
            column.lower = -Inf
            column.upper = Inf
        elseif bound_type == "BV"
            # A similar situation happens with BV bounds in leo1 and leo2.
            column.lower = 0.0
            column.upper = 1.0
            column.is_int = true
        else
            error("Invalid bound type $(bound_type): $(join(items, " "))")
        end
    else
        error("Malformed BOUNDS line: $(join(items, " "))")
    end
    return
end

end


text:
using CSFML.LibCSFML
using ModernGL
using CSyntax
using CSyntax.CSwitch

# shader sources
const vert_source = """
    #version 150 core
    in vec2 position;
    in vec3 color;
    out vec3 Color;
    void main()
    {
        Color = color;
        gl_Position = vec4(position, 0.0, 1.0);
    }"""

const frag_source = """
    #version 150 core
    in vec3 Color;
    out vec4 outColor;
    void main()
    {
        outColor = vec4(Color, 1.0);
    }"""

# SFML settings
settings = sfContextSettings(24, 8, 0, 3, 3, sfContextCore, true)
window = @c sfWindow_create(sfVideoMode(800, 600, 32), "OpenGL", sfTitlebar | sfClose, &settings)

# create Vertex Array Object
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)

# create a Vertex Buffer Object and copy the vertex data to it
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)

vertices = GLfloat[
     0.0,  0.5, 1.0, 0.0, 0.0,
     0.5, -0.5, 0.0, 1.0, 0.0,
    -0.5, -0.5, 0.0, 0.0, 1.0]

glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW)

# create and compile the vertex shader
vert_shader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vert_shader, 1, Ptr{GLchar}[pointer(vert_source)], C_NULL)
glCompileShader(vert_shader)

# create and compile the fragment shader
frag_shader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(frag_shader, 1, Ptr{GLchar}[pointer(frag_source)], C_NULL)
glCompileShader(frag_shader)

# link the vertex and fragment shader into a shader program
shader_prog = glCreateProgram()
glAttachShader(shader_prog, vert_shader)
glAttachShader(shader_prog, frag_shader)
glBindFragDataLocation(shader_prog, 0, "outColor")
glLinkProgram(shader_prog)
glUseProgram(shader_prog)

# specify the layout of the vertex data
pos_attrib = glGetAttribLocation(shader_prog, "position")
glEnableVertexAttribArray(pos_attrib)
glVertexAttribPointer(pos_attrib, 2, GL_FLOAT, GL_FALSE, 5 * sizeof(GLfloat), C_NULL)

color_attrib = glGetAttribLocation(shader_prog, "color")
glEnableVertexAttribArray(color_attrib);
glVertexAttribPointer(color_attrib, 3, GL_FLOAT, GL_FALSE, 5 * sizeof(GLfloat), Ptr{Cvoid}(2 * sizeof(GLfloat)))

window_event_ref = Ref(sfEvent(sfEvtClosed))
running = true
while running
    global running; global window_event_ref;
    while Bool(sfWindow_pollEvent(window, window_event_ref))
        @cswitch window_event_ref[].type begin
            @case sfEvtClosed
                running = false
                break
        end
    end

    # clear the screen to black
    glClearColor(0.0, 0.0, 0.0, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)

    # draw a triangle from the 3 vertices
    glDrawArrays(GL_TRIANGLES, 0, 3)

    # swap buffers
    sfWindow_display(window)
end

glDeleteProgram(shader_prog)
glDeleteShader(frag_shader)
glDeleteShader(vert_shader)

@c glDeleteBuffers(1, &vbo)

@c glDeleteVertexArrays(1, &vao)

sfWindow_close(window)
sfWindow_destroy(window)

ext: jl

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
export ODESystem, ODEFunction

using Base: RefValue

isintermediate(eq::Equation) = !(isa(eq.lhs, Operation) && isa(eq.lhs.op, Differential))

function flatten_differential(O::Operation)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(O.args[1]) || return (O.args[1], O.op.x, 1)
    (x, t, order) = flatten_differential(O.args[1])
    isequal(t, O.op.x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(O.op.x)"))
    return (x, t, order + 1)
end

struct DiffEq  # dⁿx/dtⁿ = rhs
    x::Variable
    n::Int
    rhs::Expression
end
function to_diffeq(eq::Equation)
    isintermediate(eq) && throw(ArgumentError("intermediate equation received"))
    (x, t, n) = flatten_differential(eq.lhs)
    (isa(t, Operation) && isa(t.op, Variable) && isempty(t.args)) ||
        throw(ArgumentError("invalid independent variable $t"))
    (isa(x, Operation) && isa(x.op, Variable) && length(x.args) == 1 && isequal(first(x.args), t)) ||
        throw(ArgumentError("invalid dependent variable $x"))
    return t.op, DiffEq(x.op, n, eq.rhs)
end
Base.:(==)(a::DiffEq, b::DiffEq) = isequal((a.x, a.n, a.rhs), (b.x, b.n, b.rhs))

"""
$(TYPEDEF)

A system of ordinary differential equations.

# Fields
* `eqs` - The ODEs defining the system.

# Examples

```
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = ODESystem(eqs)
```
"""
struct ODESystem <: AbstractSystem
    """The ODEs defining the system."""
    eqs::Vector{DiffEq}
    """Independent variable."""
    iv::Variable
    """Dependent (state) variables."""
    dvs::Vector{Variable}
    """Parameter variables."""
    ps::Vector{Variable}
    """
    Jacobian matrix. Note: this field will not be defined until
    [`calculate_jacobian`](@ref) is called on the system.
    """
    jac::RefValue{Matrix{Expression}}
    """
    Wfact matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact::RefValue{Matrix{Expression}}
    """
    Wfact_t matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact_t::RefValue{Matrix{Expression}}
end

function ODESystem(eqs)
    reformatted = to_diffeq.(eqs)

    ivs = unique(r[1] for r ∈ reformatted)
    length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
    iv = first(ivs)

    deqs = [r[2] for r ∈ reformatted]

    dvs = unique(deq.x for deq ∈ deqs)
    ps = filter(vars(deq.rhs for deq ∈ deqs)) do x
        x.known & !isequal(x, iv)
    end |> collect

    ODESystem(deqs, iv, dvs, ps)
end

function ODESystem(deqs::AbstractVector{DiffEq}, iv, dvs, ps)
    jac = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
    ODESystem(deqs, iv, dvs, ps, jac, Wfact, Wfact_t)
end

function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps)
    _dvs = [deq.op for deq ∈ dvs]
    _iv = iv.op
    _ps = [p.op for p ∈ ps]
    ODESystem(getindex.(to_diffeq.(deqs),2), _iv, _dvs, _ps)
end

function _eq_unordered(a, b)
    length(a) === length(b) || return false
    n = length(a)
    idxs = Set(1:n)
    for x ∈ a
        idx = findfirst(isequal(x), b)
        idx === nothing && return false
        idx ∈ idxs      || return false
        delete!(idxs, idx)
    end
    return true
end
Base.:(==)(sys1::ODESystem, sys2::ODESystem) =
    _eq_unordered(sys1.eqs, sys2.eqs) && isequal(sys1.iv, sys2.iv) &&
    _eq_unordered(sys1.dvs, sys2.dvs) && _eq_unordered(sys1.ps, sys2.ps)
# NOTE: equality does not check cached Jacobian

independent_variables(sys::ODESystem) = Set{Variable}([sys.iv])
dependent_variables(sys::ODESystem) = Set{Variable}(sys.dvs)
parameters(sys::ODESystem) = Set{Variable}(sys.ps)

function calculate_jacobian(sys::ODESystem)
    isempty(sys.jac[]) || return sys.jac[]  # use cached Jacobian, if possible
    rhs = [eq.rhs for eq ∈ sys.eqs]

    iv = sys.iv()
    dvs = [dv(iv) for dv ∈ sys.dvs]

    jac = expand_derivatives.(calculate_jacobian(rhs, dvs))
    sys.jac[] = jac  # cache Jacobian
    return jac
end

struct ODEToExpr
    sys::ODESystem
end
function (f::ODEToExpr)(O::Operation)
    if isa(O.op, Variable)
        isequal(O.op, f.sys.iv) && return O.op.name  # independent variable
        O.op ∈ f.sys.dvs        && return O.op.name  # dependent variables
        isempty(O.args)         && return O.op.name  # 0-ary parameters
        return build_expr(:call, Any[O.op.name; f.(O.args)])
    end
    return build_expr(:call, Any[Symbol(O.op); f.(O.args)])
end
(f::ODEToExpr)(x) = convert(Expr, x)

function generate_jacobian(sys::ODESystem, dvs = sys.dvs, ps = sys.ps, expression = Val{true})
    jac = calculate_jacobian(sys)
    return build_function(jac, dvs, ps, (sys.iv.name,), ODEToExpr(sys), expression)
end

function generate_function(sys::ODESystem, dvs = sys.dvs, ps = sys.ps, expression = Val{true})
    rhss = [deq.rhs for deq ∈ sys.eqs]
    dvs′ = [clean(dv) for dv ∈ dvs]
    ps′ = [clean(p) for p ∈ ps]
    return build_function(rhss, dvs′, ps′, (sys.iv.name,), ODEToExpr(sys), expression)
end

function calculate_factorized_W(sys::ODESystem, simplify=true)
    isempty(sys.Wfact[]) || return (sys.Wfact[],sys.Wfact_t[])

    jac = calculate_jacobian(sys)
    gam = Variable(:gam; known = true)()

    W = - LinearAlgebra.I + gam*jac
    Wfact = lu(W, Val(false), check=false).factors

    if simplify
        Wfact = simplify_constants.(Wfact)
    end

    W_t = - LinearAlgebra.I/gam + jac
    Wfact_t = lu(W_t, Val(false), check=false).factors
    if simplify
        Wfact_t = simplify_constants.(Wfact_t)
    end
    sys.Wfact[] = Wfact
    sys.Wfact_t[] = Wfact_t

    (Wfact,Wfact_t)
end

function generate_factorized_W(sys::ODESystem, vs = sys.dvs, ps = sys.ps, simplify=true, expression = Val{true})
    (Wfact,Wfact_t) = calculate_factorized_W(sys,simplify)
    siz = size(Wfact)
    constructor = :(x -> begin
                        A = SMatrix{$siz...}(x)
                        StaticArrays.LU(LowerTriangular( SMatrix{$siz...}(UnitLowerTriangular(A)) ), UpperTriangular(A), SVector(ntuple(n->n, max($siz...))))
                    end)

    Wfact_func   = build_function(Wfact  , vs, ps, (:gam,:t), ODEToExpr(sys), expression;constructor=constructor)
    Wfact_t_func = build_function(Wfact_t, vs, ps, (:gam,:t), ODEToExpr(sys), expression;constructor=constructor)

    return (Wfact_func, Wfact_t_func)
end

"""
$(SIGNATURES)

Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.ODEFunction{iip}(sys::ODESystem, dvs, ps;
                                     version = nothing,
                                     jac = false, Wfact = false) where {iip}
    f_oop,f_iip = generate_function(sys, dvs, ps, Val{false})

    f(u,p,t) = f_oop(u,p,t)
    f(du,u,p,t) = f_iip(du,u,p,t)

    if jac
        jac_oop,jac_iip = generate_jacobian(sys, dvs, ps, Val{false})
        _jac(u,p,t) = jac_oop(u,p,t)
        _jac(J,u,p,t) = jac_iip(J,u,p,t)
    else
        _jac = nothing
    end

    if Wfact
        tmp_Wfact,tmp_Wfact_t = generate_factorized_W(sys, dvs, ps, Val{false})
        Wfact_oop, Wfact_iip = tmp_Wfact
        Wfact_oop_t, Wfact_iip_t = tmp_Wfact_t
        _Wfact(u,p,t) = Wfact_oop(u,p,t)
        _Wfact(W,u,p,t) = Wfact_iip(W,u,p,t)
        _Wfact_t(u,p,t) = Wfact_oop_t(u,p,t)
        _Wfact_t(W,u,p,t) = Wfact_iip_t(W,u,p,t)
    else
        _Wfact,_Wfact_t = nothing,nothing
    end

    ODEFunction{iip}(f,jac=_jac,
                      Wfact = _Wfact,
                      Wfact_t = _Wfact_t)
end

function DiffEqBase.ODEFunction(sys::ODESystem, args...; kwargs...)
    ODEFunction{true}(sys, args...; kwargs...)
end

"""
$(SIGNATURES)

Generate `ODESystem`, dependent variables, and parameters from an `ODEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
    @parameters t
    vars = [Variable(:x, i)(t) for i in eachindex(prob.u0)]
    params = [Variable(:α,i; known = true)() for i in eachindex(prob.p)]
    @derivatives D'~t

    rhs = [D(var) for var in vars]

    if DiffEqBase.isinplace(prob)
        lhs = similar(vars, Any)
        prob.f(lhs, vars, params, t)
    else
        lhs = prob.f(vars, params, t)
    end

    eqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)
    de = ODESystem(eqs)

    de, vars, params
end


text:
# Generalized ADMM function based on Boyd et al's Matlab scripts:
#   http://web.stanford.edu/~boyd/papers/admm/

abstract ADMMState
abstract NonconvexADMMState <: ADMMState

function admm_loop(s::ADMMState;
    MaxIterations = 1000,
    Quiet = true,
    ABSTOL::Float64 = 1e-4,
    RELTOL::Float64 = 1e-2)

    local r_norm::Float64, s_norm::Float64, eps_pri::Float64, eps_dual::Float64
    local best_objval::Float64 = Inf # Make sure its updated in 1st iteration
    local best_solution

    for k in 1:MaxIterations
        zold = s.z
        update_xzu!(s, k)

        # Check convergence
        r_norm  = norm(s.x - s.z)
        s_norm  = norm(-s.rho*(s.z - zold))
        eps_pri = sqrt(s.n) * ABSTOL + RELTOL * max(norm(s.x), norm(-s.z))
        eps_dual= sqrt(s.n) * ABSTOL + RELTOL * norm(s.rho*s.u)
        objval = objective(s)

        if !Quiet
            println(@sprintf("%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f", k,
                 r_norm, eps_pri, s_norm, eps_dual, objval))
        end

        if objval < best_objval
            best_objval = objval
            best_solution = copy(solution(s))
        end

        if (r_norm < eps_pri) && (s_norm < eps_dual)
            break
        end
    end

    # Return the solution with best objective value if this is non-convex ADMM.
    if typeof(s) <: NonconvexADMMState
        return best_solution
    else
        return solution(s)
    end
end


text:
module aggregation_strategy

export AggregationType, parseAggStrategy, update_score_matrix!, aggregate, Aggregation

using util
using Statistics

@enum AggregationType MAX AVG_QUERY AVG_CAND AVG_SHORT AVG_LONG AVG_QUERY_CAND

function parseAggStrategy(methodName)
    if methodName == "max"
        return MAX::AggregationType
    elseif methodName == "avg_query"
        return AVG_QUERY::AggregationType
    elseif methodName == "avg_cand"
        return AVG_CAND::AggregationType
    elseif methodName == "avg_short"
        return AVG_SHORT::AggregationType
    elseif methodName == "avg_long"
        return AVG_LONG::AggregationType
    elseif methodName == "avg_query_cand"
        return AVG_QUERY_CAND::AggregationType
    end
end

struct Aggregation
    aggType::AggregationType
    query_max_vec::Vector{Float64}
    cand_max_vec::Vector{Float64}
end

Aggregation(aggType::AggregationType, q_len, c_len) = Aggregation(aggType, fill(typemin(Float64), q_len),  fill(typemin(Float64), c_len))

function update_score_matrix!(agg::Aggregation, query_idx, cand_idx, score)
    if agg.query_max_vec[query_idx] < score
        agg.query_max_vec[query_idx] = score
    end

    if agg.cand_max_vec[cand_idx] < score
        agg.cand_max_vec[cand_idx] = score
    end
end

function aggregate(agg::Aggregation)::Float64
    if agg.aggType == MAX::AggregationType
        largest = typemin(Float64)

        for v in agg.query_max_vec
            if v > largest
                largest = v
            end
        end

        for v in agg.cand_max_vec
            if v > largest
                largest = v
            end
        end

        return largest
    elseif agg.aggType == AVG_QUERY::AggregationType
        return mean(agg.query_max_vec)
    elseif agg.aggType == AVG_CAND::AggregationType
        return mean(agg.cand_max_vec)
    elseif agg.aggType == AVG_SHORT::AggregationType
        querySize = length(agg.query_max_vec)
        candSize = length(agg.cand_max_vec)

        return querySize < candSize  ? mean(agg.query_max_vec) : mean(agg.cand_max_vec)
    elseif agg.aggType == AVG_LONG::AggregationType
        querySize = length(agg.query_max_vec)
        candSize = length(agg.cand_max_vec)

        return querySize > candSize  ? mean(agg.query_max_vec) : mean(agg.cand_max_vec)
    elseif agg.aggType == AVG_QUERY_CAND::AggregationType
        return (mean(agg.query_max_vec) + mean(agg.cand_max_vec)) / 2.0
    end

    return typemax(Float64)
end

end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
using MCAnalyzer
using Test

import MCAnalyzer: code_native, code_llvm

f() = mark_start()
g() = mark_end()

let buf = IOBuffer()
    code_native(buf, f, Tuple{})
    asm = String(take!(buf))
    @test occursin("movl\t\$111, %ebx", asm)
end

let buf = IOBuffer()
    code_native(buf, g, Tuple{})
    asm = String(take!(buf))
    @test occursin("movl\t\$222, %ebx", asm)
end

function mysum(A)
    acc = zero(eltype(A))
    for i in eachindex(A)
        mark_start()
        @inbounds acc += A[i]
    end
    mark_end()
    return acc
end

let buf = IOBuffer()
    code_native(buf, mysum, Tuple{Vector{Float64}})
    asm = String(take!(buf))
    @test occursin("movl\t\$111, %ebx", asm)
    @test occursin("movl\t\$222, %ebx", asm)
end

if Sys.which("iaca") !== nothing
    @test_nowarn analyze(mysum, (Vector{Float64},))
end

using LinearAlgebra

function mynorm(out, X)
    for i in 1:length(X)
        @inbounds out[i] = X[i] / norm(X)
    end
    out
end

fname, io = mktemp()
@test_nowarn code_llvm(io, mynorm, (Vector{Float64}, Vector{Float64}))
close(io)
rm(fname)


text:
export PCGPChromo, node_genes, get_positions

# PCGP with mutated positions and constant scalar inputs and evolved params

struct PCGPChromo <: Chromosome
    genes::Array{Float64}
    nodes::Array{CGPNode}
    outputs::Array{Int64}
    nin::Int64
    nout::Int64
end

function PCGPChromo(genes::Array{Float64}, nin::Int64, nout::Int64)::PCGPChromo
    num_nodes = Int64(ceil((length(genes)-nin-nout)/5))
    nodes = Array{CGPNode}(nin+num_nodes)
    rgenes = reshape(genes[(nin+nout+1):end], (5, num_nodes))
    rgenes = sortcols(rgenes)'
    genes[1:nin] = sort(genes[1:nin], rev=true)
    # genes[nin+(1:nout)] = sort(genes[nin+(1:nout)])
    positions = [Config.input_start .* genes[1:nin]; rgenes[:, 1]]
    fc = deepcopy(hcat(zeros(2, nin), [rgenes[:, 2]'; rgenes[:, 3]']))
    e = (1.0 .- positions) .* Config.recurrency .+ positions
    for j in 1:size(fc)[1]
        fc[j, :] = (fc[j, :] .* (e .- Config.input_start)) .+ Config.input_start
    end
    if Config.recurrency == 0
        for i in (nin+1):length(positions)
            for j in 1:size(fc)[1]
                fc[j, i] = max(fc[j,i], positions[i-1])
            end
        end
    end
    connections = snap(fc, positions)
    outputs = snap(genes[nin+(1:nout)], positions)
    functions = Array{Function}(nin+num_nodes)
    functions[1:nin] = Config.f_input
    functions[(nin+1):end] = map(i->Config.index_in(Config.functions, i), rgenes[:, 4])
    params = [zeros(nin); 2.0*rgenes[:, 5]-1.0]
    active = find_active(nin, outputs, connections)
    for i in 1:(nin+num_nodes)
        nodes[i] = CGPNode(connections[:, i], functions[i], active[i], params[i])
    end
    PCGPChromo([genes[1:(nin+nout)]; rgenes'[:]], nodes, outputs, nin, nout)
end

function PCGPChromo(nin::Int64, nout::Int64)::PCGPChromo
    n_nodes = Config.static_node_size
    if Config.bloat()
        n_nodes = Config.starting_nodes
    end
    PCGPChromo(rand(nin+nout+5*n_nodes), nin, nout)
end

function PCGPChromo(c::PCGPChromo)::PCGPChromo
    gene_mutate(c)
end

function node_genes(c::PCGPChromo)
    5
end

function get_positions(c::PCGPChromo)
    [Config.input_start .* c.genes[1:c.nin]; c.genes[(c.nin+c.nout+1):5:end]]
end


text:
create_filter(x::AbstractArray{<:Real}) = "IN ($(join(x, ", ")))"
create_filter(x::AbstractArray) = "IN ('$(join(x, "', '"))')"
create_filter(x::Missing) = "IS NOT NULL"
create_filter(x::Real) = "= $x"
create_filter(x::AbstractString) = "= '$x'"

function create_filter(
    filters::Dict{String, <:Any},
    fil = ""
)
    for (key, data) in filters
        if length(fil) > 0
            fil *= " AND "
        else
            fil *= " WHERE "
        end
        fil *= "$key $(create_filter(data))"
    end
    return fil
end


text:
using WAV
using SpeechAugment: dropWav

# read a wav file
wav, fs = wavread("xxx.wav")
# simulate packet loss effects by dropping frames randomly.
# recommended ratio is ∈ [0.0, 0.05]
y = dropWav(wav, fs; ratio=0.1);

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
using Documenter, Pseudospectra, Plots

makedocs(modules = [Pseudospectra],
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == true),
         sitename = "Pseudospectra.jl",
         pages = [
             "Home" => "index.md",
             "Usage" => "usage.md",
             "Library" => Any[
                 "Public" => "lib/public.md",
                 "Internals" => "lib/internals.md",
                 "Example matrix generators" => "lib/demos.md"
             ]
         ]
         )

deploydocs(repo = "github.com/RalphAS/Pseudospectra.jl.git")


text:
"""
discretePulse(pulse, gridPot ; rpulse = false)

A pulse is discretized given by a fixed grid-array of values from a general
reference pulse.
The output is an array with the _index positions_ of the corresponding values
from the grid-array. Also, is posible to obtain a _new pulse_ with the
corresponding values.
"""
function discretePulse(pulse, gridPot ; rpulse = false)
    ΔPulse = zeros(length(pulse))
    indxΔts = zeros(Int64,length(pulse))
    for (indx, fᵢ) in enumerate(pulse)
        gaIndx = argmin((fᵢ .- gridPot).^2)
        ΔPulse[indx] = gridPot[gaIndx]
        indxΔts[indx] = gaIndx
    end
    rpulse == true ? (indxΔts, ΔPulse) : indxΔts
end

"""
getΔtsPulse(tlims, pulse, gridPot)

Computes new Δts from a discretized pulse given the general reference pulse.
The output is a new set of _Δts_ and the corresponding _indices_ for the values
in the reference pulse.
"""
function getΔtsPulse(tlims, pulse, gridPot)
    indxΔts = discretePulse(pulse, gridPot)
    tᵢ = indxΔts[1]
    setΔtij = []
    indxi = 1
    uindxΔts = [tᵢ]
    for (index, tⱼ) in enumerate(indxΔts)
        if tⱼ != tᵢ
            δtij = (index -indxi)*dt
            push!(setΔtij, δtij)
            tᵢ = tⱼ
            push!(uindxΔts, tᵢ)
            indxi = index
        elseif index==length(indxΔts) && tⱼ == tᵢ
            δtij = (index - indxi + 1)*dt
            push!(setΔtij, δtij)
        end
    end
    setΔtij, uindxΔts
end

"""
diagLightMatter(filename, gridPot, Eα, Vαβ; dimBase=600)

Diagonalization of a given discretized range of values in the vector potential.
"""
function diagLightMatter(filename, gridPot, Eα, Vαβ; dimBase=600)
    sizeGridPot = size(gridPot)[1]
    filediag = h5open(filename * ".h5", "w")
    dset_psi = d_create(filediag,"psi", datatype(Complex{Float64}), dataspace(dimBase, dimBase, sizeGridPot))
    dset_invpsi = d_create(filediag,"invpsi", datatype(Complex{Float64}), dataspace(dimBase, dimBase, sizeGridPot))
    dset_energ = d_create(filediag,"energies", datatype(Complex{Float64}), dataspace(dimBase, sizeGridPot))
    p = Progress(sizeGridPot, dt=0.2, desc="diagonalizing...", color=:blue)
    for (i, aᵢ) in enumerate(gridPot)
        H̃αβ = Diagonal(Eα) .- aᵢ*Vαβ
        Ẽα, ψ̃α = eigen!(H̃αβ)
        dset_invpsi[:, :, i] = inv(ψ̃α)
        dset_psi[:, :, i] = ψ̃α
        dset_energ[:, i] = Ẽα
        szfile = filesize(filename * ".h5")
        ProgressMeter.next!(p; showvalues = [(:Filesize, szprintp(szfile))])
        IJulia.clear_output(true) # remove this if you run it in the terminal
    end
    flush(filediag)
    close(filediag)
end
function diagLightMatterS(gridPot, Eα, Vαβ; dimBase=600)
    sizeGridPot = size(gridPot)[1]
    dset_psi = zeros(Complex{Float64}, dimBase, dimBase, sizeGridPot)
    dset_invpsi = zeros(Complex{Float64},dimBase, dimBase, sizeGridPot)
    dset_energ = zeros(dimBase,sizeGridPot)
    #p = Progress(sizeGridPot, dt=0.2, desc="diagonalizing...", color=:blue)
    for i in 1:sizeGridPot
        H̃αβ = Diagonal(Eα) .- gridPot[i]*Vαβ
        Ẽα, ψ̃α = eigen!(H̃αβ)
        dset_invpsi[:, :, i] = inv(ψ̃α)
        dset_psi[:, :, i] = ψ̃α
        dset_energ[:, i] = real.(Ẽα)
       # ProgressMeter.next!(p)
    end
    dset_psi, dset_invpsi, dset_energ
end

"""
timeEvolution(Δtij, posVecPot, ncut = 600)

Time evolution for a given state ψt₀ (provided by default) is computed,
for varying times and strength's potentials. Before using this function
please run _getΔtsPulse_.
"""
function timeEvolution(Δtij, posVecPot, fψ̃α, finvψ̃α, fEα, ncut = 600)
    ψt₀ = zeros(Complex{Float64}, ncut)
    ψt₀[1] = 1
    for (index, pᵢ) in enumerate(posVecPot)
        sψ̃α = @view fψ̃α[:,:,pᵢ]
        sinvψ̃α = @view finvψ̃α[:,:,pᵢ]
        sEα = @view fEα[:,pᵢ]
        sΔt = @view Δtij[index]
        # time evolution
        ψtᵢ = sψ̃α * (exp.(-1im*real(sEα) .* sΔt) .* (sinvψ̃α * ψt₀))
        # be careful with the order in operations and the change from * to .*
        ψt₀ .= ψtᵢ
    end
    ψt₀
end

"""
timeEvolutionCN(pulse_ex, Vαβ, Eα, dt, ncut = 600)

Applying the Crank-nicolson method th time evolution for a given state ψt₀
(provided by default) is computed,for varying times and strength's potentials.
"""
function timeEvolutionCN(pulse, Vαβ, Eα, dt, ncut = 600)
    ψt₀ = zeros(Complex{Float64}, ncut)
    ψt₀[1] = 1
    p = Progress(size(pulse)[1], dt=0.01, desc="propagating...", color=:blue)
    for aᵢ in pulse
        H̃αβ = Diagonal(Eα) .- aᵢ*Vαβ
        Uf = -1im .* dt .* H̃αβ./2 + I
        Ub = 1im .* dt .* H̃αβ./2 + I
        #evolution
        y = Uf * ψt₀
        ψt₀ =  Ub\y             # solves Ub*x = y, new ψt0
        ProgressMeter.next!(p)
        #sleep(0.05)
        IJulia.clear_output(true)
    end
    ψt₀
end

"""
gaussKernel(Eα_r, Eαi, ΔEw)

Gaussian kernel. Cented at Eαi and width ΔEw.
"""
gaussKernel(Eα_r, Eαi, ΔEw) = exp.(-(Eα_r .- Eαi).^2 ./ΔEw.^2)/(ΔEw*sqrt(π))

"""
getSpectrum(coeffs, ΔEw, Eα, Eα_r)

Computes the power spectrum from the photon-electron interaction. With a gaussian kernel.
Where _ΔEw_ is the width for the gaussian kernel, _Eα_ is the spectrum from the interaction region and
_coeffs_ are the coefficients after the evolution of an initial state ψt₀.
"""
function getSpectrum(coeffs, ΔEw, Eα, Eα_r)
    nf = length(Eα)
    spectrum = zeros(nf)
    square_coefs = abs.(coeffs).^2
    for i in 1:nf
        spectrum += square_coefs[i] .* gaussKernel(Eα_r, Eα[i], ΔEw)
    end
    spectrum
end

function getSpectrumCs(square_coefs, ΔEw, Eα, Eα_r)
    nf = length(Eα)
    spectrum = zeros(length(Eα_r))
    for i in 1:nf
        spectrum += square_coefs[i] .* gaussKernel(Eα_r, Eα[i], ΔEw)
    end
    spectrum
end

χ0(xs, Ωs) = (Ωs/π)^(1/4) .* exp.(- Ωs*xs.^2.0./2)
χ1(xs, Ωs) = (4Ωs^3/π)^(1/4) .* xs .* exp.(- Ωs*xs.^2.0./2)
function χn(x; m = 2, t0 = 3, ω₀ = 0.57, Eα0 = -0.9)
    nr = m - 1
    T0 = 41.34*t0
    Ωs = T0^2/(4*log(2))
    x₀ = Eα0 .+ 2ω₀
    xs = (x .- x₀)
    a, b = χ1(xs, Ωs), χ0(xs, Ωs)
    if m<0
        return println("m must be greater than zero")
    elseif m==0
        return b
    elseif  m==1
        return a
    else
        for n in 1:nr
            a, b = √(2Ωs/(n+1)) .* xs .* a .-  √(n/(n+1)) .* b, a
        end
    end
    return a
end

"""
posnegE(Eαn)

Returns: Int64, bound position of positive values

"""
function posnegE(Eαn; valor = 0)
    nbound = 0
    for i in 1:length(Eαn)
        if Eαn[i] >valor
            nbound = i - 1
            break
        end
    end
    nbound
end

"""
randEαVαβ(Eα, Vαβ; γ = 0.88, upperlimit = 1)

Returns: Synthetic elements for new Hamiltonians (synthetic hamiltonians).
"""
function randEαVαβ(Eα, Vαβ; γ = 0.88, lowerlimit = -1, upperlimit = 1)
    Eαn = copy(Eα)
    Vαβn = copy(Vαβ)
    posEnerg = posnegE(Eαn)
    χ2, χ3, χ4 = 2*rand(3) .+ lowerlimit
    χ1 = (upperlimit - lowerlimit) * rand() .+ lowerlimit

    Eαn[2:posEnerg] .= 3^(χ1 - γ) * Eαn[2:posEnerg]
    Vαβn[1,1:posEnerg] .= 3^(χ2)* Vαβn[1,1:posEnerg]
    Vαβn[1:posEnerg,1] .= 3^(χ2) * Vαβn[1:posEnerg,1]
    Vαβn[1:posEnerg,posEnerg+1:end] .= 3^(χ3)* Vαβn[1:posEnerg,posEnerg+1:end]
    Vαβn[posEnerg+1:end,1:posEnerg] .= 3^(χ3)* Vαβn[posEnerg+1:end,1:posEnerg]
    Vαβn[posEnerg+1:end,posEnerg+1:end] .= 3^(χ4)* Vαβn[posEnerg+1:end,posEnerg+1:end]
    Eαn, Vαβn
end

function findCns(Eα_r, spectra, δE, Ntot, tt0, omega0; Eg = Eαr[1])
    Cᵢs = zeros(Ntot)
    for nn in 0:Ntot-1
        Cᵢs[nn+1] = sum(sqrt.(spectra) .* χn(Eα_r; m = nn, t0 = tt0, ω₀ = omega0, Eα0 = Eg))*δE
    end
    Cᵢs
end

function fitSpectro(Eα_r, coefref, tt0, omega0; Eg = Eαr[1])
    ncoefi = length(coefref)
    spectFitref = zeros(length(Eα_r))
    for i in 1:ncoefi
        spectFitref += coefref[i] * χn(Eα_r; m =i-1, t0=tt0,ω₀ =omega0,
            Eα0= Eg)
    end
    abs.(spectFitref).^2
end

function noisyStates(fψ̃α, finvψ̃α, fEα; npulses=1000, gref=false, tt0 = 3,
    τ0 = 0.5, omegac = 0.77, fflu = 1.0, vseed=1021)
    fstates = zeros(600, npulses)
    Random.seed!(vseed)
    for pulso in 1:npulses
        tlims_ex, pulse_ex =  FELpulse(t, Nt; T = tt0, τ = τ0, ω₀ = omegac,
        gauss_ref = gref, fluence = fflu)
        Δtij, posVecPot = getΔtsPulse(tlims_ex, pulse_ex, gridPot)
        ψt₀ = timeEvolution(Δtij, posVecPot, fψ̃α, finvψ̃α, fEα)
        fstates[:,pulso] = abs.(ψt₀).^2
    end
    fstates
end
function stateRefs(fψ̃α, finvψ̃α, fEα; tt0 = 3, omegac = 0.77, gaFlux = 1.0)
    tlims_ex, pulse_ex = pulseGauss(t, tt0, 1e10, omegac, gaFlux)
    Δtij, posVecPot = getΔtsPulse(tlims_ex, pulse_ex, gridPot)
    ψt₀ = timeEvolution(Δtij, posVecPot, fψ̃α, finvψ̃α, fEα)
    abs.(ψt₀).^2
end

function stateRefsDoubles(fψ̃α, finvψ̃α, fEα; tt0 = 3, omegac = 0.77,gaFlux = 1.0)
    tlims_ex, pulse_ex = pulsedouble(t, tt0,1e10, omegac, gaFlux)
    Δtij, posVecPot = getΔtsPulse(tlims_ex, pulse_ex, gridPot)
    ψt₀ = timeEvolution(Δtij, posVecPot, fψ̃α, finvψ̃α, fEα)
    abs.(ψt₀).^2
end

function stateRefsDoubles(fψ̃α, finvψ̃α, fEα, separation; tt0 = 3, omegac = 0.77,
    gaFlux = 1.0)
    tlims_ex, pulse_ex = pulsesSep(t,tt0,separation,omegac,gaFlux)
    Δtij, posVecPot = getΔtsPulse(tlims_ex, pulse_ex, gridPot)
    ψt₀ = timeEvolution(Δtij, posVecPot, fψ̃α, finvψ̃α, fEα)
    abs.(ψt₀).^2
end

function fbeststar(Eα_r, refspectra, specFitted,besT₀, Cns, difFit,
     δE, ncoef, cfreq)
    for i in 0.01:0.01:4
        Cnsx = findCns(Eα_r, refspectra, δE, ncoef, i, cfreq)
        specFittedx = fitSpectro(Eα_r, Cnsx, i, cfreq)
        difFitx = sum(abs.(refspectra .- specFittedx))*δE
        if difFitx <= difFit
            besT₀ = i
            Cns = Cnsx
            difFit = difFitx
            specFitted = specFittedx
        end
    end
    besT₀, Cns, difFit, specFitted
end

max_stars_repo_name: lazarusA/LightMatterInteraction

max_issues_repo_name: lazarusA/LightMatterInteraction

max_forks_repo_name: lazarusA/LightMatterInteraction


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule tree_sitter_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("tree_sitter")
JLLWrappers.@generate_main_file("tree_sitter", UUID("695101d5-1b81-5434-9529-70430e482717"))
end  # module tree_sitter_jll


text:
"""
    MOKernel

Abstract type for kernels with multiple outpus.
"""
abstract type MOKernel <: Kernel end

function _kernelmatrix_kron_helper(::Type{<:MOInputIsotopicByFeatures}, Kfeatures, Koutputs)
    return kron(Kfeatures, Koutputs)
end

function _kernelmatrix_kron_helper(::Type{<:MOInputIsotopicByOutputs}, Kfeatures, Koutputs)
    return kron(Koutputs, Kfeatures)
end

if VERSION >= v"1.6"
    function _kernelmatrix_kron_helper!(
        K, ::Type{<:MOInputIsotopicByFeatures}, Kfeatures, Koutputs
    )
        return kron!(K, Kfeatures, Koutputs)
    end

    function _kernelmatrix_kron_helper!(
        K, ::Type{<:MOInputIsotopicByOutputs}, Kfeatures, Koutputs
    )
        return kron!(K, Koutputs, Kfeatures)
    end
end


text:
## Signal operations

function ensemble_parameters(strategy::Function, returns::Returns, ensembleRange::UnitRange; method::Function = mean)::Weights
    weights = @pipe ensembleRange |>
        map(x -> strategy(returns, x), _) |>
        TimeSeries.merge(_...) |>
        TimeSeries.map((timestamp, values) -> (timestamp, repeat([method(values)], length(ensembleRange))), _)
    return weights[colnames(weights)[1]]
end

function merge_weights(weights::Array{Weights})::Weights
    common = TimeSeries.merge(weights...)
    target = common[colnames(common)[1]]
    for col in colnames(common)
        if col == colnames(common)[1] continue end # skip the first col, as it's already present
        target = target .+ common[col]
    end
    return target ./ length(weights)
end


text:
struct Chip

    number_of_mixing_compartments::Int64
    number_of_input_streams::Int64
    number_of_output_streams::Int64

    function Chip(compartments=0,inputs=0,outputs=0)
        this = new(compartments,inputs,outputs)
    end
end

struct Biochemistry
end

max_stars_repo_name: varnerlab/N-Compartment-Chip-Model

max_issues_repo_name: varnerlab/N-Compartment-Chip-Model

max_forks_repo_name: varnerlab/N-Compartment-Chip-Model


text:
using Random
using MLDatasets

# download mnist train and test dataset
Xtrain, Ytrain = MNIST.traindata()
Xtest,  Ytest  = MNIST.testdata()

# construct a simple mlp network
topology = [784, 256, 256, 256, 128, 10]
operator = [relu, relu, relu, relu, relu]
mlpmodel = MLP(topology, operator)

# or equivalently construct a mlp like this:
mlpmodel = Chain(
    Dense(784,256,relu),
    Dense(256,256,relu),
    Dense(256,256,relu),
    Dense(256,128,relu),
    Dense(128, 10,relu)
)

# training function
function train(Xtrain, Ytrain, Xtest, Ytest, mlpmodel)
    parameter = paramsof(mlpmodel)
    optimizer = Adam(parameter;lr=1e-3,b1=0.9,b2=0.996)
    epoch     = 2
    batchsize = 32
    filesnums = size(Xtrain,3)
    batchnums = div(filesnums,batchsize)
    lossvalue = zeros(batchnums,1)
    idx   = collect(1:filesnums)
    loss  = zeros(epoch,1)
    accc  = zeros(epoch,1)
    for e = 1:epoch
        shuffle!(idx)
        for b = 1:batchnums
            shift = batchsize*(b-1)
            input = zeros(784,0)
            label = zeros(10,0)
            for i = 1:batchsize
                feats = convert.(Float64, reshape(Xtrain[:,:,idx[i+shift]],784,1))
                numbs = tolabel(Ytrain[idx[i+shift]])
                input = hcat(input,feats)
                label = hcat(label,numbs)
            end
            outs = forward(mlpmodel, Variable(input))
            prob = softmax(outs; dims=1)
            Cost = crossEntropyCost(prob, Variable(label))
            backward()
            if b%8==0
                update!(optimizer,parameter;clipvalue=50)
                zerograds!(parameter)
            end
            lossvalue[b] = Cost / batchsize
        end
        accc[e] = test(Xtest, Ytest, mlpmodel)
        loss[e] = sum(lossvalue) / batchnums
        println("epoch[$(e)]","loss: ",loss[e] , " accuracy: ",accc[e])
    end
    return loss,accc
end

# testing function
function test(Xtest, Ytest, mlpmodel)
    acc = 0.0
    num = size(Xtest,3)
    for i = 1:num
        feats = convert.(Float64, reshape(Xtest[:,:,i],784,1))
        label = convert(Int,Ytest[i]) + 1
        result = argmax(softmax(predict(mlpmodel, feats),dims=1))[1]
        acc += (label==result) ? 1.0 : 0.0
    end
    return acc/num
end

function tolabel(x)
    y = zeros(10,1)
    y[x+1] = 1.0
    return y
end

loss, acc = train(Xtrain, Ytrain, Xtest, Ytest, mlpmodel)

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
module NewPkgEval

import Pkg.TOML
using Pkg
using Base: UUID
using Dates

downloads_dir(name) = joinpath(@__DIR__, "..", "deps", "downloads", name)
julia_path(ver) = joinpath(@__DIR__, "..", "deps", "julia-$ver")
versions_file() = joinpath(@__DIR__, "..", "deps", "Versions.toml")
registry_path(name) = joinpath(first(DEPOT_PATH), "registries", name)
registries_file() = joinpath(@__DIR__, "..", "deps", "Registries.toml")
builds_file() = joinpath(@__DIR__, "..", "deps", "Builds.toml")
log_path(julia) = joinpath(@__DIR__, "..", "logs/logs-$julia")

"""
    read_versions() -> Dict

Parse the `deps/Versions.toml` file containing download and verification information for
various versions of Julia.
"""
read_versions() = TOML.parse(read(versions_file(), String))

"""
    read_registries() -> Dict

Parse the `deps/Registries.toml` file containing a URL and packages to skip for listed
registries.
"""
read_registries() = TOML.parsefile(registries_file())

"""
    read_builds() -> Dict

Parse the `deps/Builds.toml` file containing download information for various Julia builds.
"""
read_builds() = TOML.parsefile(builds_file())

include("registry.jl")
include("julia.jl")
include("run.jl")

end # module


text:
using LinearAlgebra
export CER

struct CER
    mode::Matrix{Bool}
    α::Real
    directed::Bool
    self_loops::Bool
end

function CER(
    mode::Matrix{Int},
    α::Float64;
    directed::Bool=!issymmetric(mode),
    self_loops::Bool=any(diag(mode).>0)
    )
    @assert 0 < α < 1 "α must be within interval (0,1)."

    # If intger matrix given with just 0 and 1s make boolean
    if prod(x->x∈[0,1], mode)
        return CER(convert.(Bool,mode), α, directed, self_loops)
    else
        error("Entries must be boolean or 0/1 integers.")
    end
end

function CER(
    mode::Matrix{Bool},
    α::Float64;
    directed::Bool=!issymmetric(mode),
    self_loops::Bool=any(diag(mode))
    )
    @assert 0 < α < 1 "α must be within interval (0,1)."

    # If intger matrix given with just 0 and 1s make boolean
    return CER(mode, α, directed, self_loops)
end

# Sampling functions

function draw_sample!(
    out::Vector{Matrix},
    model::CER
    )
    for A in out
        copy!(A, model.mode)
        for c in eachcol(A)
            for i in eachindex(c)
                if rand() < model.α
                    c[i] = !c[i]
                end
            end
        end

    end
end

function draw_sample_no_copy!(
    out::Vector{Matrix{Bool}},
    model::CER
    )
    for A in out
        for c in eachcol(A)
            for i in eachindex(c)
                if rand() < model.α
                    c[i] = !c[i]
                end
            end
        end
    end
end

function draw_sample(
    model::CER,
    n::Int
    )
    out = [copy(model.mode) for i in 1:n]
    draw_sample_no_copy!(out, model)
    return out
end

max_stars_repo_path: src/Graph Models/CER/CER.jl

max_issues_repo_path: src/Graph Models/CER/CER.jl

max_forks_repo_path: src/Graph Models/CER/CER.jl


text:
module Bifurcation

using ArgParse

include("BifurcationLib.jl")

function main(args)

  s = ArgParseSettings("A demonstration of bifurcation and the resulting " *
    "chaos in the equilibrium points of the logistic map. This program " *
    "plots the equilibrium point of repeated evaluations of the logistic " *
    "map for various r.")

  @add_arg_table! s begin
    "x"
      arg_type = Float64
      default = 0.5
      help = "value to be used for first iteration of the logistic map"

    "r0"
      arg_type = Float64
      default = 0.0
      help = "initial r: the starting point for the plot"

    "n"
      arg_type = Int
      default = 100
      range_tester = (n->n>0)
      help = "number of iterations to execute logistic map before result is " *
        "considered to be at equilibrium"

    "m"
      arg_type = Int
      default = 100
      range_tester = (m->m>0)
      help = "number of iterations after equilibrium to plot"

    "r_inc"
      arg_type = Float64
      default = 0.01
      help = "amount to increment r by between tests"

    "r_final"
      arg_type = Float64
      default = 4.0
      help = "point at which we will stop testing and plotting"

  end

  parsed_args = parse_args(args, s)

  BifurcationLib.generate_and_plot(
    parsed_args["x"],
    parsed_args["r0"],
    parsed_args["n"],
    parsed_args["m"],
    parsed_args["r_inc"],
    parsed_args["r_final"],
  )
end

main(ARGS)

end # module


text:
## Load the packages.
using Tullio
using LoopVectorization
using Statistics
using BenchmarkTools
using PyPlot
using Interpolations

## Create the refine function.
function refine_field_nn!(hi, lo, n_hi, n_lo)
    # if n_hi ÷ n_lo == 2
    #     @tullio hi[2(i-1)+l] = lo[i] (l in 1:2)
    # elseif n_hi ÷ n_lo == 3
    #     @tullio hi[3(i-1)+l] = lo[i] (l in 1:3)
    # else
    #     throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    # end

    if n_hi ÷ n_lo == 2
        @tturbo for i in 0:2:size(hi, 1)-3
            i_lo = i÷2+2
            i_hi = i+2
            hi[i_hi  ] = lo[i_lo]
            hi[i_hi+1] = lo[i_lo]
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for i in 0:3:size(hi, 1)-3
            i_lo = i÷3+2
            i_hi = i+2
            hi[i_hi  ] = lo[i_lo]
            hi[i_hi+1] = lo[i_lo]
            hi[i_hi+2] = lo[i_lo]
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

function refine_field_int!(hi, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for i in 0:2:size(hi, 1)-3
            i_lo = i÷2+2
            i_hi = i+2
            hi[i_hi  ] = 3//4*lo[i_lo] + 1//4*lo[i_lo-1]
            hi[i_hi+1] = 3//4*lo[i_lo] + 1//4*lo[i_lo+1]
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for i in 0:3:size(hi, 1)-3
            i_lo = i÷3+2
            i_hi = i+2
            hi[i_hi  ] = 2//3*lo[i_lo] + 1//3*lo[i_lo-1]
            hi[i_hi+1] = lo[i_lo]
            hi[i_hi+2] = 2//3*lo[i_lo] + 1//3*lo[i_lo+1]
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

function refine_field_int2!(a_hi_gc, a_lo_gc, n_hi, n_lo)
    dx_hi = 1/n_hi; dx_lo = 1/n_lo
    x_hi = -dx_hi/2:dx_hi:1+dx_hi |> collect
    x_lo = -dx_lo/2:dx_lo:1+dx_lo |> collect
    interp = LinearInterpolation(x_lo, a_lo_gc)
    a_hi_gc[2:end-1] = interp(x_hi[2:end-1])
end

function refine_field_int3!(hi, hi_tmp, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for i in 0:2:size(hi, 1)-3
            i_lo = i÷2+2
            i_hi = i+2
            hi_tmp[i_hi  ] = lo[i_lo]
            hi_tmp[i_hi+1] = lo[i_lo]
        end

        hi_tmp[1] = hi_tmp[end-1]
        hi_tmp[end] = hi_tmp[2]

        @tturbo for i in 2:size(hi, 1)-1
            hi[i] = 1//4*hi_tmp[i-1] + 1//2*hi_tmp[i] + 1//4*hi_tmp[i+1]
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for i in 0:3:size(hi, 1)-3
            i_lo = i÷3+2
            i_hi = i+2
            hi_tmp[i_hi  ] = lo[i_lo]
            hi_tmp[i_hi+1] = lo[i_lo]
            hi_tmp[i_hi+2] = lo[i_lo]
        end

        hi_tmp[1] = hi_tmp[end-1]
        hi_tmp[end] = hi_tmp[2]

        @tturbo for i in 2:size(hi, 1)-1
            hi[i] = 1//3*hi_tmp[i-1] + 1//3*hi_tmp[i] + 1//3*hi_tmp[i+1]
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

## Set up the grids.
n_hi = 48
n_lo = n_hi ÷ 3

a_lo = rand(n_lo)
a_hi = zeros(n_hi)

a_hi_gc = zeros(n_hi+2)
a_lo_gc = zeros(n_lo+2)

a_hi_gc[2:n_hi+1] = a_hi[:]
a_lo_gc[2:n_lo+1] = a_lo[:]

a_hi_gc[1] = a_hi_gc[end-1]
a_hi_gc[end] = a_hi_gc[2]
a_lo_gc[1] = a_lo_gc[end-1]
a_lo_gc[end] = a_lo_gc[2]

dx_hi = 1/n_hi; dx_lo = 1/n_lo
x_hi = dx_hi/2:dx_hi:1 |> collect
x_lo = dx_lo/2:dx_lo:1 |> collect

## Compute.
@btime refine_field_nn!(a_hi_gc, a_lo_gc, n_hi, n_lo)
a_hi_nn = a_hi_gc[2:end-1]

@btime refine_field_int!(a_hi_gc, a_lo_gc, n_hi, n_lo)
a_hi_int = a_hi_gc[2:end-1]

@btime refine_field_int2!(a_hi_gc, a_lo_gc, n_hi, n_lo)
a_hi_int2 = a_hi_gc[2:end-1]

a_hi_tmp = zeros(size(a_hi_gc))
@btime refine_field_int3!(a_hi_gc, a_hi_tmp, a_lo_gc, n_hi, n_lo)
a_hi_int3 = a_hi_gc[2:end-1]

println(mean(a_lo) ≈ mean(a_hi_nn))
println(mean(a_lo) ≈ mean(a_hi_int))
println(mean(a_lo) ≈ mean(a_hi_int2))
println(mean(a_lo) ≈ mean(a_hi_int3))

println(mean(a_lo), " ", mean(a_hi_int), " ", mean(a_hi_int3))
println(a_hi_int3 .- a_hi_int)

## Plot the output.
figure()
#plot(x_hi, a_hi_nn, "C0-o")
#plot(x_hi, a_hi_int, "C1-o")
#plot(x_hi, a_hi_int2, "C2-^")
plot(x_hi, a_hi_int3, "C3-+")
plot(x_lo, a_lo, "k:")
display(gcf())
show()


text:
# Module for creating tables
module Table

"""
Generate table typeset in LaTeX

A table that looks like

        a | b |
        --|---|
        0 | - |
        0 |   |

is created by passing header = ["a", "b"], a = [0, 0], b = [Inf, 0],
variables = [a, b] to the function. The generated table is returned as a string.
"""
function latex(header, variables; format="c")
        @assert length(header) == length(variables)
        n = length(header)
        fmt = repeat(format, n)
        table = []
        txt = "\\begin{tabular}{$fmt}\n"
        append!(table, txt)
        m = length(variables[1])
        for j in 1:n
                @assert length(variables[j]) == m
                hi = header[j]
                if j != n
                        txt = "$hi & "
                else
                        txt = "$hi \\\\ \n"
                end
                append!(table, txt)
        end
        txt = "\\hline\n"
        append!(table, txt)

        for i in 1:m
                for j in 1:n
                        var = variables[j][i]
                        if var != Inf
                                val = "\$$var\$"
                        else
                                val = "-"
                        end
                        if j != n
                                txt = "$val & "
                        else
                                txt = "$val"
                        end
                        append!(table, txt)
                end

                if i != m
                        txt = " \\\\ \n"
                else
                        txt = "\n"
                end
                append!(table, txt)
        end

        txt = "\\end{tabular}\n"
        append!(table, txt)
        return table

end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# 2D FWI gradient test with 4 sources
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

parsed_args = parse_commandline()

nlayer = parsed_args["nlayer"]
tti = parsed_args["tti"]
fs =  parsed_args["fs"]

### Model
model, model0, dm = setup_model(parsed_args["tti"], 4)
q, srcGeometry, recGeometry, info = setup_geom(model)
dt = srcGeometry.dt[1]

###################################################################################################

@testset "TWRI gradient test w.r.t m with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
	optw = TWRIOptions(grad_corr=false, comp_alpha=false, weight_fun=nothing, eps=0, params=:m)
	# Gradient test
	h = 5f-2
	maxiter = 5
	err1 = zeros(maxiter)
	err2 = zeros(maxiter)
	h_all = zeros(maxiter)
	modelH = deepcopy(model0)

	# Observed data
	opt = Options(sum_padding=true, free_surface=parsed_args["fs"])
	F = judiModeling(info, model, srcGeometry, recGeometry; options=opt)
	F0 = judiModeling(info, model0, srcGeometry, recGeometry; options=opt)
	d = F*q
	d0 = F0*q
	y = 2.5f0*(d0 - d)

	# FWI gradient and function value for m0
	Jm0, gradm = twri_objective(model0, q, d, y; options=opt, optionswri=optw)

	dJ = dot(gradm, dm)
	@printf("Perturbation size is %2.2e and reference objective function is %2.4e \n", dJ, Jm0)

	for j=1:maxiter
		# FWI gradient and function falue for m0 + h*dm
		modelH.m = model0.m + h*dm
		Jm, _ = twri_objective(modelH, q, d, y; options=opt, optionswri=optw)

		# Check convergence
		err1[j] = abs(Jm - Jm0)
		err2[j] = abs(Jm - Jm0 - h*dJ)
		j == 1 ? prev = 1 : prev = j - 1
		@printf("h = %2.2e, phi = %2.4e, e1 = %2.2e, rate = %2.2e", h, Jm, err1[j], err1[prev]/err1[j])
		@printf(", e2  = %2.2e, rate = %2.2e \n", err2[j], err2[prev]/err2[j])
		h_all[j] = h
		h = h * .8f0
	end
	# CHeck convergence rates
	rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
	rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)
	@test isapprox(rate_1, 1.25f0; rtol=5f-2)
	@test isapprox(rate_2, 1.5625f0; rtol=5f-2)
end

@testset "TWRI gradient test w.r.t y with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
	optw = TWRIOptions(grad_corr=false, comp_alpha=false, weight_fun=nothing, eps=0, params=:y)
	# Gradient test
	h = 5f-2
	maxiter = 5
	err1 = zeros(maxiter)
	err2 = zeros(maxiter)
	h_all = zeros(maxiter)
	modelH = deepcopy(model0)

	# Observed data
	opt = Options(sum_padding=true, free_surface=parsed_args["fs"])
	F = judiModeling(info, model, srcGeometry, recGeometry; options=opt)
	F0 = judiModeling(info, model0, srcGeometry, recGeometry; options=opt)
	d = F*q
	d0 = F0*q
	y = .1f0 * (d0 - d)
	dy = .05f0*(d0 - d)

	# FWI gradient and function value for m0
	Jm0, grady = twri_objective(model0, q, d, y; options=opt, optionswri=optw)

	dJ = dot(grady, dy)
	@printf("Perturbation size is %2.4e and reference objective function is %2.4e \n", dJ, Jm0)

	for j=1:maxiter
		# FWI gradient and function falue for m0 + h*dm
		yloc = y + h * dy
		Jm, _ = twri_objective(model0, q, d, yloc;options=opt, optionswri=optw)

		# Check convergence
		err1[j] = abs(Jm - Jm0)
		err2[j] = abs(Jm - Jm0 - h*dJ)
		j == 1 ? prev = 1 : prev = j - 1
		@printf("h = %2.2e, phi = %2.4e, e1 = %2.2e, rate = %2.2e", h, Jm, err1[j], err1[prev]/err1[j])
		@printf(", e2  = %2.2e, rate = %2.2e \n", err2[j], err2[prev]/err2[j])
		h_all[j] = h
		h = h * .8f0
	end
	# CHeck convergence rates
	rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
	rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)
	@test isapprox(rate_1, 1.25f0; rtol=5f-2)
	@test isapprox(rate_2, 1.5625f0; rtol=5f-2)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
module GaussianMixture
    using DataModel
    using MLUtilities

    export GaussianMixtureModel

    type GaussianMixtureModel <: AbstractDataModel
        #currently just 1D
        mu::Array{Float64}
        var::Array{Float64}
        weights::Array{Float64}
    end

    normldf(x,mu,var) =  -.5*log(2.0*pi*var) -.5*((x-mu).*(x-mu))/var
    normpdf(x,mu,var) = exp(normldf(x,mu,var))

    function DataModel.getgrad(dm::GaussianMixtureModel)
        function gradient(x)
            z = exp(logsumexp([log(dm.weights[i])+normldf(x,dm.mu[i],dm.var[i]) for i in 1:length(dm.mu)]...))
            d = sum([dm.weights[i].*normpdf(x,dm.mu[i],dm.var[i]).*(dm.mu[i]-x)./dm.var[i] for i in 1:length(dm.mu)])
            return d./z
        end
        return gradient
    end

    function DataModel.getllik(dm::GaussianMixtureModel)
        function logdensity(x)
            logsumexp([log(dm.weights[i])+normldf(x,dm.mu[i],dm.var[i]) for i in 1:length(dm.mu)]...)[1]
        end
        return logdensity
    end
end


text:
function gradient(var::CellVar, grid::TreeGrid; filter = active)
    grad = FaceVar(zeros(eltype(var), size(FaceVar, grid)))
    @inbounds for face ∈ faces(grid, filter = filter)
        val = zero(eltype(var))
        for cell ∈ FullyThreadedTree.cells(face)
            val -= FullyThreadedTree.face_area(cell, face) * normal_sign(cell, face) * var[cell]
        end
        grad[face] = val /= FullyThreadedTree.cell_volume(face)
    end
    return grad
end

function divergence(var::FaceVar, grid::TreeGrid; filter = active)
    div = CellVar(zeros(eltype(var), size(CellVar, grid)))
    @inbounds for cell ∈ cells(grid, filter = filter)
        val = zero(eltype(var))
        for face ∈ FullyThreadedTree.faces(cell)
            val += FullyThreadedTree.area(face) * normal_sign(cell, face) * var[face]
        end
        div[cell] = val /= FullyThreadedTree.volume(cell)
    end
    return div
end

function gradient(var::CellVar, grid::CartesianGrid)

end


text:
module MyGEMM

include("types.jl")
include("utils.jl")
include("transpose_no_packing.jl")
include("transpose_multithreading.jl")

end

max_stars_repo_name: bischtob/MyGEMM

max_issues_repo_name: bischtob/MyGEMM

max_forks_repo_name: bischtob/MyGEMM


text:
# License is MIT: https://github.com/genkuroki/MyUtils.jl/blob/main/LICENSE

"""
    savevar(fn, x)

saves the value of `x` to the file `fn`, where `fn` is the filename string of the file.
"""
savevar(fn, x) = write(fn, string(x))

"""
    loadvar(fn)

loads the file `fn` (the filename string of the file) and `Meta.parse |> eval`.
"""
loadvar(fn) = read(fn, String) |> Meta.parse |> eval

"""
    dir_savevar[]

is the default directory to which `@savevar` saves the values of variables.
"""
const dir_savevar = Ref(".")

"""
    fn_savevar(x::Symbol)

is the filename string to which `@savevar` saves the value of a variable.
"""
fn_savevar(x::Symbol) = joinpath(dir_savevar[], string(x) * ".txt")

"""
    @savevar(args...)

saves the variables in args to the corresponding textfiles.

Example: `@savevar A B C` saves the variables `A`, `B`, `C` to textfiles.
"""
macro savevar(args...)
    A = [:(savevar($(fn_savevar(x)), $(esc(x)))) for x in args]
    quote $(A...); nothing end
end

"""
    @loadvar(args...)

loads the values from the textfiles corresponding to `args`.
If `length(args)` is greater than 1, then it returns the tuple of the values.

Example: `a, b, c = @loadvar A B C` loads
the values of `A`, `B`, `C` in textfiles to the variables `a`, `b`, `c`.
"""
macro loadvar(args...)
    if length(args) == 1
        x = args[1]
        :(loadvar($(fn_savevar(x))))
    else
        A = [:(loadvar($(fn_savevar(x)))) for x in args]
        :(($(A...),))
    end
end


text:
using Distributions
using Cytof5, Random
using JLD2, FileIO
import BSON
include("salso.jl")

# Plotting
using PyCall
matplotlib = pyimport("matplotlib")
plt = pyimport("matplotlib.pyplot")
matplotlib.use("Agg")
# Load current dir
pushfirst!(PyVector(pyimport("sys")."path"), "../vb")
plot_yz = pyimport("plot_yz")
blue2red = pyimport("blue2red")

# Path to scratch
RESULTS_DIR = "/scratchdata/alui2/cytof/results/"

# y heatmap colorbar range
VLIM = (-4, 4)

# get posterior samples by symbol
getPosterior(sym::Symbol, monitor) = [m[sym] for m in monitor]

# read a single object from JLD2
function loadSingleObj(objPath)
  data = load(objPath)
  return data[collect(keys(data))[1]]
end

# plot y/z
function make_yz(y, Zs, Ws, lams, imgdir; w_thresh=.01, lw=3,
                 Z_true=nothing, fs_y=22, fs_z=22, fs_ycbar=22, fs_zcbar=22)
  mkpath(imgdir)
  I = length(y)
  for i in 1:I
    idx_best = estimate_ZWi_index(Zs, Ws, i)

    Zi = Int.(Zs[idx_best])
    Wi = Float64.(Ws[idx_best][i, :])
    lami = Int64.(lams[idx_best][i])
    yi = Float64.(y[i])

    # plot Yi, lami
    plt.figure(figsize=(8, 8))
    plot_yz.plot_y(yi, Wi, lami, vlim=VLIM, cm=blue2red.cm(9), lw=lw,
                   fs_xlab=fs_y, fs_ylab=fs_y, fs_lab=fs_y, fs_cbar=fs_ycbar)
    plt.savefig("$(imgdir)/y$(i).pdf", bbox_inches="tight")
    plt.close()

    # plot Zi, Wi
    plt.figure(figsize=(8, 8))
    plot_yz.plot_Z(Zi, Wi, lami, w_thresh=w_thresh, add_colorbar=false,
                   fs_lab=fs_z, fs_celltypes=fs_z, fs_markers=fs_z, fs_cbar=fs_zcbar)
    plt.savefig("$(imgdir)/Z$(i).pdf", bbox_inches="tight")
    plt.close()
  end

  if Z_true != nothing
    # plot Z true
    plt.figure(figsize=(8, 8))
    plot_yz.plot_Z_only(Z_true, fs=fs_z,
                        xlab="cell phenotype", ylab="markers", rotate_xticks=false)
    plt.savefig("$(imgdir)/Z_true.pdf", bbox_inches="tight")
    plt.close()

    # plot ZT true
    plt.figure(figsize=(8, 8))
    plot_yz.plot_Z_only(Z_true', fs=fs_z,
                        xlab="markers", ylab="cell phenotype")
    plt.savefig("$(imgdir)/ZT_true.pdf", bbox_inches="tight")
    plt.close()
  end
end

function make_yz(path_to_output)
  if occursin(".jld2", path_to_output)
    output = load(path_to_output)
    path_to_data = "$(dirname(path_to_output))/reduced_data/reduced_cb.jld2"
    y = loadSingleObj(path_to_data)
    out = output["out"]
  elseif occursin(".bson", path_to_output)
    output = BSON.load(path_to_output)
    y = output[:simdat][:y]
    out = output[:out]
  else
    println("Neither jld2 nor bson?!")
  end

  results_dir = "$(dirname(path_to_output))/img"
  results_dir = replace(results_dir, RESULTS_DIR => "results/")

  Zs = getPosterior(:Z, out[1])
  Ws = getPosterior(:W, out[1])
  lams = getPosterior(:lam, out[1])

  if haskey(output, :simdat)
    Z_true = output[:simdat][:Z]
  else
    Z_true = nothing
  end
  make_yz(y, Zs, Ws, lams, results_dir, Z_true=Z_true)
end

# MAIN
# path_to_output = "$(RESULTS_DIR)/cb/best/output.jld2"
@time for (root, dirs, files) in walkdir(RESULTS_DIR)
  for file in files
    if occursin("output.jld2", file) || occursin("output.bson", file)
    # if occursin("output.bson", file)
      path_to_output = "$(root)/$(file)"
      println("Current: $(path_to_output)")
      make_yz(path_to_output)
    end
  end
end


text:
const BAM_FPAIRED        =    1 # the read is paired in sequencing, no matter whether it is mapped in a pair
const BAM_FPROPER_PAIR   =    2 # the read is mapped in a proper pair
const BAM_FUNMAP         =    4 # the read itself is unmapped; conflictive with BAM_FPROPER_PAIR
const BAM_FMUNMAP        =    8 # the mate is unmapped
const BAM_FREVERSE       =   16 # the read is mapped to the reverse strand
const BAM_FMREVERSE      =   32 # the mate is mapped to the reverse strand
const BAM_FREAD1         =   64 # this is read1
const BAM_FREAD2         =  128 # this is read2
const BAM_FSECONDARY     =  256 # not primary alignment
const BAM_FQCFAIL        =  512 # QC failure
const BAM_FDUP           = 1024 # optical or PCR duplicate
const BAM_FSUPPLEMENTARY = 2048 # supplementary alignment

const CIGAR_M  = 0
const CIGAR_I  = 1
const CIGAR_D  = 2
const CIGAR_N  = 3
const CIGAR_S  = 4
const CIGAR_H  = 5
const CIGAR_P  = 6
const CIGAR_Eq = 7
const CIGAR_X  = 8


text:
@testset "Type Stability" begin
    function rosenbrock(x::Vector{T}) where T
        o = one(T)
        c = convert(T,100)
        return (o - x[1])^2 + c * (x[2] - x[1]^2)^2
    end

    function rosenbrock_gradient!(storage::Vector{T}, x::Vector{T}) where T
        o = one(T)
        c = convert(T,100)
        storage[1] = (-2*o) * (o - x[1]) - (4*c) * (x[2] - x[1]^2) * x[1]
        storage[2] = (2*c) * (x[2] - x[1]^2)
    end

    function rosenbrock_hessian!(storage::Matrix{T}, x::Vector{T}) where T
        o = one(T)
        c = convert(T,100)
        f = 4*c
        storage[1, 1] = (2*o) - f * x[2] + 3 * f * x[1]^2
        storage[1, 2] = -f * x[1]
        storage[2, 1] = -f * x[1]
        storage[2, 2] = 2*c
    end

    for method in (NelderMead(),
                   SimulatedAnnealing(),
                   BFGS(),
                   ConjugateGradient(),
                   GradientDescent(),
                   MomentumGradientDescent(),
                   AcceleratedGradientDescent(),
                   LBFGS(),
                   Newton())
        for T in (Float32, Float64)
            result = optimize(rosenbrock, rosenbrock_gradient!, rosenbrock_hessian!, fill(zero(T), 2), method)
            @test eltype(Optim.minimizer(result)) == T
        end
    end
end


text:
@testset "Display" begin
    function count_filled_rect(svg::AbstractString)
        n = 0
        i = firstindex(svg)
        while true
            r = findnext(r"<rect[^>]*\sfill=\"#[0-9A-F]{6}\"[^>]+>", svg, i)
            r === nothing && return n
            n += 1
            i = last(r)
        end
    end

    function count_colored_stroke(svg::AbstractString)
        n = 0
        i = firstindex(svg)
        while true
            r = findnext(r"<path[^>]*\sstroke=\"#[0-9A-F]{6}\"[^>]+>", svg, i)
            r === nothing && return n
            n += 1
            i = last(r)
        end
    end

    # test ability to add to previously shown array of colors - issue #328
    a = [colorant"white", colorant"red"]
    buf = IOBuffer()
    show(buf, "image/svg+xml", a)
    take!(buf)
    push!(a, colorant"blue")
    @test length(a) == 3

    # the following tests depend on the constants in "src/display.jl"

    # single color
    # ------------
    show(buf, "image/svg+xml", colorant"hsl(120, 100%, 50%)")
    single = String(take!(buf))
    @test occursin(r"<svg[^>]+\swidth=\"25mm\"", single)
    @test occursin(r"<svg[^>]+\sheight=\"25mm\"", single)
    @test occursin(r"\sfill=\"#00FF00\"", single)
    @test occursin(r"</svg>$", single)

    # vectors
    # -------
    # square swatches
    v3 = [colorant"red", colorant"green", colorant"blue"]
    show(buf, "image/svg+xml", v3)
    vector3 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"75mm\"", vector3)
    @test occursin(r"<svg[^>]*\sheight=\"25mm\"", vector3)
    @test occursin(r"<rect[^>]*\swidth=\".96\"", vector3) # 0.96 = 24mm/25mm
    @test occursin(r"<rect[^>]*\sheight=\".96\"", vector3) # height = width
    @test occursin(r"<rect[^>]*\sfill=\"#FF0000\"", vector3)
    @test occursin(r"<rect[^>]*\sfill=\"#008000\"", vector3)
    @test occursin(r"<rect[^>]*\sfill=\"#0000FF\"", vector3)
    @test occursin(r"</svg>$", vector3)

    # rectangle swatches
    show(buf, "image/svg+xml", colormap("RdBu", 100))
    vec100 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"180mm\"", vec100) # max_width
    @test occursin(r"<svg[^>]*\sheight=\"25mm\"", vec100)
    @test occursin(r"<rect[^>]*\swidth=\"1\"", vec100) # no padding
    @test occursin(r"<rect[^>]*\sheight=\".96\"", vec100) # 0.96 = 24mm/25mm
    @test count_filled_rect(vec100) == 100

    # issue #317
    # limits max swatches
    # the size warning should be suppressed
    show(buf, MIME"image/svg+xml"(), v3, max_swatches=2) # not string but MIME
    lim_vector3 = String(take!(buf))
    # since `max_swatches`(=2) is less practical,
    # the width/height is greater than `max_swatch_size`
    @test occursin(r"<svg[^>]*\swidth=\"180mm\"", lim_vector3) # max_width
    @test occursin(r"<svg[^>]*\sheight=\"60mm\"", lim_vector3) # 180mm/3 * 1
    # decimation factor `d` should be 2
    # viewBox is based on the original size (=3x1)
    @test occursin(r"<svg[^>]*\sviewBox=\"0 -1 3 1\"", lim_vector3)
    # therefore, each swatch is contrarily enlarged due to the swatch reduction.
    # stroke width is equivalent to swatch height which is equal to `d`(=2)
    @test occursin(r"<svg[^>]*\sstroke-width=\"2\"", lim_vector3)
    # swatch width is equal to `d`(=2)
    @test occursin(r"<path[^>]*\sd=\"[^\"]+h2\"", lim_vector3) # no padding
    # the mean RGB of first 2 elements (red and green)
    @test occursin(r"stroke=\"#804000\"[^/]*/>\s*<path\s", lim_vector3)
    # the RGB of last element (blue)
    @test occursin(r"stroke=\"#0000FF\"[^/]*/>\s*</svg>", lim_vector3)
    @test count_colored_stroke(lim_vector3) == 2 # <= max_swatches

    # m * n matrices
    # --------------
    # n * max_swatch_size <= max_width &&
    # m * max_swatch_size <= max_height
    # square swatches
    m2x2 = [colorant"white" colorant"blue"
            colorant"black" colorant"cyan"]
    show(buf, "image/svg+xml", m2x2)
    mat2x2 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"50mm\"", mat2x2)
    @test occursin(r"<svg[^>]*\sheight=\"50mm\"", mat2x2)
    @test occursin(r"<rect[^>]*\swidth=\".96\"", mat2x2) # 0.96 = 24mm/25mm
    @test occursin(r"<rect[^>]*\sheight=\".96\"", mat2x2) # height = width
    @test occursin(r"<rect[^>]*\sfill=\"#FFFFFF\"", mat2x2)
    @test occursin(r"<rect[^>]*\sfill=\"#0000FF\"", mat2x2)
    @test occursin(r"<rect[^>]*\sfill=\"#000000\"", mat2x2)
    @test occursin(r"<rect[^>]*\sfill=\"#00FFFF\"", mat2x2)
    @test occursin(r"</svg>$", vector3)

    # n * max_swatch_size > max_width &&
    # m * max_swatch_size <= max_height
    # square swatches
    show(buf, "image/svg+xml", rand(RGB, 2, 9))
    mat2x9 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"180mm\"", mat2x9) # max_width
    @test occursin(r"<svg[^>]*\sheight=\"40mm\"", mat2x9) # 40mm = 180mm/9 * 2
    @test occursin(r"<rect[^>]*\swidth=\".95\"", mat2x9) # 0.95 = 19mm/20mm
    @test occursin(r"<rect[^>]*\sheight=\".95\"", mat2x9) # height = width
    @test count_filled_rect(mat2x9) == 18

    # n * max_swatch_size <= max_width &&
    # m * max_swatch_size > max_height
    # square swatches
    show(buf, "image/svg+xml", rand(RGB, 10, 2))
    mat10x2 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"30mm\"", mat10x2) # 30mm = 150mm/10 * 2
    @test occursin(r"<svg[^>]*\sheight=\"150mm\"", mat10x2) # max_height
    @test occursin(r"<rect[^>]*\swidth=\".93\"", mat10x2) # 0.93 = 14mm/15mm
    @test occursin(r"<rect[^>]*\sheight=\".93\"", mat10x2) # height = width
    @test count_filled_rect(mat10x2) == 20

    # issue #341
    # n * max_swatch_size > max_width &&
    # m * max_swatch_size > max_height &&
    # max_width / n * m >= max_swatch_size
    # square swatches (strokes)
    show(buf, "image/svg+xml", rand(RGB, 28, 181))
    mat28x181 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"180mm\"", mat28x181) # max_width
    @test occursin(r"<svg[^>]*\sheight=\"27.85mm\"", mat28x181) # 180mm/181 * 28
    @test occursin(r"<path[^>]*\sd=\"[^\"]+h1\"", mat28x181) # no padding
    @test count_colored_stroke(mat28x181) == 28*181

    # issue #341
    # n * max_swatch_size > max_width &&
    # m * max_swatch_size > max_height &&
    # max_height / m * n < max_swatch_size
    # square swatches (strokes)
    show(buf, "image/svg+xml", rand(RGB, 200, 28))
    mat200x28 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"21mm\"", mat200x28) # 150mm/200*28
    @test occursin(r"<svg[^>]*\sheight=\"150mm\"", mat200x28) # max_height
    @test occursin(r"<path[^>]*\sd=\"[^\"]+h1\"", mat200x28) # no padding
    @test count_colored_stroke(mat200x28) == 200*28

    # issue #317
    # implicitly limits max swatches (but explicitly warns)
    # m * n > default_max_swatches
    # square swatches (strokes) with reduction
    m184x360 = [HSV(h, s / 183, 0.8) for s = 0:183, h = 0:359]
    # warning of size
    warning = r"the large size \(184×360\)"
    @test_logs (:warn, warning) show(buf, "image/svg+xml", m184x360)
    mat184x360 = String(take!(buf))
    @test occursin(r"<svg[^>]*\swidth=\"180mm\"", mat184x360) # max_width
    @test occursin(r"<svg[^>]*\sheight=\"92mm\"", mat184x360) # 180mm/360 * 184
    # decimation factor should be 3
    @test occursin(r"<svg[^>]*\sviewBox=\"0 -1.5[0]* 360 184\"", mat184x360)
    @test occursin(r"<svg[^>]*\sstroke-width=\"3\"", mat184x360) # `d`(=3)
    @test occursin(r"<path[^>]*\sd=\"[^\"]+h3\"", mat184x360) # `d`(=3)
    # the mean RGB of 3 (not 9) elements in the bottom-right corner
    # "CC0007" == hex(RGB(0.8, 0.0, (0.12/3 + 0.08/3 + 0.04/3)/3))
    @test occursin(r"stroke=\"#CC0007\"[^/]*/>\s*</svg>", mat184x360)
    @test count_colored_stroke(mat184x360) == ceil(184/3)*ceil(360/3)

    # issue #317
    # limits max swatches
    # the size warning should be suppressed
    show(buf, MIME"image/svg+xml"(), m2x2, max_swatches=3) # not string but MIME
    lim_mat2x2= String(take!(buf))
    # since `max_swatches`(=3) is less practical,
    # the width/height is greater than `max_swatch_size`
    @test occursin(r"<svg[^>]*\swidth=\"150mm\"", lim_mat2x2) # 180mm/2 * 2
    @test occursin(r"<svg[^>]*\sheight=\"150mm\"", lim_mat2x2) # max_height
    # decimation factor `d` should be 2
    @test occursin(r"<svg[^>]*\sviewBox=\"0 -1 2 2\"", lim_mat2x2)
    @test occursin(r"<svg[^>]*\sstroke-width=\"2\"", lim_mat2x2) # `d`(=2)
    @test occursin(r"<path[^>]*\sd=\"[^\"]+h2\"", lim_mat2x2) # `d`(=2)
    # the mean RGB of 2x2 elements
    @test occursin(r"stroke=\"#4080BF\"[^/]*/>\s*</svg>", lim_mat2x2)
    @test count_colored_stroke(lim_mat2x2) == 1 # <= max_swatches
end