text:
module MPS

import ..FileFormats

import MathOptInterface
const MOI = MathOptInterface

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

MOI.supports(::Model, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = false

struct Options
    warn::Bool
end

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

"""
    Model(; kwargs...)

Create an empty instance of FileFormats.MPS.Model.

Keyword arguments are:

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

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

@enum(VType, VTYPE_CONTINUOUS, VTYPE_INTEGER, VTYPE_BINARY)

# ==============================================================================
#
#   Base.write
#
# ==============================================================================

"""
    Base.write(io::IO, model::FileFormats.MPS.Model)

Write `model` to `io` in the MPS file format.
"""
function Base.write(io::IO, model::Model)
    options = get_options(model)
    FileFormats.create_unique_names(
        model;
        warn = options.warn,
        replacements = Function[s -> replace(s, ' ' => '_')]
    )
    ordered_names = String[]
    names = Dict{MOI.VariableIndex, String}()
    for x in MOI.get(model, MOI.ListOfVariableIndices())
        n = MOI.get(model, MOI.VariableName(), x)
        push!(ordered_names, n)
        names[x] = n
    end
    write_model_name(io, model)
    write_rows(io, model)
    discovered_columns = write_columns(io, model, ordered_names, names)
    write_rhs(io, model)
    write_ranges(io, model)
    write_bounds(io, model, discovered_columns, ordered_names, names)
    write_sos(io, model, names)
    println(io, "ENDATA")
    return
end

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

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

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

const SET_TYPES = (
    (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, model, S, sense_char)
    for index in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}()
    )
        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
    return
end

function write_rows(io::IO, model::Model)
    println(io, "ROWS\n N  OBJ")
    for (set_type, sense_char) in SET_TYPES
        _write_rows(io, model, set_type, sense_char)
    end
    return
end

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

function _list_of_integer_variables(model, names, integer_variables, S)
    for index in MOI.get(
        model, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}()
    )
        v_index = MOI.get(model, MOI.ConstraintFunction(), index)
        push!(integer_variables, names[v_index.variable])
    end
    return
end

function list_of_integer_variables(model::Model, names)
    integer_variables = Set{String}()
    for S in (MOI.ZeroOne, MOI.Integer)
        _list_of_integer_variables(model, names, integer_variables, S)
    end
    return integer_variables
end

function extract_terms(
    v_names::Dict{MOI.VariableIndex, String},
    coefficients::Dict{String, Vector{Tuple{String, Float64}}},
    row_name::String,
    func::MOI.ScalarAffineFunction,
    discovered_columns::Set{String},
    multiplier::Float64 = 1.0,
)
    for term in func.terms
        variable_name = v_names[term.variable_index]
        push!(
            coefficients[variable_name],
            (row_name, multiplier * term.coefficient)
        )
        push!(discovered_columns, variable_name)
    end
    return
end

function _collect_coefficients(
    model,
    S,
    v_names::Dict{MOI.VariableIndex, String},
    coefficients::Dict{String, Vector{Tuple{String, Float64}}},
    discovered_columns::Set{String},
)
    for index in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}(),
    )
        row_name = MOI.get(model, MOI.ConstraintName(), index)
        func = MOI.get(model, MOI.ConstraintFunction(), index)
        extract_terms(
            v_names, coefficients, row_name, func, discovered_columns
        )
    end
    return
end

function write_columns(io::IO, model::Model, ordered_names, names)
    coefficients = Dict{String, Vector{Tuple{String, Float64}}}(
        n => Tuple{String, Float64}[] for n in ordered_names
    )
    # 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}()
    # Build constraint coefficients
    for (S, _) in SET_TYPES
        _collect_coefficients(model, S, names, coefficients, discovered_columns)
    end
    # Build objective
    # MPS doesn't support maximization so we flip the sign on the objective
    # coefficients.
    s = MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE ? -1.0 : 1.0
    obj_func = MOI.get(
        model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
    )
    extract_terms(names, coefficients, "OBJ", obj_func, discovered_columns, s)
    integer_variables = list_of_integer_variables(model, names)
    println(io, "COLUMNS")
    int_open = false
    for variable in ordered_names
        is_int = variable in integer_variables
        if is_int && !int_open
            println(io, "    MARKER    'MARKER'                 'INTORG'")
            int_open = true
        elseif !is_int && int_open
            println(io, "    MARKER    'MARKER'                 'INTEND'")
            int_open = false
        end
        for (constraint, coefficient) in coefficients[variable]
            print(io, "     ")
            print(io, rpad(variable, 8))
            print(io, " ")
            print(io, rpad(constraint, 8))
            print(io, " ")
            Base.Grisu.print_shortest(io, coefficient)
            println(io)
        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, model, S, sense_char)
    for index in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}()
    )
        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

function write_rhs(io::IO, model::Model)
    println(io, "RHS")
    for (set_type, sense_char) in SET_TYPES
        _write_rhs(io, model, set_type, sense_char)
    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::Model)
    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
        println(io, " FR bounds    ", var_name)
    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(x::Tuple{Float64, Float64, VType}, set::MOI.GreaterThan)
    return (max(x[1], set.lower), x[2], x[3])
end

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

function update_bounds(x::Tuple{Float64, Float64, VType}, set::MOI.Interval)
    return (set.lower, set.upper, x[3])
end

function update_bounds(x::Tuple{Float64, Float64, VType}, set::MOI.EqualTo)
    return (set.value, set.value, x[3])
end

function update_bounds(x::Tuple{Float64, Float64, VType}, set::MOI.ZeroOne)
    return (x[1], x[2], VTYPE_BINARY)
end

function _collect_bounds(bounds, model, S, names)
    for index in MOI.get(
        model, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}()
    )
        func = MOI.get(model, MOI.ConstraintFunction(), index)
        set = MOI.get(model, MOI.ConstraintSet(), index)::S
        name = names[func.variable]
        bounds[name] = update_bounds(bounds[name], set)
    end
    return
end

function write_bounds(
    io::IO, model::Model, discovered_columns::Set{String}, ordered_names, names
)
    println(io, "BOUNDS")
    bounds = Dict{String, Tuple{Float64, Float64, VType}}(
        n => (-Inf, Inf, VTYPE_CONTINUOUS) for n in ordered_names
    )
    for S in (
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
        MOI.Interval{Float64},
        MOI.ZeroOne,
    )
        _collect_bounds(bounds, model, S, names)
    end
    for var_name in ordered_names
        lower, upper, vtype = bounds[var_name]
        if !(var_name in discovered_columns)
            @warn(
                "Variable $var_name is mentioned in BOUNDS, but is not " *
                "mentioned in the COLUMNS section. We are ignoring it."
            )
            continue
        elseif vtype == VTYPE_BINARY
            println(io, " BV bounds    ", var_name)
        end
        write_single_bound(io, var_name, lower, upper)
    end
    return
end

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

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

function write_sos(io::IO, model::Model, names)
    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, names)
                idx += 1
            end
        end
    end
    return
end

# ==============================================================================
#
#   Base.read!
#
# 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
# ==============================================================================

@enum(Sense, SENSE_N, SENSE_G, SENSE_L, SENSE_E, SENSE_UNKNOWN)
function Sense(s::String)
    if s == "G"
        return SENSE_G
    elseif s == "L"
        return SENSE_L
    elseif s == "E"
        return SENSE_E
    elseif s == "N"
        return SENSE_N
    else
        return SENSE_UNKNOWN
    end
end
mutable struct TempMPSModel
    name::String
    obj_name::String
    c::Vector{Float64}
    col_lower::Vector{Float64}
    col_upper::Vector{Float64}
    row_lower::Vector{Float64}
    row_upper::Vector{Float64}
    sense::Vector{Sense}
    A::Vector{Vector{Tuple{Int, Float64}}}
    vtype::Vector{VType}
    name_to_col::Dict{String, Int}
    col_to_name::Vector{String}
    name_to_row::Dict{String, Int}
    row_to_name::Vector{String}
    intorg_flag::Bool  # A flag used to parse COLUMNS section.
end

function TempMPSModel()
    return TempMPSModel(
        "",
        "",
        Float64[],  # c
        Float64[],  # col_lower
        Float64[],  # col_upper
        Float64[],  # row_lower
        Float64[],  # row_upper
        Sense[],   # sense
        Vector{Vector{Tuple{Int, Float64}}}[],  # A
        Bool[],
        Dict{String, Int}(),
        String[],
        Dict{String, Int}(),
        String[],
        false,
    )
end

@enum(
    Headers,
    HEADER_NAME,
    HEADER_ROWS,
    HEADER_COLUMNS,
    HEADER_RHS,
    HEADER_RANGES,
    HEADER_BOUNDS,
    HEADER_SOS,
    HEADER_ENDATA,
    HEADER_UNKNOWN,
)

# Headers(s) gets called _alot_ (on every line), so we try very hard to be
# efficient.]
function Headers(s::AbstractString)
    N = length(s)
    if N > 7 || N < 3
        return HEADER_UNKNOWN
    elseif N == 3
        x = first(s)
        if (x == 'R' || x == 'r') && uppercase(s) == "RHS"
            return HEADER_RHS
        elseif (x == 'S' || x == 's') && uppercase(s) == "SOS"
            return HEADER_SOS
        end
    elseif N == 4
        x = first(s)
        if (x == 'R' || x == 'r') && uppercase(s) == "ROWS"
            return HEADER_ROWS
        end
    elseif N == 5
        return HEADER_UNKNOWN
    elseif N == 6
        x = first(s)
        if (x == 'R' || x == 'r') && uppercase(s) == "RANGES"
            return HEADER_RANGES
        elseif (x == 'B' || x == 'b') && uppercase(s) == "BOUNDS"
            return HEADER_BOUNDS
        elseif (x == 'E' || x == 'e') && uppercase(s) == "ENDATA"
            return HEADER_ENDATA
        end
    elseif N == 7
        x = first(s)
        if (x == 'C' || x == 'c') && (uppercase(s) == "COLUMNS")
            return HEADER_COLUMNS
        end
    end
    return HEADER_UNKNOWN
end

function line_to_items(line)
    items = split(line, " "; keepempty = false)
    return String.(items)
end

"""
    Base.read!(io::IO, model::FileFormats.MPS.Model)

Read `io` in the MPS file format and store the result in `model`.
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end
    data = TempMPSModel()
    header = HEADER_NAME
    multi_objectives = String[]
    while !eof(io)
        line = string(strip(readline(io)))
        if isempty(line) || first(line) == '*'
            continue  # Skip blank lines and comments.
        end
        h = Headers(line)
        if h != HEADER_UNKNOWN
            header = h
            continue
        else
            # Carry on with the previous header
        end
        # TODO: split into hard fields based on column indices.
        items = line_to_items(line)
        if header == HEADER_NAME
            parse_name_line(data, items)
        elseif header == HEADER_ROWS
            multi_obj = parse_rows_line(data, items)
            if multi_obj !== nothing
                push!(multi_objectives, multi_obj)
            end
        elseif header == HEADER_COLUMNS
            parse_columns_line(data, items, multi_objectives)
        elseif header == HEADER_RHS
            parse_rhs_line(data, items)
        elseif header == HEADER_RANGES
            parse_ranges_line(data, items)
        elseif header == HEADER_BOUNDS
            parse_bounds_line(data, items)
        elseif header == HEADER_SOS
            error("TODO(odow): implement parsing of SOS constraints.")
        else
            @assert header == HEADER_ENDATA
            break
        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  # free variable
end

function copy_to(model::Model, data::TempMPSModel)
    MOI.set(model, MOI.Name(), data.name)
    variable_map = Dict{String, MOI.VariableIndex}()
    for (i, name) in enumerate(data.col_to_name)
        _add_variable(model, data, variable_map, i, name)
    end
    _add_objective(model, data, variable_map)
    for (j, c_name) in enumerate(data.row_to_name)
        set = bounds_to_set(data.row_lower[j], data.row_upper[j])
        if set !== nothing
            _add_linear_constraint(model, data, variable_map, j, c_name, set)
        end
        # `else` is a free constraint. Don't add it.
    end
    return
end

function _add_variable(model, data, variable_map, i, name)
    x = MOI.add_variable(model)
    variable_map[name] = x
    MOI.set(model, MOI.VariableName(), x, name)
    set = bounds_to_set(data.col_lower[i], data.col_upper[i])
    if set !== nothing
        MOI.add_constraint(model, MOI.SingleVariable(x), set)
    end
    if data.vtype[i] == VTYPE_INTEGER
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
    elseif data.vtype[i] == VTYPE_BINARY
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    end
    return
end

function _add_objective(model, data, variable_map)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(data.c[i], variable_map[v])
                for (i, v) in enumerate(data.col_to_name) if !iszero(data.c[i])
            ],
            0.0,
        )
    )
    return
end

function _add_linear_constraint(model, data, variable_map, j, c_name, set)
    terms = MOI.ScalarAffineTerm{Float64}[
        MOI.ScalarAffineTerm(coef, variable_map[data.col_to_name[i]])
        for (i, coef) in data.A[j]
    ]
    c = MOI.add_constraint(model, MOI.ScalarAffineFunction(terms, 0.0), set)
    MOI.set(model, MOI.ConstraintName(), c, c_name)
    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 = Sense(items[1]), items[2]
    if haskey(data.name_to_row, name)
        error("Duplicate row encountered: $(join(items, " ")).")
    elseif sense == SENSE_UNKNOWN
        error("Invalid row sense: $(join(items, " "))")
    end
    if sense == SENSE_N
        if data.obj_name != ""
            return name  # Detected a duplicate objective. Skip it.
        end
        data.obj_name = name
        return
    end
    if name == data.obj_name
        error("Found row with same name as objective: $(join(items, " ")).")
    end
    # Add some default bounds for the constraints.
    push!(data.row_to_name, name)
    row = length(data.row_to_name)
    data.name_to_row[name] = row
    push!(data.sense, sense)
    push!(data.A, Tuple{Int, Float64}[])
    if sense == SENSE_G
        push!(data.row_lower, 0.0)
        push!(data.row_upper, Inf)
        data.row_upper[row] = Inf
    elseif sense == SENSE_L
        push!(data.row_lower, -Inf)
        push!(data.row_upper, 0.0)
    else
        @assert sense == SENSE_E
        push!(data.row_lower, 0.0)
        push!(data.row_upper, 0.0)
    end
    return
end

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

function parse_single_coefficient(
    data, row_name::String, column::Int, value::Float64
)
    if row_name == data.obj_name
        data.c[column] += value
        return
    end
    row = get(data.name_to_row, row_name, nothing)
    if row === nothing
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    push!(data.A[row], (column, value))
    return
end

function _add_new_column(data, column_name)
    if haskey(data.name_to_col, column_name)
        return
    end
    push!(data.col_to_name, column_name)
    data.name_to_col[column_name] = length(data.col_to_name)
    push!(data.c, 0.0)
    push!(data.col_lower, 0.0)
    push!(data.col_upper, Inf)
    push!(data.vtype, VTYPE_CONTINUOUS)
    return
end

function _set_intorg(data, column, column_name)
     if data.intorg_flag
        data.vtype[column] = VTYPE_INTEGER
    elseif data.vtype[column] != VTYPE_CONTINUOUS
        error(
            "Variable $(column_name) appeared in COLUMNS outside an " *
            "`INT` marker after already being declared as integer."
        )
    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 row_name == "'MARKER'"
            if value == "'INTORG'"
                data.intorg_flag = true
                return
            elseif value == "'INTEND'"
                data.intorg_flag = false
                return
            end
        elseif row_name in multi_objectives
            return
        end
        _add_new_column(data, column_name)
        column = data.name_to_col[column_name]
        parse_single_coefficient(data, row_name, column, parse(Float64, value))
        _set_intorg(data, column, column_name)
    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 row_name_1 in multi_objectives || row_name_2 in multi_objectives
            return
        end
        _add_new_column(data, column_name)
        column = data.name_to_col[column_name]
        parse_single_coefficient(
            data, row_name_1, column, parse(Float64, value_1)
        )
        parse_single_coefficient(
            data, row_name_2, column, parse(Float64, value_2)
        )
        _set_intorg(data, column, column_name)
    else
        error("Malformed COLUMNS line: $(join(items, " "))")
    end
    return
end

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

function parse_single_rhs(
    data, row_name::String, value::Float64, items::Vector{String}
)
    row = get(data.name_to_row, row_name, nothing)
    if row === nothing
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    if data.sense[row] == SENSE_E
        data.row_upper[row] = value
        data.row_lower[row] = value
    elseif data.sense[row] == SENSE_G
        data.row_lower[row] = value
    elseif data.sense[row] == SENSE_L
        data.row_upper[row] = value
    else
        @assert data.sense[row] == 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, parse(Float64, value), items)
    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, parse(Float64, value_1), items)
        parse_single_rhs(data, row_name_2, parse(Float64, value_2), items)
    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::String, value::Float64)
    row = get(data.name_to_row, row_name, nothing)
    if row === nothing
        error("ROW name $(row_name) not recognised. Is it in the ROWS field?")
    end
    if data.sense[row] == SENSE_G
        data.row_upper[row] = data.row_lower[row] + abs(value)
    elseif data.sense[row] == SENSE_L
        data.row_lower[row] = data.row_upper[row] - abs(value)
    elseif data.sense[row] == SENSE_E
        if value > 0.0
            data.row_upper[row] = data.row_lower[row] + value
        else
            data.row_lower[row] = data.row_upper[row] + 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]
        _, row_name, value = items
        parse_single_range(data, row_name, parse(Float64, value))
    elseif length(items) == 5
        # [rhs name] [row name] [value] [row name 2] [value 2]
        _, row_name_1, value_1, row_name_2, value_2 = items
        parse_single_range(data, row_name_1, parse(Float64, value_1))
        parse_single_range(data, row_name_2, parse(Float64, value_2))
    else
        error("Malformed RANGES line: $(join(items, " "))")
    end
    return
end

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

function _parse_single_bound(data, column_name::String, bound_type::String)
    col = get(data.name_to_col, column_name, nothing)
    if col === nothing
        error("Column name $(column_name) not found.")
    end
    if bound_type == "PL"
        data.col_upper[col] = Inf
    elseif bound_type == "MI"
        data.col_lower[col] = -Inf
    elseif bound_type == "FR"
        data.col_lower[col] = -Inf
        data.col_upper[col] = Inf
    elseif bound_type == "BV"
        data.col_lower[col] = -Inf
        data.col_upper[col] = Inf
        data.vtype[col] = VTYPE_BINARY
    else
        error("Invalid bound type $(bound_type): $(column_name)")
    end
end

function _parse_single_bound(
    data, column_name::String, bound_type::String, value::Float64
)
    col = get(data.name_to_col, column_name, nothing)
    if col === nothing
        error("Column name $(column_name) not found.")
    end
    if bound_type == "FX"
        data.col_lower[col] = value
        data.col_upper[col] = value
    elseif bound_type == "UP"
        data.col_upper[col] = value
    elseif bound_type == "LO"
        data.col_lower[col] = value
    elseif bound_type == "LI"
        data.col_lower[col] = value
        data.vtype[col] = VTYPE_INTEGER
    elseif bound_type == "UI"
        data.col_upper[col] = value
        data.vtype[col] = VTYPE_INTEGER
    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.
        data.col_lower[col] = -Inf
        data.col_upper[col] = Inf
    elseif bound_type == "BV"
        # A similar situation happens with BV bounds in leo1 and leo2.
        data.col_lower[col] = -Inf
        data.col_upper[col] = Inf
        data.vtype[col] = VTYPE_BINARY
    else
        error("Invalid bound type $(bound_type): $(column_name)")
    end
end

function parse_bounds_line(data::TempMPSModel, items::Vector{String})
    if length(items) == 3
        bound_type, _, column_name = items
        _parse_single_bound(data, column_name, bound_type)
    elseif length(items) == 4
        bound_type, _, column_name, value = items
        _parse_single_bound(data, column_name, bound_type, parse(Float64, value))
    else
        error("Malformed BOUNDS line: $(join(items, " "))")
    end
    return
end

end


text:
Base.rand(d::AbstractMeasure) = rand(Random.GLOBAL_RNG, Float64, d)

Base.rand(T::Type, μ::AbstractMeasure) = rand(Random.GLOBAL_RNG, T, μ)

Base.rand(rng::AbstractRNG, d::AbstractMeasure) = rand(rng, Float64, d)

Base.rand(rng::AbstractRNG, T::Type, d::ParameterizedMeasure) = rand(rng, distproxy(d))

@inline Random.rand!(d::AbstractMeasure, arr::AbstractArray) = rand!(GLOBAL_RNG, d, arr)

# struct ArraySlot{A,I}
#     arr::A
#     i::I
# end

# function rand!(rng::AbstractRNG, d::AbstractMeasure, x::ArraySlot)
#     x.arr[x.i...] = rand(rng, d)
# end

ext: jl


text:
function greedy_grid(F::AbstractFiber, grid, generator::StartValueGenerator,
    initial_queue,
    bitmap=empty(grid), queued=copy(bitmap.data); callback=_do_nothing, options=MembershipTestOptions())
    # queued = copy(bitmap.data)
    queue = copy(initial_queue)

    if Threads.nthreads() > 1
        Fs = [F]
        for _=2:Threads.nthreads()
            push!(Fs, deepcopy(F))
        end
        partition = partition_work(length(queue))
        Threads.@threads for tid = 1:Threads.nthreads()
            t_queue = queue[partition[tid]]
            G = Fs[tid]
            greedy_grid_kernel!(bitmap, queued, t_queue, G, grid, generator, options, callback)
        end
    else
        greedy_grid_kernel!(bitmap, queued, queue, F, grid, generator, options, callback)
    end

    bitmap
end

function greedy_grid_kernel!(bitmap, queued, queue, F::AbstractFiber, grid, generator::StartValueGenerator, options, callback)
    while !isempty(queue)
        k = pop!(queue)
        if !bitmap[k]
            w = grid[k]
            result = membershiptest(F, w, generator, options)
            if result.successfull
                bitmap[k] = true
                callback(result, k, bitmap)
                addneighbours!(queue, queued, k)
            end
        end
    end
end

_do_nothing(result, k, bitmap) = nothing

function addneighbours!(queue, queued::BitMatrix, ij)
    i, j = Tuple(ij)
    i0 = max(i - 1, 1)
    i1 = min(i + 1, size(queued, 1))

    j0 = max(j - 1, 1)
    j1 = min(j + 1, size(queued, 2))
    for jj = j0:j1, ii = i0:i1
        if !queued[ii, jj]
            queued[ii, jj] = true
            push!(queue,  CartesianIndex(ii, jj))
        end
    end
end

function addneighbours!(queue, queued::BitArray{3}, ijk)
    i, j, k = Tuple(ijk)
    i0 = max(i - 1, 1)
    i1 = min(i + 1, size(queued, 1))

    j0 = max(j - 1, 1)
    j1 = min(j + 1, size(queued, 2))

    k0 = max(k - 1, 1)
    k1 = min(k + 1, size(queued, 3))
    for kk=k0:k1, jj = j0:j1, ii = i0:i1
        if !queued[ii, jj, kk]
            queued[ii, jj, kk] = true
            push!(queue, CartesianIndex(ii, jj, kk))
        end
    end
end

function initial_neighbour_queue(B::BitArray{N}) where N
    queued = copy(B)
    initial_queue = CartesianIndex{N}[]
    for (ind, val) in pairs(B)
        if val
            addneighbours!(initial_queue, queued, ind)
        end
    end
    initial_queue, queued
end

# compared to `greedy_grid ` does this uses previous neighboring startvalues.
function greedy_grid_memorized(F::AbstractFiber, grid, generator::StartValueGenerator,
    initial_queue::Vector,
    startvalues::Array,
    bitmap=empty(grid), queued=copy(bitmap.data); options=MembershipTestOptions())
    # queued = copy(bitmap.data)

    Fs = [F]
    for _=2:Threads.nthreads()
        push!(Fs, deepcopy(F))
    end

    partition = partition_work(length(initial_queue))
    Threads.@threads for tid = 1:Threads.nthreads()
        t_queue = initial_queue[partition[tid]]
        G = Fs[tid]
        greedy_grid_memorized_kernel!(bitmap, queued, t_queue, startvalues, G, grid, generator, options)
    end

    bitmap
end
function greedy_grid_memorized_kernel!(bitmap, queued, queue, startvalues, F::AbstractFiber, grid, generator::StartValueGenerator, options)
    tmp_neighbor_startvalues = Vector{eltype(startvalues)}()
    while !isempty(queue)
        k = pop!(queue)
        # @show length(queue)
        if !bitmap[k]
            w = grid[k]
            neighbor_startvalues!(tmp_neighbor_startvalues, F, startvalues, bitmap, k)
            result = membershiptest(F, w, generator, tmp_neighbor_startvalues, options)
            if result.successfull
                startvalues[k] = result.solution
                bitmap[k] = true
                addneighbours!(queue, queued, bitmap, k)
            end
        end
    end
end

function addneighbours!(queue, queued::BitMatrix, bitmap::Bitmap2D, k)
    i, j = Tuple(k)
    i0 = max(i - 1, 1)
    i1 = min(i + 1, size(queued, 1))

    j0 = max(j - 1, 1)
    j1 = min(j + 1, size(queued, 2))
    for ii = i0:i1, jj = j0:j1
        if !queued[ii, jj] || !bitmap[ii, jj]
            queued[ii, jj] = true
            push!(queue, CartesianIndex(ii, jj))
        end
    end
end

function addneighbours!(queue, queued::BitArray{3}, bitmap::Bitmap3D, l)
    i, j, k = Tuple(l)
    i0 = max(i - 1, 1)
    i1 = min(i + 1, size(queued, 1))

    j0 = max(j - 1, 1)
    j1 = min(j + 1, size(queued, 2))

    k0 = max(k - 1, 1)
    k1 = min(k + 1, size(queued, 3))
    for ii = i0:i1, jj = j0:j1, kk=k0:k1
        if !queued[ii, jj, kk] || !bitmap[ii, jj, kk]
            queued[ii, jj, kk] = true
            push!(queue, CartesianIndex(ii, jj, kk))
        end
    end
end
function neighbor_startvalues!(out, F, startvalues, B::Bitmap2D, k)
    empty!(out)
    i, j = Tuple(k)
    offset = 1
    for jj in j-offset:j+offset, ii in i-offset:i+offset
        if safe_getindex(B, ii, jj) && !isnan(startvalues[ii, jj])
            predict_startvalue!(out, F, startvalues[ii, jj], B.grid[ii, jj], B.grid[k])
        end
    end
end

function neighbor_startvalues!(out, F, startvalues, B::Bitmap3D, index)
    empty!(out)
    i, j, k = Tuple(index)
    offset = 1
    for kk in k-offset:k+offset, jj in j-offset:j+offset, ii in i-offset:i+offset
        if safe_getindex(B, ii, jj, kk) && !isnan(startvalues[ii, jj, kk])
            predict_startvalue!(out, F, startvalues[ii, jj, kk], B.grid[ii, jj, kk], B.grid[index])
        end
    end
    out
end

function predict_startvalue!(out, F, x_known, y_known, y_target)
    push!(out, x_known)
end
function predict_startvalue!(out, F::ImaginaryFiber2D, x_known, y_known, y_target)
    xy = SVector(x_known[1], x_known[2], y_known[1], y_known[2])
    U = F.U
    @inbounds begin
        ∇f_re = SP.gradient(F.f_re, xy)
        for i=1:4
            U[1, i] = ∇f_re[i]
        end
        U[2, 3] = U[1, 1]
        U[2, 4] = U[1, 2]
        U[2, 1] = -U[1, 3]
        U[2, 2] = -U[1, 4]
    end

    ∂xy = SMatrix{2,2}(U[1,1], U[2,1], U[1, 2], U[2, 2])
    ∂t = SMatrix{2,2}(U[1,3], U[2,3], U[1, 4], U[2, 4]) *
        SVector(y_target[1] - y_known[1], y_target[2] - y_known[2])

    Δx, ok = mysolve(∂xy, ∂t)
    if ok
        push!(out, x_known - Δx)
        # push!(out, x_known - 0.5Δx)
        push!(out, x_known)

    else
        @show "not ok"
        push!(out, x_known)
    end
end


text:
using Distributions
using Plots
using PDMats
using BenchmarkTools
using LinearAlgebra
using StatsPlots
using DataFrames
using RCall

include("probit.reg.jl")
include("designs.jl")

# Size of full domain
N = 100

# Size of initial sample
n = 10

# Parameters
β = [0.0, 0.5]

# Data
X = hcat(fill(1.0, N), sort(rand(Normal(0, 3), N)))
p = cdf.(Normal(), X * β)
y = rand.(Bernoulli.(p))

# Selecting observations
obs = sample(1:N, n, replace = false)

# Fitting model
chain = probit_reg(y[obs], X[obs, :], X, 500000)

# Optimal design ==========================================================

# Multiple imputation
ppred = mean(chain.y, dims = 2)
ypred = hcat(fill(0, N), fill(1, N))

# Designs
D = findall(x -> !(x in obs), 1:N)
L = length(D)

# PPRB with post-hoc MI
scores_pprb = designs_pprb(ypred, D, chain.p, chain.β)

vars = sum(scores_pprb.vars, dims = 1)[1, :, :]
meanvar = ppred[D] .* vars[:, 2] + (1.0 .- ppred[D]) .* vars[:, 1]

# Output data for plotting in R ================================================

data = DataFrame(loc = 1:100, y = y, x = X[:, 2], K1 = 0, score = NaN)
data[obs, :K1] = 1
data[D, :score] = meanvar

R"
saveRDS($(chain.p), 'output/chain.rds')
saveRDS($(data), 'output/data.rds')
"

max_stars_repo_name: clint-leach/optimal-design

max_issues_repo_name: clint-leach/optimal-design

max_forks_repo_name: clint-leach/optimal-design


text:
using PersistentArrays
using Base.Test

default_array = PersistentArray(4)
default_array[2] = 3.5

@test version(default_array) == 1

commit!(default_array)
@test version(default_array) == 2

@test default_array[1] == 0.0
@test default_array[2] == 3.5

default_array[2] = 2.5
commit!(default_array)
@test version(default_array) == 3

@test default_array[:] == [0.0, 2.5, 0.0, 0.0]
@test lookup(default_array, 1, :) == [0.0, 3.5, 0.0, 0.0]

# Test that we can still overwrite previous version of indices if necessary
update!(default_array, 4.5, 1, 2)
@test lookup(default_array, 1, :) == [0.0, 4.5, 0.0, 0.0]

# Run show
show(default_array)

# Check that everying still works for more than 1 dimension
ndim_array = PersistentArray(3,2)
ndim_array[1,1] = 3.5
commit!(ndim_array)

tmp_array = zeros(3,2)
tmp_array[1,1] = 3.5
@test ndim_array[:,:] == tmp_array

commit!(ndim_array)

# Check that regular array assignments also still work
tmp_array[1,:] = ones(2)
ndim_array[1,:] = ones(2)

@test ndim_array[:,:] == tmp_array
@test lookup(ndim_array, 1, :,:) != tmp_array[:,:]

# Check that size still works
@test size(ndim_array) == (3,2)
@test size(ndim_array, 1) == 3

println(ndim_array)

# Check that sub works
sub_array = sub(ndim_array, 1:2, :)
sub_tmp_array = sub(tmp_array, 1:2, :)
println(sub_array)
println(sub_tmp_array)
@test sub_array[1,1] == sub_tmp_array[1,1]


text:
# Find a Position in a Genome Minimizing the Skew

function get_skewness_argmins(genome)
    skew = 0
    min_skew_value = 0
    min_skew_indices = [0]

    for (ix, nucleotide) in enumerate(genome)
        if nucleotide == 'C'
            skew -= 1
        elseif nucleotide == 'G'
            skew += 1
        end

        if skew == min_skew_value
            push!(min_skew_indices, ix)
        elseif skew < min_skew_value
            min_skew_value = skew
            min_skew_indices = [ix]
        end
    end
    min_skew_indices
end

function main()
    genome = readline()
    result = get_skewness_argmins(genome)
    println(join(result, " "))
end

main()


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/opt.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_set_string3,
    av_set_double,
    av_set_q,
    av_set_int,
    av_get_double,
    av_get_q,
    av_get_int,
    av_get_string,
    av_next_option,
    av_opt_show2,
    av_opt_set_defaults,
    av_opt_set_defaults2,
    av_set_options_string,
    av_opt_set_from_string,
    av_opt_free,
    av_opt_flag_is_set,
    av_opt_set_dict,
    av_opt_set_dict2,
    av_opt_get_key_value,
    av_opt_eval_flags,
    av_opt_eval_int,
    av_opt_eval_int64,
    av_opt_eval_float,
    av_opt_eval_double,
    av_opt_eval_q,
    av_opt_find,
    av_opt_find2,
    av_opt_next,
    av_opt_child_next,
    av_opt_child_class_next,
    av_opt_set,
    av_opt_set_int,
    av_opt_set_double,
    av_opt_set_q,
    av_opt_set_bin,
    av_opt_set_image_size,
    av_opt_set_pixel_fmt,
    av_opt_set_sample_fmt,
    av_opt_set_video_rate,
    av_opt_set_channel_layout,
    av_opt_set_dict_val,
    av_opt_get,
    av_opt_get_int,
    av_opt_get_double,
    av_opt_get_q,
    av_opt_get_image_size,
    av_opt_get_pixel_fmt,
    av_opt_get_sample_fmt,
    av_opt_get_video_rate,
    av_opt_get_channel_layout,
    av_opt_get_dict_val,
    av_opt_ptr,
    av_opt_freep_ranges,
    av_opt_query_ranges,
    av_opt_copy,
    av_opt_query_ranges_default

function av_set_string3(obj,name,val,alloc::Integer,o_out)
    ccall((:av_set_string3,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint,Ptr{Ptr{AVOption}}),obj,name,val,alloc,o_out)
end

function av_set_double(obj,name,n::Cdouble)
    ccall((:av_set_double,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{Uint8},Cdouble),obj,name,n)
end

function av_set_q(obj,name,n::AVRational)
    ccall((:av_set_q,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{Uint8},AVRational),obj,name,n)
end

function av_set_int(obj,name,n::Int64)
    ccall((:av_set_int,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{Uint8},Int64),obj,name,n)
end

function av_get_double(obj,name,o_out)
    ccall((:av_get_double,libavutil),Cdouble,(Ptr{Void},Ptr{Uint8},Ptr{Ptr{AVOption}}),obj,name,o_out)
end

function av_get_q(obj,name,o_out)
    ccall((:av_get_q,libavutil),AVRational,(Ptr{Void},Ptr{Uint8},Ptr{Ptr{AVOption}}),obj,name,o_out)
end

function av_get_int(obj,name,o_out)
    ccall((:av_get_int,libavutil),Int64,(Ptr{Void},Ptr{Uint8},Ptr{Ptr{AVOption}}),obj,name,o_out)
end

function av_get_string(obj,name,o_out,buf,buf_len::Integer)
    ccall((:av_get_string,libavutil),Ptr{Uint8},(Ptr{Void},Ptr{Uint8},Ptr{Ptr{AVOption}},Ptr{Uint8},Cint),obj,name,o_out,buf,buf_len)
end

function av_next_option(obj,last)
    ccall((:av_next_option,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{AVOption}),obj,last)
end

function av_opt_show2(obj,av_log_obj,req_flags::Integer,rej_flags::Integer)
    ccall((:av_opt_show2,libavutil),Cint,(Ptr{Void},Ptr{Void},Cint,Cint),obj,av_log_obj,req_flags,rej_flags)
end

function av_opt_set_defaults(s)
    ccall((:av_opt_set_defaults,libavutil),Void,(Ptr{Void},),s)
end

function av_opt_set_defaults2(s,mask::Integer,flags::Integer)
    ccall((:av_opt_set_defaults2,libavutil),Void,(Ptr{Void},Cint,Cint),s,mask,flags)
end

function av_set_options_string(ctx,opts,key_val_sep,pairs_sep)
    ccall((:av_set_options_string,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Ptr{Uint8}),ctx,opts,key_val_sep,pairs_sep)
end

function av_opt_set_from_string(ctx,opts,shorthand,key_val_sep,pairs_sep)
    ccall((:av_opt_set_from_string,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Ptr{Uint8}},Ptr{Uint8},Ptr{Uint8}),ctx,opts,shorthand,key_val_sep,pairs_sep)
end

function av_opt_free(obj)
    ccall((:av_opt_free,libavutil),Void,(Ptr{Void},),obj)
end

function av_opt_flag_is_set(obj,field_name,flag_name)
    ccall((:av_opt_flag_is_set,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8}),obj,field_name,flag_name)
end

function av_opt_set_dict(obj,options)
    ccall((:av_opt_set_dict,libavutil),Cint,(Ptr{Void},Ptr{Ptr{AVDictionary}}),obj,options)
end

function av_opt_set_dict2(obj,options,search_flags::Integer)
    ccall((:av_opt_set_dict2,libavutil),Cint,(Ptr{Void},Ptr{Ptr{AVDictionary}},Cint),obj,options,search_flags)
end

function av_opt_get_key_value(ropts,key_val_sep,pairs_sep,flags::Integer,rkey,rval)
    ccall((:av_opt_get_key_value,libavutil),Cint,(Ptr{Ptr{Uint8}},Ptr{Uint8},Ptr{Uint8},Uint32,Ptr{Ptr{Uint8}},Ptr{Ptr{Uint8}}),ropts,key_val_sep,pairs_sep,flags,rkey,rval)
end

function av_opt_eval_flags(obj,o,val,flags_out)
    ccall((:av_opt_eval_flags,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{Cint}),obj,o,val,flags_out)
end

function av_opt_eval_int(obj,o,val,int_out)
    ccall((:av_opt_eval_int,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{Cint}),obj,o,val,int_out)
end

function av_opt_eval_int64(obj,o,val,int64_out)
    ccall((:av_opt_eval_int64,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{Int64}),obj,o,val,int64_out)
end

function av_opt_eval_float(obj,o,val,float_out)
    ccall((:av_opt_eval_float,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{Cfloat}),obj,o,val,float_out)
end

function av_opt_eval_double(obj,o,val,double_out)
    ccall((:av_opt_eval_double,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{Cdouble}),obj,o,val,double_out)
end

function av_opt_eval_q(obj,o,val,q_out)
    ccall((:av_opt_eval_q,libavutil),Cint,(Ptr{Void},Ptr{AVOption},Ptr{Uint8},Ptr{AVRational}),obj,o,val,q_out)
end

function av_opt_find(obj,name,unit,opt_flags::Integer,search_flags::Integer)
    ccall((:av_opt_find,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint,Cint),obj,name,unit,opt_flags,search_flags)
end

function av_opt_find2(obj,name,unit,opt_flags::Integer,search_flags::Integer,target_obj)
    ccall((:av_opt_find2,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint,Cint,Ptr{Ptr{Void}}),obj,name,unit,opt_flags,search_flags,target_obj)
end

function av_opt_next(obj,prev)
    ccall((:av_opt_next,libavutil),Ptr{AVOption},(Ptr{Void},Ptr{AVOption}),obj,prev)
end

function av_opt_child_next(obj,prev)
    ccall((:av_opt_child_next,libavutil),Ptr{Void},(Ptr{Void},Ptr{Void}),obj,prev)
end

function av_opt_child_class_next(parent,prev)
    ccall((:av_opt_child_class_next,libavutil),Ptr{AVClass},(Ptr{AVClass},Ptr{AVClass}),parent,prev)
end

function av_opt_set(obj,name,val,search_flags::Integer)
    ccall((:av_opt_set,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint),obj,name,val,search_flags)
end

function av_opt_set_int(obj,name,val::Int64,search_flags::Integer)
    ccall((:av_opt_set_int,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Int64,Cint),obj,name,val,search_flags)
end

function av_opt_set_double(obj,name,val::Cdouble,search_flags::Integer)
    ccall((:av_opt_set_double,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cdouble,Cint),obj,name,val,search_flags)
end

function av_opt_set_q(obj,name,val::AVRational,search_flags::Integer)
    ccall((:av_opt_set_q,libavutil),Cint,(Ptr{Void},Ptr{Uint8},AVRational,Cint),obj,name,val,search_flags)
end

function av_opt_set_bin(obj,name,val,size::Integer,search_flags::Integer)
    ccall((:av_opt_set_bin,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint,Cint),obj,name,val,size,search_flags)
end

function av_opt_set_image_size(obj,name,w::Integer,h::Integer,search_flags::Integer)
    ccall((:av_opt_set_image_size,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Cint,Cint),obj,name,w,h,search_flags)
end

function av_opt_set_pixel_fmt(obj,name,fmt::AVPixelFormat,search_flags::Integer)
    ccall((:av_opt_set_pixel_fmt,libavutil),Cint,(Ptr{Void},Ptr{Uint8},AVPixelFormat,Cint),obj,name,fmt,search_flags)
end

function av_opt_set_sample_fmt(obj,name,fmt::AVSampleFormat,search_flags::Integer)
    ccall((:av_opt_set_sample_fmt,libavutil),Cint,(Ptr{Void},Ptr{Uint8},AVSampleFormat,Cint),obj,name,fmt,search_flags)
end

function av_opt_set_video_rate(obj,name,val::AVRational,search_flags::Integer)
    ccall((:av_opt_set_video_rate,libavutil),Cint,(Ptr{Void},Ptr{Uint8},AVRational,Cint),obj,name,val,search_flags)
end

function av_opt_set_channel_layout(obj,name,ch_layout::Int64,search_flags::Integer)
    ccall((:av_opt_set_channel_layout,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Int64,Cint),obj,name,ch_layout,search_flags)
end

function av_opt_set_dict_val(obj,name,val,search_flags::Integer)
    ccall((:av_opt_set_dict_val,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Ptr{AVDictionary},Cint),obj,name,val,search_flags)
end

function av_opt_get(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Ptr{Uint8}}),obj,name,search_flags,out_val)
end

function av_opt_get_int(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get_int,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Int64}),obj,name,search_flags,out_val)
end

function av_opt_get_double(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get_double,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Cdouble}),obj,name,search_flags,out_val)
end

function av_opt_get_q(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get_q,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{AVRational}),obj,name,search_flags,out_val)
end

function av_opt_get_image_size(obj,name,search_flags::Integer,w_out,h_out)
    ccall((:av_opt_get_image_size,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Cint},Ptr{Cint}),obj,name,search_flags,w_out,h_out)
end

function av_opt_get_pixel_fmt(obj,name,search_flags::Integer,out_fmt)
    ccall((:av_opt_get_pixel_fmt,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{AVPixelFormat}),obj,name,search_flags,out_fmt)
end

function av_opt_get_sample_fmt(obj,name,search_flags::Integer,out_fmt)
    ccall((:av_opt_get_sample_fmt,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{AVSampleFormat}),obj,name,search_flags,out_fmt)
end

function av_opt_get_video_rate(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get_video_rate,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{AVRational}),obj,name,search_flags,out_val)
end

function av_opt_get_channel_layout(obj,name,search_flags::Integer,ch_layout)
    ccall((:av_opt_get_channel_layout,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Int64}),obj,name,search_flags,ch_layout)
end

function av_opt_get_dict_val(obj,name,search_flags::Integer,out_val)
    ccall((:av_opt_get_dict_val,libavutil),Cint,(Ptr{Void},Ptr{Uint8},Cint,Ptr{Ptr{AVDictionary}}),obj,name,search_flags,out_val)
end

function av_opt_ptr(avclass,obj,name)
    ccall((:av_opt_ptr,libavutil),Ptr{Void},(Ptr{AVClass},Ptr{Void},Ptr{Uint8}),avclass,obj,name)
end

function av_opt_freep_ranges(ranges)
    ccall((:av_opt_freep_ranges,libavutil),Void,(Ptr{Ptr{AVOptionRanges}},),ranges)
end

function av_opt_query_ranges(arg1,obj,key,flags::Integer)
    ccall((:av_opt_query_ranges,libavutil),Cint,(Ptr{Ptr{AVOptionRanges}},Ptr{Void},Ptr{Uint8},Cint),arg1,obj,key,flags)
end

function av_opt_copy(dest,src)
    ccall((:av_opt_copy,libavutil),Cint,(Ptr{Void},Ptr{Void}),dest,src)
end

function av_opt_query_ranges_default(arg1,obj,key,flags::Integer)
    ccall((:av_opt_query_ranges_default,libavutil),Cint,(Ptr{Ptr{AVOptionRanges}},Ptr{Void},Ptr{Uint8},Cint),arg1,obj,key,flags)
end


text:
export Simplex

struct Consensus{T <: Real}

end

function project!(s::Consensus{<:Real}, v::Matrix{T}, y::Matrix{T}) where {T <: Real}
  return repeat(mean(Y, dims = 2), inner = [1, size(Y, 2)])
end

function project(s::Consensus{<:Real}, y::Matrix{<:Real})
    project!(s, similar(y), y)
end


text:
include(joinpath(@__DIR__, "config.jl"))
include(joinpath(@__DIR__, "..", "subspace_util.jl"))

using OneClassActiveLearning, SVDD, JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

data_dirs = ["HeartDisease", "Stamps"]
data_info = "01-v-comparison"

classify_precision = SVDD.OPT_PRECISION

### learning scenario ###
initial_pool_strategy = [("Pu", Dict())]
split_strategy = [("Sf", Dict())]

#### models ####

steps = collect(10.0.^range(-3, stop=-0.05, length=10))

models = []
for v_out in steps
    for v_in in (steps.^-1)
        push!(models, Dict(:type => :SubSVDD,
                           :param => Dict{Symbol, Any}(
                                      :weight_update_strategy => SVDD.FixedWeightStrategy(v_in, v_out))))
    end
end

init_strategies = [SimpleSubspaceStrategy(SVDD.WangGammaStrategy(solver, range(0.1, stop=20, length=10), 1.0), FixedCStrategy(0.45), gamma_scope=Val(:Subspace))]

subspace_generators = [generate_randsub_from_data(10; min_ss_size=5, max_ss_size=8)]

#### oracle ####
oracle_param = Dict{Symbol, Any}(
    :type => PoolOracle,
    :param => Dict{Symbol, Any}()
)

#### query strategies ####
query_strategies = [
    Dict(:type => SubspaceQs{DecisionBoundaryPQs},
           :param => Dict{Symbol, Any}(:scale_fct => min_max_normalize,
                                       :combination_fct => + ))]

num_al_iterations = 50
num_resamples_initial_pool = 1
num_resamples_subspaces = 3
data_output_root = joinpath(data_root, "output", "evaluation-part2")

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
#!/usr/bin/env julia

using StochasticCubicDescent
using Base.Test

@testset "gradient descent on a cubic problem" begin
    # define a test function
    H0 = [1. 2.; 3. -4.]
    A = H0' * H0 - 10 * eye(2)  # the eigenvalues of A are -6.18 and 16.18
    b = [-1.6, -2.8]
    ρ = 3.3
    val(x) = 0.5x' * A * x + b' * x + (ρ / 3) * norm(x)^3
    grad!(x, grad) = grad[:] = A * x + b + ρ * norm(x) * x
    cubic_objective = ObjectiveFunction(val, grad!)

    # verify that a gradient is correct
    x = [5.0, 7.0]
    tmp = zeros(2)
    finite_diff_grad_sum = (val(x + 1e-4) - val(x)) / 1e-4
    analytic_grad_sum = sum(grad!(x, tmp))
    @test finite_diff_grad_sum ≈ analytic_grad_sum atol=1e-1

    # Below the step size and initial iterate are set according to
    #     "Gradient Descent Efficiently Finds the Cubic-Rgularized Non-Convex
    #      Newton Step" by Yair Carmon and John Duchi. 2017.
    β = norm(A)
    R = β / (2ρ) + sqrt((β / 2ρ)^2 + norm(b) / ρ)
    η = (8(β + ρ * R))^(-1)

    Rc_lower_bound = -β / ρ + R
    x0 = -Rc_lower_bound / norm(b) * b

    solver = GradientDescent(η, 1000, true, GDState[])
    state = GDState(x0)
    optimize!(solver, state, cubic_objective)

    @test norm(state.grad) ≈ 0.0 atol=1e-8

    # should hold by Lemma 2.2
    iterate_norms = [norm(st.x) for st in solver.history]
    iterate_increases = iterate_norms[2:end] - iterate_norms[1:end-1]
    @test all(iterate_increases .> 0)
end

@testset "Stochastic Cubic Descent" begin
    val(x, ξ) = 0.1x[1]^2 + exp(x[1] + 1) + x[1] * ξ[1]
    grad!(x, ξ, grad) = begin
        grad[1] = 0.2x[1] + exp(x[1] + 1) + ξ[1]
    end
    hv!(x, ξ, v, hv) = begin
        @assert x[1] < 5
        hv[1] = (0.2 + exp(x[1] + 1)) * v[1]
    end
    # very low noise
    ξ!(x, ξ) = ξ[:] = (rand(1) - 0.5) * 1e-4
    # exp(5) isn't a global Lipschitz constant, but it should hold
    # at every point that the algorithm visits
    f = ObjectiveFunction(val, grad!, hv!, ξ!, exp(5))

    solver = SCD(100)
    state = SCDState([-3.])

    optimize!(solver, state, f)

    @test norm(state.grad) ≈ 0.0 atol=1e-2
end


text:
"Evaluates a rational filter, exploiting numerical optimizations provided by JULIA."
function phi(x::Float64, z, alpha)
    y = zero(Complex{Float64});
    @inbounds @fastmath for i = eachindex(z) # @fastmath can be dangerous but no impact has been noticeable in various tests
        y = y + alpha[i] / (x[1] - z[i]) + conj(alpha[i]) / (x[1] - conj(z[i])) - alpha[i] / (x[1] + z[i]) - conj(alpha[i]) / (x[1] + conj(z[i]));
    end
    return real(y)
end

# For compatibility with JULIA versions => 0.7
function Linspace(left, right, points)
    return left:(right-left)/(points-1):right
end

"""
    wcr( z, a, G; interval_fevals=500000, upper=10., middle=2.)

Obtain the WCR of a rational filter.

Samples thousands points, exploiting numerical optimizations provided by JULIA
in order to obtain the WCR.

In terms of both computational time and, most importantly, validity, the employed
approach is the only feasible choice, as opposed to sampling the vector-based phi()
function or using the eigenvalues of the companion matrix of the nominator of the
rational filter or of its derivative, whilst both Float64 and BigFloat were tested
but could not overcome the occuring large condition numbers. This has not produced
any large derivations from the expected WCR on a set of 60 'difficult' rational
filters. More such filters are to be obtained throghout the subsequent tests in
filter computation.

Underlying assumptions: 1. Rational filters are continuous, thus limited in terms
of derivation. 2. We study at worst G=0.99998, thus we want to be able to consider
lenghts less than 0.99998^-1-1, which requires at least (0.99998^-1-1)^-1 = 49999
evaluations in (1,2)."""
function wcr(z,a,G; interval_fevals=500000, upper=10., middle=2., lowerGap = G, mintol = Inf)
    # Compute maximizer outside (-1,1) with respect to the gap
    # Proceed in three steps, contributing different precision to each interval
    maximizer = 0
    for x in Linspace(G^-1,1.1+(G^-1-1),interval_fevals)
        w = abs(phi(x,z,a))
        if (w > maximizer)
            maximizer = w
        end
    end
    for x in Linspace(1.1+(G^-1-1),middle+(G^-1-1),interval_fevals)
        w = abs(phi(x,z,a))
        if (w > maximizer)
            maximizer = w
        end
    end
    for x in Linspace(middle+(G^-1-1),upper+(G^-1-1),interval_fevals)
        w = abs(phi(x,z,a))
        if (w > maximizer)
            maximizer = w
        end
    end

    # Compute the minimizer inside (-1,1)
    minimizer = mintol
    for x in Linspace(0,lowerGap,interval_fevals)
        w = abs(phi(x,z,a))
        if (w < minimizer)
            minimizer = w
        end
    end

    return maximizer/minimizer
end


text:
"""
    BitCount

type used to hold bit counts and shift amounts
"""
const BitCount = Int16

"""
    BitField

The structured BitField.

fields:
    nbits
    shift
    maskof1s
    maskof0s
    name::Symbol
"""
struct BitField{U} <: AbstractBitField
    nbits::BitCount
    shift::BitCount
    maskof1s::U
    maskof0s::U
    name::Symbol
end

Base.eltype(::Type{BitField{U}}) where {U} = U
Base.eltype(x::BitField{U}) where {U} = U

name(x::BitField{U}) where {U<:UBits} = x.name
nbits(x::BitField{U}) where {U<:UBits} = x.nbits
shift(x::BitField{U}) where {U<:UBits} = x.shift
maskof1s(x::BitField{U}) where {U<:UBits} = x.maskof1s
maskof0s(x::BitField{U}) where {U<:UBits} = x.maskof0s

const BitField8 = BitField{UInt8}
const BitField16 = BitField{UInt16}
const BitField32 = BitField{UInt32}
const BitField64 = BitField{UInt64}
const BitField128 = BitField{UInt128}

BitField8(; name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt8, bitspan, bitshift, name)
BitField16(; name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt16, bitspan, bitshift, name)
BitField32(; name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt32, bitspan, bitshift, name)
BitField64(; name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt64, bitspan, bitshift, name)
BitField128(; name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt128, bitspan, bitshift, name)

function BitField(::Type{U}, bitspan::Int, bitshift::Int, name::Symbol) where {U<:UBits}
    span  = BitCount(bitspan)
    shift = BitCount(bitshift)
    validate(U, span, shift)

    maskof1s = onebits(U, span, shift)
    maskof0s = ~maskof1s
    return BitField{U}(span, shift, maskof1s, maskof0s, name)
end

BitField(::Type{U}, name::Symbol, bitspan::Int, bitshift::Int) where {U<:UBits} =
    BitField(U, bitspan, bitshift, name)

BitField(bitspan::Int, bitshift::Int, name::Symbol) = BitField(UInt64, bitspan, bitshift, name)

BitField(name::Symbol, bitspan::Int, bitshift::Int) = BitField(UInt64, bitspan, bitshift, name)

@inline function isolate(bitfield::BitField{U}, source::U) where {U<:UBits}
    return source & bitfield.maskof1s
end

@inline function filter(bitfield::BitField{U}, source::U) where {U<:UBits}
    return source & bitfield.maskof0s
end

@inline function Base.get(bitfield::BitField{U}, source::U) where {U<:UBits}
    return isolate(bitfield, source) >> bitfield.shift
end

@inline function Base.get(bitfield::BitField{U}, source::Base.RefValue{U}) where {U<:UBits}
    return isolate(bitfield, source[]) >> bitfield.shift
end

@inline function Base.get(bitfield::BitField{U}, source::ByRef{S,U}) where {S,U<:UBits}
    return isolate(bitfield, source.ref[]) >> bitfield.shift
end

@inline function put(bitfield::BitField{U}, value::U) where {U<:UBits}
    value << bitfield.shift
end

set(bitfield::BitField{U}, value::B, target::U) where {U<:UBits, B<:Union{UBits,IBits}} =
    set(bitfield, value%U, target)

set!(bitfield::BitField{U}, value::B, target::Base.RefValue{U})  where {U<:UBits, B<:Union{UBits,IBits}} =
    set!(bitfield, value%U, target)

set!(bitfield::BitField{U}, value::B, target::ByRef{S,U})  where {S, U<:UBits, B<:Union{UBits,IBits}} =
    set!(bitfield, value%U, target)

@inline function set(bitfield::BitField{U}, value::U, target::U)  where {U<:UBits}
    return filter(bitfield, target) | put(bitfield, value)
end

@inline function set!(bitfield::BitField{U}, value::U, target::Base.RefValue{U})  where {U<:UBits}
    target[] = set(bitfield, value, target[])
    return target
end

@inline function set!(bitfield::BitField{U}, value::U, target::ByRef{S, U})  where {S, U<:UBits}
    target.ref[] = set(bitfield, value, target.ref[])
    return target
end

function validate(::Type{U}, bitspan::Integer, bitshift::Integer) where {U<:UBits}
    if !(bitspan > 0 && bitspan + bitshift <= bitsof(U))
        if bitspan <= 0
            throw(ErrorException("bitspan ($bitspan) must be > 0"))
        else
            throw(ErrorException("bitspan + bitshift > totalbits: $bitspan + $bitshift > $(bitsizeof(U))"))
        end
    end
    return nothing
end

onebits(::Type{T}, bitspan::I, bitshift::I) where {T<:UBits, I<:Integer} =
    ~(~zero(T) << bitspan) << bitshift


text:
using DataStructures

CI = CartesianIndex

slurp_grid(filename) = parse.(Int, permutedims(hcat(collect.(readlines(filename))...)))

get_neighbors(idx) = [idx + Δ for Δ in [CI(0, 1), CI(1, 0), CI(-1, 0), CI(0, -1)]]

function in_bounds(grid, idx)
    y_max, x_max = size(grid)
    y, x = Tuple(idx)
    1 <= y <= y_max && 1 <= x <= x_max
end

get_neighbors(grid, idx) = [idx for idx in get_neighbors(idx) if in_bounds(grid, idx)]

function a_star(start, neighbors_fn, heuristic, step_cost)
    frontier = PriorityQueue()
    enqueue!(frontier, start, heuristic(start))
    previous = Dict{CI,Union{Nothing,CI}}(start => nothing)
    path_cost = Dict(start => 0)
    path = function (idx)
        result = []
        while idx !== nothing
            push!(result, idx)
            idx = previous[idx]
        end
        reverse(result)
    end

    while !isempty(frontier)
        next = dequeue!(frontier)
        if heuristic(next) == 0
            return path(next)
        end
        for idx in neighbors_fn(next)
            g = path_cost[next] + step_cost(next, idx)
            if !haskey(path_cost, idx) || g < path_cost[idx]
                frontier[idx] = g + heuristic(idx)
                path_cost[idx] = g
                previous[idx] = next
            end
        end
    end
end

function a_star_path(grid, start, goal)
    neighbors_fn(idx) = get_neighbors(grid, idx)
    heuristic(idx) = sum(Tuple(goal - idx))
    step_cost(from, to) = grid[to]
    a_star(start, neighbors_fn, heuristic, step_cost)
end

plus_mod(n) = x -> mod1(x + n, 9)

expand_map(grid) = [
    plus_mod(0).(grid) plus_mod(1).(grid) plus_mod(2).(grid) plus_mod(3).(grid) plus_mod(4).(grid)
    plus_mod(1).(grid) plus_mod(2).(grid) plus_mod(3).(grid) plus_mod(4).(grid) plus_mod(5).(grid)
    plus_mod(2).(grid) plus_mod(3).(grid) plus_mod(4).(grid) plus_mod(5).(grid) plus_mod(6).(grid)
    plus_mod(3).(grid) plus_mod(4).(grid) plus_mod(5).(grid) plus_mod(6).(grid) plus_mod(7).(grid)
    plus_mod(4).(grid) plus_mod(5).(grid) plus_mod(6).(grid) plus_mod(7).(grid) plus_mod(8).(grid)
]

function solve2(filename)
    grid = expand_map(slurp_grid(filename))
    CIs = CartesianIndices(grid)
    path = a_star_path(grid, first(CIs), last(CIs))
    reduce(+, [grid[idx] for idx in path if idx != first(CIs)])
end

solve2("15eg.txt") == 315
solve2("15.txt")


text:
# mass:        g / mol
# length:      nm
# time:        ps
# momentum:    g nm / ps mol
# energy:      kJ / mol
# temperature: K

const to_kilo = 1e-3
const to_pico = 1e12

# 2018 CODATA recommended values from NIST.
const N_A = 6.022_140_76e23 # 1 / mol
const planck = 6.626_070_15e-34 * to_kilo * to_pico * N_A # kJ ps / mol
const hbar = planck / 2pi
const kB = 1.380_649e-23 * to_kilo * N_A # kJ / mol K


text:
using OpenQuantumBase, Test
import LinearAlgebra: Diagonal

funcs = [(x) -> x, (x) -> 1 - x]
test_diag_operator = OpenQuantumBase.DiagonalOperator(funcs)
@test test_diag_operator(0.5) == OpenQuantumBase.Diagonal([0.5, 0.5])

test_geometric_operator = OpenQuantumBase.GeometricOperator(((x) -> -1.0im * x))
@test test_geometric_operator(0.5) == [0 0.5im; -0.5im 0]
@test_throws ArgumentError OpenQuantumBase.GeometricOperator((x) -> x, (x) -> 1 - x)

H = AdiabaticFrameHamiltonian(funcs, [])
@test H(2, 0.5) ≈ Diagonal([π, π])
H = AdiabaticFrameHamiltonian(funcs, nothing)
@test H(2, 0.0) ≈ Diagonal([0, 2π])
H = AdiabaticFrameHamiltonian((s)->[s, 1 - s], nothing)
@test H(2, 0.5) ≈ Diagonal([π, π])

# in_place update for matrices
du = [1.0 + 0.0im 0; 0 0]
u = PauliVec[2][1]
ρ = u * u'
hres = Diagonal([0, 2π])
H(du, ρ, 10, 0.0)
@test du ≈ -1.0im * (hres * ρ - ρ * hres)

dθ = (s) -> π / 2
gap = (s) -> (cos(2 * π * s) + 1) / 2
H = AdiabaticFrameHamiltonian([(x) -> -gap(x), (x) -> gap(x)], [dθ])
u = PauliVec[2][1]
ρ = u * u'

@test get_cache(H) == zeros(eltype(H), 2, 2)
@test size(H) == (2, 2)
@test H(10, 0.5) ≈ π * σx / 20
@test H(5, 0.0) ≈ π * σx / 10 - 2π * σz

# in_place update for vector
cache = get_cache(H)
update_cache!(cache, H, 10, 0.0)
@test cache ≈ -1.0im * (π * σx / 20 - 2π * σz)
update_cache!(cache, H, 10, 0.5)
@test cache ≈ -1.0im * (π * σx / 20)

# in_place update for matrices
du = [1.0 + 0.0im 0; 0 0]
hres = π * σx / 20 - 2π * σz
H(du, ρ, 10, 0.0)
@test du ≈ -1.0im * (hres * ρ - ρ * hres)


text:
################################################################################
# Script to generate julia code from the tox header files.                     #
# However these files need to be modified a little to get our final result.    #
################################################################################

using Clang.wrap_c

tox_headers = ["tox.h", "tox_old.h", "toxav.h", "toxdns.h"]

function header_wrapped(top_header, cursor_header)
      for tox_header in tox_headers
            if contains(cursor_header, tox_header)
                  return true
            end
      end

      return false
end

function cursor_wrapped(cursorname, cursor)
      return true
end

function header_outputfile(header_name)
      return "$(header_name).jl"
end

context = wrap_c.init(
            headers                 = ["cheader/tox.h", "cheader/toxav.h", "cheader/toxdns.h"],
            index                   = None,
            output_file             = "output.jl",
            common_file             = "common.jl",
            output_dir              = Pkg.dir("Toxcore", "generation", "cheader"),
            clang_args              = ["-pedantic", "-v"],
            clang_includes          = ["/usr/lib/gcc/x86_64-linux-gnu/4.8/include/"], # this is a hack to find the standard c libraries with the clang that is shipped with llvm and julia
            clang_diagnostics       = true,
            header_wrapped          = header_wrapped,
            header_library          = x->"libtoxcore",
            header_outputfile       = header_outputfile,
            cursor_wrapped          = cursor_wrapped,
            rewriter                = x -> x
        )

context.options.wrap_structs = false
run(context)


text:
# ....................... Header and Data input ..................

function _PRIMARY_input(dataobject::FITS_data)

    dataobject.hdutype == "PRIMARY" || error("FitsError: FITS_data is not of hdutype 'PRIMARY'")

    r::Array{String,1} = []

         E = Base.eltype(dataobject.data)
    nbytes = E.size
     nbits = 8 * nbytes
    bitpix = E <: AbstractFloat ? -abs(nbits) : nbits
    bitpix = Base.lpad(bitpix,20)
      dims = Base.size(dataobject.data)
     ndims = dims == (0,) ? 0 : Base.length(dims)
     naxis = Base.lpad(ndims,20)
      dims = ndims > 0 ? [Base.lpad(dims[i],20) for i=1:ndims] : 0
     bzero = Base.lpad(string(_fits_bzero(E)),20)

    incl = ndims > 0 ? true : false

           Base.push!(r,"SIMPLE  =                    T / file does conform to FITS standard             ")
    incl ? Base.push!(r,"BITPIX  = "   * bitpix  *    " / number of bits per data pixel                  ") : 0
           Base.push!(r,"NAXIS   = "   * naxis   *    " / number of data axes                            ")
    incl ? [Base.push!(r,"NAXIS$i  = " * dims[i] *    " / length of data axis " * rpad(i,27) ) for i=1:ndims] : 0
    incl ? Base.push!(r,"BZERO   = "   * bzero   *    " / offset data range to that of unsigned integer  ") : 0
    incl ? Base.push!(r,"BSCALE  =                  1.0 / default scaling factor                         ") : 0
           Base.push!(r,"EXTEND  =                    T / FITS dataset may contain extensions            ")
           Base.push!(r,"COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg               ")
           Base.push!(r,"END                                                                             ")

    return r          # r: Array{String,1} of records

end

function _IMAGE_input(data::Array{T,N} where {T <: Real,N})

    eltype(data) <: Real || error("FitsError: Array of real numbers expected")

         E = Base.eltype(data)
    nbytes = E.size
     nbits = 8 * nbytes
    bitpix = E <: AbstractFloat ? -abs(nbits) : nbits
    bitpix = Base.lpad(bitpix,20)
      dims = Base.size(data)
     ndims = dims == (0,) ? 0 : Base.length(dims)
      dims = ndims > 0 ? [Base.lpad(dims[i],20) for i=1:ndims] : 0
     naxis = Base.lpad(ndims,20)
     bzero = Base.lpad(string(_fits_bzero(E)),20)

    r::Array{String,1} = [];  incl = ndims > 0 ? true : false

           Base.push!(r,"XTENSION= 'IMAGE   '           / FITS standard extension                        ")
    incl ? Base.push!(r,"BITPIX  = "   * bitpix  *    " / number of bits per data pixel                  ") : 0
           Base.push!(r,"NAXIS   = "   * naxis   *    " / number of data axes                            ")
    incl ? [Base.push!(r,"NAXIS$i  = " * dims[i] *    " / length of data axis " * rpad(i,27) ) for i=1:ndims] : 0
    incl ? Base.push!(r,"PCOUNT  =                    0 / number of data axes                            ") : 0
    incl ? Base.push!(r,"GCOUNT  =                    1 / number of data axes                            ") : 0
    incl ? Base.push!(r,"BZERO   = "   * bzero   *    " / offset data range to that of unsigned integer  ") : 0
    incl ? Base.push!(r,"BSCALE  =                  1.0 / default scaling factor                         ") : 0
           Base.push!(r,"COMMENT    Extended FITS HDU   / http://fits.gsfc.nasa.gov/iaufwg               ")
           Base.push!(r,"END                                                                             ")

    return (r, data)           # r: Array{String,1} of records; data: same as input

end

function _table_data_types(cols::Vector{Vector{T} where T})

    ncols = length(cols)
    nrows = length(cols[1])
        f = Array{String,1}(undef,ncols)                                        # format specifier Xw.d

    for i ∈ eachindex(f)
        E = eltype(cols[i][1])
        x = E <: Integer ? "I" : E <: Real ? "E" : E == Float64 ? "D" : E <: Union{String,Char} ? "A" : "X"
        w = string(maximum([length(string(cols[i][j])) for j=1:nrows]))

        E <: Union{Char,String} ? (isascii(join(cols[i])) || error("FitsError: non-ASCII character in table $i")) : 0

        if E <: Union{Float16,Float32,Float64}
            v = string(cols[i][1])
            x = (('e' ∉ v) & ('p' ∉ v)) ? 'F' : x
            v = 'e' ∈ v ? split(v,'e')[1] : 'p' ∈ v ? split(v,'p')[1] : v
            d = !isnothing(findfirst('.',v)) ? string(length(split(v,'.')[2])) : '0'
        end

        f[i] = E <: Union{Float16,Float32,Float64} ? (x * w * '.' * d) : x * w
    end

    return f

end

function _TABLE_input(cols::Vector{Vector{T} where T})     # input array of table columns

    pcols = 1                                              # pointer to starting position of column in table row
    ncols = length(cols)                                   # number of columns
    nrows = length(cols[1])
    ncols < 1 && error("FitsError: a minimum of one column is mandatory")
    ncols = ncols < 999 ? ncols : 999
    ncols == 999 && println("FitsWarning: maximum number of columns exceeded (truncated at 999)")
    lcols = [length(cols[i]) for i=1:ncols]                          # length of columns (number of rows)
     pass = (sum(.!(lcols .== fill(nrows, ncols))) == 0)             # equal colum length test
     pass || error("FitsError: cannot create ASCII table (columns not of equal length)")

        w = [maximum([length(string(cols[i][j])) + 1  for j=1:nrows])  for i=1:ncols]
     data = [join([rpad(string(cols[i][j]),w[i])[1:w[i]] for i=1:ncols]) for j=1:nrows]
    tbcol = [pcols += w[i] for i=1:(ncols-1)]                        # field pointers (first column)
    tbcol = prepend!(tbcol,1)
    tbcol = [Base.lpad(tbcol[i],20) for i=1:ncols]
    tform = _table_data_types(cols)
    tform = ["'" * Base.rpad(tform[i],8) * "'"  for i=1:ncols]
    tform = [Base.rpad(tform[i],20) for i=1:ncols]
    ttype = ["HEAD$i"  for i=1:ncols]
    ttype = ["'" * Base.rpad(ttype[i],18) * "'" for i=1:ncols]          # default column headers
    nrows = Base.lpad(nrows,20)
    wcols = Base.lpad(sum(w),20)
    tcols = Base.lpad(ncols,20)

    r::Array{String,1} = []

    Base.push!(r,"XTENSION= 'TABLE   '           / FITS standard extension                        ")
    Base.push!(r,"BITPIX  =                    8 / number of bits per data pixel                  ")
    Base.push!(r,"NAXIS   =                    2 / number of data axes                            ")
    Base.push!(r,"NAXIS1  = "  * wcols    *    " / number of bytes/row                            ")
    Base.push!(r,"NAXIS2  = "  * nrows    *    " / number of rows                                 ")
    Base.push!(r,"PCOUNT  =                    0 / number of bytes in supplemetal data area       ")
    Base.push!(r,"GCOUNT  =                    1 / data blocks contain single table               ")
    Base.push!(r,"TFIELDS = "  * tcols    *    " / number of data fields (columns)                ")
    Base.push!(r,"COLSEP  =                    1 / number of spaces in column separator           ")
   [Base.push!(r,"TTYPE$i  = " * ttype[i] *    " / header of column " * rpad(i,30) ) for i=1:ncols]
   [Base.push!(r,"TBCOL$i  = " * tbcol[i] *    " / pointer to column " * rpad(i,29) ) for i=1:ncols]
   [Base.push!(r,"TFORM$i  = " * tform[i] *    " / data type of column " * rpad(i,27) ) for i=1:ncols]
   [Base.push!(r,"TDISP$i  = " * tform[i] *    " / data type of column " * rpad(i,27) ) for i=1:ncols]
    Base.push!(r,"COMMENT    Extended FITS HDU   / http://fits.gsfc.nasa.gov/iaufwg               ")
    Base.push!(r,"END                                                                             ")

    return (r, data)         # r: Array{String,1} of records; data: Array{String,1} of rows of table

end


text:
# ExoplanetsSysSim/src/summary_statistics.jl
## (c) 2015 Eric B. Ford

using Statistics

mutable struct CatalogSummaryStatistics
  stat::Dict{String,Any}          # For storing summary statistics
  cache::Dict{String,Any}         # For caching data that's not a summary statistic
end

function CatalogSummaryStatistics()    # Constructor for empty CatalogSummaryStatistics objects for creating globals to be used by closure
  CatalogSummaryStatistics( Dict{String,Any}(), Dict{String,Any}() )
end

function calc_summary_stats_sim_pass_one_demo(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, param::SimParam )      # Version for simulated data, since includes cat_phys
  ssd = Dict{String,Any}()
  cache = Dict{String,Any}()

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = findall(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = findall(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = findall(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  cache["num_tranets"] = num_tranets
  cache["idx_tranets"] = idx_tranets                                   # We can save lists of indices to summary stats for pass 2, even though we won't use these for computing a distance or probability
  #cache["idx_n_tranets"] = idx_n_tranets

  expected_num_detect = 0.0
  expected_num_sys_n_tranets = zeros(max_tranets_in_sys)
  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      p_tr_and_det = prob_detect(cat_obs.target[i].prob_detect,j)   # WARNING: Check why not using cat_phys here?
      expected_num_detect += p_tr_and_det
    end
    for k in 1:max_tranets_in_sys
      expected_num_sys_n_tranets[k] += prob_detect_n_planets(cat_obs.target[i].prob_detect,k)  # WARNING: Check why not use cat_phys here?
    end
  end
  ssd["expected planets detected"] = expected_num_detect
  ssd["num_sys_tranets"] = expected_num_sys_n_tranets
  ssd["num targets"] = get_int(param,"num_targets_sim_pass_one")
  #println("expected planets = ",expected_num_detect,", num_sys_tranets = ",expected_num_sys_n_tranets,", num targets = ",ssd["num targets"])

  # Arrays to store values for each tranet
  period_list = zeros(num_tranets)
  depth_list = zeros(num_tranets)
  weight_list = zeros(num_tranets)

  tr_id = 1   # tranet id
  for i in idx_tranets     # For each target with at least one tranet
     targ = cat_obs.target[i]
     for j in 1:min(length(targ.obs),max_tranets_in_sys)   # For each tranet around that target (but truncated if too many tranets in one system)
         #println("# i= ",i," j= ",j," tr_id= ",tr_id)
         period_list[tr_id] = targ.obs[j].period
         depth_list[tr_id] = targ.obs[j].depth
         # (s,p) = targ.phys_id[j]
         # ptr = calc_transit_prob_single(cat_phys.target[i],s,p)   # WARNING: Could access physical catalog, rather than observed catalog, but obviously that's dangerous for observations.
	       weight_list[tr_id] = prob_detect(cat_obs.target[i].prob_detect,j)
         tr_id += 1
      end
   end
  ssd["P list"] = period_list                                     # We can store whole lists, e.g., if we want to compute K-S distances
  ssd["depth list"] = depth_list
  ssd["weight list"] = weight_list

  idx_good = Bool[ period_list[i]>0.0 && depth_list[i]>0.0 && weight_list[i]>0.0 for i in 1:length(period_list) ]
  log_period_list = log10.(period_list[idx_good])
  log_depth_list = log10.(depth_list[idx_good])
  weight_list = weight_list[idx_good]
  weight_sum = sum(weight_list)
  ssd["mean log10 P"]  =  mean_log_P = sum( weight_list .* log_period_list) / weight_sum                           # TODO TEST: Check that these four weighted mean and stddevs are working properly
  ssd["std log10 P"]  =  sum( weight_list .* (log_period_list.-mean_log_P).^2 ) / weight_sum
  ssd["mean log10 depth"]  =  mean_log_depth = sum( weight_list .* log_depth_list) / weight_sum
  ssd["std log10 depth"]  =  sum( weight_list .* (log_depth_list.-mean_log_depth).^2 ) / weight_sum

  return CatalogSummaryStatistics(ssd, cache)
end

function calc_summary_stats_obs_demo(cat_obs::KeplerObsCatalog, param::SimParam )      # Version for observed data, thus no use of cat_phys
  ssd = Dict{String,Any}()
  cache = Dict{String,Any}()

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")                  # Demo that simulation parameters can specify how to evalute models, too
  idx_tranets = findall(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet
  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ [] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = findall(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = findall(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )
  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  if ( length( findall(x::KeplerTargetObs-> length(x.obs) > max_tranets_in_sys, cat_obs.target[idx_tranets] ) ) > 0)   # Make sure max_tranets_in_sys is at least big enough for observed systems
    warn("Observational data has more transiting planets in one systems than max_tranets_in_sys allows.")
  end
  num_tranets = Int64(num_tranets)           # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations
  #println("# num_tranets= ",num_tranets)

  # QUERY:  Is there any reason to cache anything for the real observations?  We only need to do this once, so might as well use one pass to simplicity.
  #cache["num_tranets"] = num_tranets
  #cache["idx_tranets"] = idx_tranets                                   # We can save lists of indices to summary stats for pass 2, even though we won't use these for computing a distance or probability
  cache["idx_n_tranets"] = idx_n_tranets

  ssd["planets detected"] = num_tranets                                 # WARNING: Note that we'll comparing two different things for simulated and real data during pass 1 (expected planets detected)
  num_sys_tranets = zeros(max_tranets_in_sys)                           # Since observed data, don't need to calculate probabilities.
  for n in 1:max_tranets_in_sys                                         # Make histogram of N-tranet systems
    num_sys_tranets[n] = length(idx_n_tranets[n])
  end
  ssd["num_sys_tranets"] = num_sys_tranets
  ssd["num targets"] = get_int(param,"num_kepler_targets")

  # Arrays to store values for each tranet
  period_list = zeros(num_tranets)
  depth_list = zeros(num_tranets)
  weight_list = ones(num_tranets)

   i = 1   # tranet id
   for targ in cat_obs.target[idx_tranets]                        # For each target with at least one tranet
     for j in 1:min(length(targ.obs),max_tranets_in_sys)          # For each tranet around that target (but truncated if too many tranets in one system)
         #println("# i= ",i," j= ",j)
         period_list[i] = targ.obs[j].period
         depth_list[i] = targ.obs[j].depth
         #weight_list[i] = 1.0
         i = i+1
      end
   end

  ssd["P list"] = period_list                                     # We can store whole lists, e.g., if we want to compute K-S distances
  ssd["depth list"] = depth_list
  ssd["weight list"] = weight_list

  idx_good = Bool[ period_list[i]>0.0 && depth_list[i]>0.0 for i in 1:length(period_list) ]
  log_period_list = log10.(period_list[idx_good])
  log_depth_list = log10.(depth_list[idx_good])
  ssd["mean log10 P"]  =  mean_log_P = mean(log_period_list)
  ssd["mean log10 depth"]  =  mean_log_depth = mean(log_depth_list)
  ssd["std log10 P"]  =  stdm(log_period_list,mean_log_P)
  ssd["std log10 depth"]  =  stdm(log_depth_list,mean_log_depth)

  return CatalogSummaryStatistics(ssd, cache)
end

# Just returns summary statistics passed, but provides a demo/hook for computing more expensive summary statistics if a model is good enouguh to be worth the extra time.
function calc_summary_stats_sim_pass_two_demo(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, ss::CatalogSummaryStatistics, param::SimParam )
  return ss
end

function test_summary_statistics(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, sim_param::SimParam)
  ss = calc_summary_stats_sim_pass_one_demo(cat_obs,cat_phys,sim_param)
  #println("len (ss pass 1)= ",length(collect(keys(ss1.stat))))
  ss = calc_summary_stats_sim_pass_two_demo(cat_obs,cat_phys,ss,sim_param)
  #println("len (ss pass 1)= ",length(collect(keys(ss1.stat))), "...   len (ss pass 2)= ",length(collect(keys(ss2.stat))) )
  ss = calc_summary_stats_obs_demo(cat_obs,sim_param)   # So tests can compare to simulated observed catalog
  return ss
end


text:
workspace()

using IntervalArithmetic
using EAGODomainReduction

Y = Interval{Float64}[Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0)]
ymult_lo = [50, 1.0, 2.0, 3.0]
ymult_hi = [0, 1.0, 2.0, 3.0]
yLBD = 1.0
yUBD = 3.0
Variable_DR!(Y,ymult_lo,ymult_hi,yLBD,yUBD)
println("Y post full duality-based bound tightening: ", Y)

X = Interval{Float64}[Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0)]
mult_lo = [50, 1.0]
mult_hi = [-1.0, 0]
LBD = 1.0
UBD = 2.0
nx = 2
Variable_DR_Imp!(X,mult_lo,mult_hi,LBD,UBD,nx)
println("X post implicit duality-based bound tightening: ", X)


text:
include("berryphase.jl")

using BerryPhase

const L = 4
const delta = 0.5
const rungs = 1.2:0.1:1.7

const M = (length(ARGS)>0 ? parse(Int, ARGS[1]) : 10)

bp_leg(rung) = berryphase(altladder(1, L, 1.0, 1.0, delta, rung, rung, :plus_leg), M)
bp_rung(rung) = berryphase(altladder(1, L, 1.0, 1.0, delta, rung, rung, :rung), M)

open("res.dat", "w") do io
  for r in rungs
    bl = bp_leg(r)
    br = bp_rung(r)
    println(r, " ", bl/pi, " ", br/pi)
    println(io, r, " ", bl/pi, " ", br/pi)
  end
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
struct TermCombination
    terms::Set{Dict{Int, Int}} # idx => pow
end

@eval Base.zero(::Type{TermCombination}) = $(TermCombination(Set{Dict{Int,Int}}()))
@eval Base.one(::Type{TermCombination}) = $(TermCombination(Set([Dict{Int,Int}()])))

function Base.:(==)(comb1::TermCombination, comb2::TermCombination)
    comb1.terms == comb2.terms && return true

    n1 = reduce(max, (k for (k,_) in Iterators.flatten(comb1.terms)), init=0)
    n2 = reduce(max, (k for (k,_) in Iterators.flatten(comb2.terms)), init=0)
    n = max(n1, n2)

    _sparse(comb1, n) == _sparse(comb2, n)
end

function Base.:+(comb1::TermCombination, comb2::TermCombination)
    if isone(comb1) && !iszero(comb2)
        return comb2
    elseif isone(comb2) && !iszero(comb1)
        return comb1
    elseif comb1 === comb2
        return comb1
    end
    TermCombination(union(comb1.terms, comb2.terms))
end

Base.:+(comb1::TermCombination) = comb1

function _merge(dict1, dict2)
    d = copy(dict1)
    for (k, v) in dict2
        d[k] = min(2, get(dict1, k, 0) + v)
    end
    d
end

function Base.:*(comb1::TermCombination, comb2::TermCombination)
    if isone(comb1)
        return comb2
    elseif isone(comb2)
        return comb1
    elseif comb1 === comb2 # squaring optimization
        terms = comb1.terms
        # turns out it's enough to track
        # a^2*b^2
        # and a^2 + b^2 + ab
        # have the same hessian sparsity
        t = Dict(k=>2 for (k,_) in
                 Iterators.flatten(terms))
        TermCombination(Set([t]))
        #=
        # square each term
        t1 = [Dict(k=>2 for (k,_) in dict)
              for dict in comb1.terms]
        # multiply each term
        t2 = Dict{Int,Int}[]
        for i in 1:length(terms)
            for j in i+1:length(terms)
                push!(t2, _merge(terms[i], terms[j]))
            end
        end
        TermCombination(union(t1, t2))
        =#
    else
        Set([_merge(dict1, dict2)
             for dict1 in comb1.terms,
             dict2 in comb2.terms]) |> TermCombination
    end
end
Base.:*(comb1::TermCombination) = comb1
Base.iszero(c::TermCombination) = isempty(c.terms)
Base.isone(c::TermCombination) = all(isempty, c.terms)

function _sparse(t::TermCombination, n)
    I = Int[]
    J = Int[]
    for dict in t.terms
        kv = collect(pairs(dict))
        for i in 1:length(kv)
            k, v = kv[i]
            if v>=2
                push!(I, k)
                push!(J, k)
            end
            for j in i+1:length(kv)
                if v >= 1 && kv[j][2] >= 1
                    push!(I, k)
                    push!(J, kv[j][1])
                end
            end
        end
    end
    s1 = sparse(I,J,fill!(BitVector(undef, length(I)), true),n,n)
    s1 .| s1'
end


text:
# Test the differentiate function
using Test
using SAC

@testset "Differentiation" begin
    let sample_data = SAC.sample()

        # Don't accept wrong arguments
        for i in (0, 1, 4, 6)
            @test_throws ArgumentError SAC.differentiate(sample_data, i)
        end

        # Compare against output from SAC
        for (i, str) in zip((2, 3, 5), ("two", "three", "five"))
            local file = joinpath(dirname(@__FILE__), "data", "test_diff_data_points_$(str).sac")
            isfile(file) || error("Test data file '$(file)' for differentiate does not exist")
            local reference_data = SAC.read(file)
            @test SAC.differentiate(sample_data, i).t ≈ reference_data.t
        end
    end
end


text:
using ITensors,
  ITensorGPU,
  LinearAlgebra, # For tr()
  Combinatorics, # For permutations()
  Random,
  CUDA,
  Test

# gpu tests!
@testset "cuITensor, Dense{$SType} storage" for SType in (Float64, ComplexF64)
  mi, mj, mk, ml, ma = 2, 3, 4, 5, 6, 7
  i = Index(mi, "i")
  j = Index(mj, "j")
  k = Index(mk, "k")
  l = Index(ml, "l")
  a = Index(ma, "a")
  @testset "Constructor" begin
    A = cuITensor(one(SType), i, j, k)
    @test collect(CuArray(A, i, j, k)) == ones(SType, dim(i), dim(j), dim(k))
    A = randomCuITensor(IndexSet(i, j, k))
    @test inds(A) == IndexSet(i, j, k)
    @test ITensorGPU.storage(A) isa ITensorGPU.CuDense
    Aarr = rand(SType, dim(i) * dim(j) * dim(k))
    @test cpu(ITensor(Aarr, i, j, k)) == cpu(cuITensor(Aarr, i, j, k))
    @test cuITensor(SType, i, j, k) isa ITensor
    @test storage(cuITensor(SType, i, j, k)) isa ITensorGPU.CuDense{SType}
    @test vec(collect(CuArray(ITensor(Aarr, i, j, k), i, j, k))) == Aarr
  end
  @testset "Test permute(cuITensor,Index...)" begin
    CA = randomCuITensor(SType, i, k, j)
    permCA = permute(CA, k, j, i)
    permA = cpu(permCA)
    @test k == inds(permA)[1]
    @test j == inds(permA)[2]
    @test i == inds(permA)[3]
    A = cpu(CA)
    for ii in 1:dim(i), jj in 1:dim(j), kk in 1:dim(k)
      @test A[k => kk, i => ii, j => jj] == permA[i => ii, j => jj, k => kk]
    end
  end
  @testset "Test permute(cuITensor,Index...) for large tensors" begin
    inds = [Index(2) for ii in 1:14]
    A = randomITensor(SType, IndexSet(inds))
    CA = cuITensor(A)
    for shuffle_count in 1:20
      perm_inds = shuffle(inds)
      permCA = permute(CA, perm_inds...)
      permA = cpu(permCA)
      pA = permute(A, perm_inds...)
      for ci in CartesianIndices(pA)
        @test pA[ci] == permA[ci]
      end
    end
  end
  #=@testset "Test scalar(cuITensor)" begin
    x = SType(34)
    A = randomCuITensor(a)
    @test x==scalar(A)
  end=#
  @testset "Test CuVector(cuITensor)" begin
    v = CuVector(ones(SType, dim(a)))
    A = cuITensor(v, a)
    @test v == CuVector(A)
  end
  @testset "Test CuMatrix(cuITensor)" begin
    v = CuMatrix(ones(SType, dim(a), dim(l)))
    A = cuITensor(vec(v), a, l)
    @test v == CuMatrix(A, a, l)
    A = cuITensor(vec(v), a, l)
    @test v == CuMatrix(A)
    A = cuITensor(vec(v), a, l)
    @test v == CuArray(A, a, l)
    @test v == CuArray(A)
  end
  @testset "Test norm(cuITensor)" begin
    A = randomCuITensor(SType, i, j, k)
    B = dag(A) * A
    @test norm(A) ≈ sqrt(real(scalar(B)))
  end
  @testset "Test complex(cuITensor)" begin
    A = randomCuITensor(SType, i, j, k)
    cA = complex(A)
    @test complex.(CuArray(A)) == CuArray(cA)
  end
  #@testset "Test exp(cuITensor)" begin
  #  A  = randomCuITensor(SType,i,i')
  #  @test CuArray(exp(A,i,i')) ≈ exp(CuArray(A))
  #end
  @testset "Test add cuITensors" begin
    dA = randomCuITensor(SType, i, j, k)
    dB = randomCuITensor(SType, k, i, j)
    A = cpu(dA)
    B = cpu(dB)
    C = cpu(dA + dB)
    @test CuArray(permute(C, i, j, k)) ==
      CuArray(permute(A, i, j, k)) + CuArray(permute(B, i, j, k))
    for ii in 1:dim(i), jj in 1:dim(j), kk in 1:dim(k)
      @test C[i => ii, j => jj, k => kk] ==
        A[j => jj, i => ii, k => kk] + B[i => ii, k => kk, j => jj]
    end
  end

  @testset "Test factorizations of a cuITensor" begin
    A = randomCuITensor(SType, i, j, k, l)
    @testset "Test SVD of a cuITensor" begin
      U, S, V = svd(A, (j, l))
      u = commonind(U, S)
      v = commonind(S, V)
      @test cpu(A) ≈ cpu(U * S * V)
      @test cpu(U * dag(prime(U, u))) ≈ δ(SType, u, u') rtol = 1e-14
      @test cpu(V * dag(prime(V, v))) ≈ δ(SType, v, v') rtol = 1e-14
    end

    A = randomCuITensor(SType, i, j, k, l)
    @testset "Test SVD consistency between CPU and GPU" begin
      U_gpu, S_gpu, V_gpu = svd(A, (j, l))
      U_cpu, S_cpu, V_cpu = svd(cpu(A), (j, l))
      @test cpu(U_gpu) * cpu(S_gpu) * cpu(V_gpu) ≈ U_cpu * S_cpu * V_cpu
    end

    #=@testset "Test SVD truncation" begin
        M = randn(4,4)
        (U,s,V) = svd(M)
        ii = Index(4)
        jj = Index(4)
        S = Diagonal(s)
        T = cuITensor(vec(CuArray(U*S*V')),IndexSet(ii,jj))
        (U,S,V) = svd(T,ii;maxm=2)
        @test norm(U*S*V-T)≈sqrt(s[3]^2+s[4]^2)
    end=#

    @testset "Test QR decomposition of a cuITensor" begin
      Q, R = qr(A, (i, l))
      q = commonind(Q, R)
      @test cpu(A) ≈ cpu(Q * R)
      @test cpu(Q * dag(prime(Q, q))) ≈ δ(SType, q, q') atol = 1e-14
    end

    #=@testset "Test polar decomposition of a cuITensor" begin
      U,P = polar(A,(k,l))
      @test cpu(A)≈cpu(U*P)
    end=#

  end # End ITensor factorization testset
end # End Dense storage test