text:
import MathOptInterface
import SparseArrays

const MOI = MathOptInterface

MOI.Utilities.@product_of_sets(
    _LPProductOfSets,
    MOI.EqualTo{T},
    MOI.LessThan{T},
    MOI.GreaterThan{T},
    MOI.Interval{T},
)

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Cint,
            MOI.Utilities.ZeroBasedIndexing,
        },
        MOI.Utilities.Hyperrectangle{Float64},
        _LPProductOfSets{Float64},
    },
}

const SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Ptr{Cvoid}
    solver_options::Ptr{Cvoid}
    options_set::Set{String}
    optimize_called::Bool
    solve_time::Float64
    # Work-around for upstream bug in Clp:
    maximumSeconds::Float64

    """
        Optimizer()

    Create a new Optimizer object.

    Set optimizer attributes using `MOI.RawOptimizerAttribute` or
    `JuMP.set_optimizer_atttribute`. For a list of supported parameter names,
    see `Clp.SUPPORTED_PARAMETERS`.

    ## Example

    ```julia
    using JuMP, Clp
    model = JuMP.Model(Clp.Optimizer)
    set_optimizer_attribute(model, "LogLevel", 0)
    ```
    """
    function Optimizer(; kwargs...)
        model =
            new(Clp_newModel(), ClpSolve_new(), Set{String}(), false, 0.0, -1.0)
        if length(kwargs) > 0
            @warn("""Passing optimizer attributes as keyword arguments to
            Clp.Optimizer is deprecated. Use
                MOI.set(model, MOI.RawOptimizerAttribute("key"), value)
            or
                JuMP.set_optimizer_attribute(model, "key", value)
            instead.
            """)
        end
        for (key, value) in kwargs
            MOI.set(model, MOI.RawOptimizerAttribute(String(key)), value)
        end
        finalizer(model) do m
            Clp_deleteModel(m)
            ClpSolve_delete(m.solver_options)
            return
        end
        return model
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, model::Optimizer) = model
function Base.unsafe_convert(::Type{Ptr{Cvoid}}, model::Optimizer)
    return model.inner::Ptr{Cvoid}
end

# ====================
#   empty functions
# ====================

function MOI.is_empty(model::Optimizer)
    # A problem is empty if it has no variable and no linear constraints
    return (Clp_getNumRows(model) == 0) && (Clp_getNumCols(model) == 0)
end

function MOI.empty!(model::Optimizer)
    # Copy parameters from old model into new model
    old_options = Dict(
        key => MOI.get(model, MOI.RawOptimizerAttribute(key)) for
        key in model.options_set
    )
    empty!(model.options_set)
    Clp_deleteModel(model)
    model.inner = Clp_newModel()
    model.optimize_called = false
    model.solve_time = 0.0
    for (key, value) in old_options
        MOI.set(model, MOI.RawOptimizerAttribute(key), value)
    end
    # Work-around for maximumSeconds
    Clp_setMaximumSeconds(model, model.maximumSeconds)
    return
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Clp"

# If you update this list, remember to update the README documentation!
const SUPPORTED_PARAMETERS = (
    "PrimalTolerance",
    "DualTolerance",
    "DualObjectiveLimit",
    "MaximumIterations",
    "MaximumSeconds",
    "LogLevel",
    "Scaling",
    "Perturbation",
    # TODO(odow): `Algorithm` is excluded from the README because it isn't
    # apparent what it controls. Use `SolveType` instead.
    "Algorithm",
    "PresolveType",
    "SolveType",
    "InfeasibleReturn",
)

function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute)
    return param.name in SUPPORTED_PARAMETERS
end

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    name = String(param.name)
    push!(model.options_set, name)
    if name == "PrimalTolerance"
        Clp_setPrimalTolerance(model, value)
    elseif name == "DualTolerance"
        Clp_setDualTolerance(model, value)
    elseif name == "DualObjectiveLimit"
        Clp_setDualObjectiveLimit(model, value)
    elseif name == "MaximumIterations"
        Clp_setMaximumIterations(model, value)
    elseif name == "MaximumSeconds"
        Clp_setMaximumSeconds(model, value)
    elseif name == "LogLevel"
        Clp_setLogLevel(model, value)
    elseif name == "Scaling"
        Clp_scaling(model, value)
    elseif name == "Perturbation"
        Clp_setPerturbation(model, value)
    elseif name == "Algorithm"
        Clp_setAlgorithm(model, value)
    elseif name == "PresolveType"
        ClpSolve_setPresolveType(model.solver_options, value, -1)
    elseif name == "SolveType"
        ClpSolve_setSolveType(model.solver_options, value, -1)
    elseif name == "InfeasibleReturn"
        ClpSolve_setInfeasibleReturn(model.solver_options, value)
    else
        throw(MOI.UnsupportedAttribute(param))
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    name = String(param.name)
    if name == "PrimalTolerance"
        return Clp_primalTolerance(model)
    elseif name == "DualTolerance"
        return Clp_dualTolerance(model)
    elseif name == "DualObjectiveLimit"
        return Clp_dualObjectiveLimit(model)
    elseif name == "MaximumIterations"
        # TODO(odow): Clp doesn't follow convention with maximumIterations!
        return maximumIterations(model)
    elseif name == "MaximumSeconds"
        return Clp_maximumSeconds(model)
    elseif name == "LogLevel"
        return Clp_logLevel(model)
    elseif name == "Scaling"
        return Clp_scalingFlag(model)
    elseif name == "Perturbation"
        return Clp_perturbation(model)
    elseif name == "Algorithm"
        return Clp_algorithm(model)
    elseif name == "PresolveType"
        return ClpSolve_getPresolveType(model.solver_options)
    elseif name == "SolveType"
        return ClpSolve_getSolveType(model.solver_options)
    elseif name == "InfeasibleReturn"
        return ClpSolve_infeasibleReturn(model.solver_options)
    end
    return throw(MOI.UnsupportedAttribute(param))
end

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    push!(model.options_set, "LogLevel")
    Clp_setLogLevel(model, value ? 0 : 1)
    return
end

MOI.get(model::Optimizer, ::MOI.Silent) = Clp_logLevel(model) == 0

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    push!(model.options_set, "MaximumSeconds")
    value = value === nothing ? -1.0 : value
    Clp_setMaximumSeconds(model, value)
    model.maximumSeconds = value
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    # TODO(odow): replace with `Clp_maximumSeconds(model)` when upstream
    # is fixed.
    return model.maximumSeconds
end

MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = false

# ========================================
#   Supported constraints and objectives
# ========================================

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{<:Union{MOI.VariableIndex,MOI.ScalarAffineFunction{Float64}}},
    ::Type{<:SCALAR_SETS},
)
    return true
end

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

# =======================
#   `copy_to` function
# =======================

function _index_map(
    src::OptimizerCache,
    index_map,
    ::Type{F},
    ::Type{S},
) where {F,S}
    for ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        row = MOI.Utilities.rows(src.constraints, ci)
        index_map[ci] = MOI.ConstraintIndex{F,S}(row)
    end
    return
end

function _index_map(
    src::OptimizerCache,
    index_map,
    F::Type{MOI.VariableIndex},
    ::Type{S},
) where {S}
    for ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        col = index_map[MOI.VariableIndex(ci.value)].value
        index_map[ci] = MOI.ConstraintIndex{F,S}(col)
    end
    return
end

"""
    _index_map(src::OptimizerCache)

Create an `IndexMap` mapping the variables and constraints in `OptimizerCache`
to their corresponding 1-based columns and rows.
"""
function _index_map(src::OptimizerCache)
    index_map = MOI.IndexMap()
    for (i, x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
        index_map[x] = MOI.VariableIndex(i)
    end
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        _index_map(src, index_map, F, S)
    end
    return index_map
end

function MOI.copy_to(dest::Optimizer, src::OptimizerCache)
    @assert MOI.is_empty(dest)
    A = src.constraints.coefficients
    row_bounds = src.constraints.constants
    obj =
        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    c = zeros(A.n)
    for term in obj.terms
        c[term.variable.value] += term.coefficient
    end
    Clp_setObjectiveOffset(dest, -obj.constant)
    Clp_loadProblem(
        dest,
        A.n,
        A.m,
        A.colptr,
        A.rowval,
        A.nzval,
        src.variables.lower,
        src.variables.upper,
        c,
        row_bounds.lower,
        row_bounds.upper,
    )
    sense = MOI.get(src, MOI.ObjectiveSense())
    if sense == MOI.MIN_SENSE
        Clp_setObjSense(dest, 1)
    elseif sense == MOI.MAX_SENSE
        Clp_setObjSense(dest, -1)
    else
        @assert sense == MOI.FEASIBILITY_SENSE
        Clp_setObjSense(dest, 0)
    end
    return _index_map(src)
end

function MOI.copy_to(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache},
)
    return MOI.copy_to(dest, src.model)
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    src_cache = MOI.copy_to(cache, src)
    cache_dest = MOI.copy_to(dest, cache)
    index_map = MOI.IndexMap()
    for (src_x, cache_x) in src_cache.var_map
        index_map[src_x] = cache_dest[cache_x]
    end
    for (src_ci, cache_ci) in src_cache.con_map
        index_map[src_ci] = cache_dest[cache_ci]
    end
    return index_map
end

# ===============================
#   Optimize and post-optimize
# ===============================

function MOI.optimize!(model::Optimizer)
    t = time()
    Clp_initialSolveWithOptions(model, model.solver_options)
    model.solve_time = time() - t
    model.optimize_called = true
    return
end

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = Clp_getNumCols(model)

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return Clp_getObjValue(model)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if !model.optimize_called
        return MOI.OPTIMIZE_NOT_CALLED
    end
    st = Clp_status(model)
    if st == 0
        return MOI.OPTIMAL
    elseif st == 1
        return MOI.INFEASIBLE
    elseif st == 2
        return MOI.DUAL_INFEASIBLE
    elseif st == 3
        # No more granular information that "some limit is reached"
        return MOI.OTHER_LIMIT
    end
    return MOI.OTHER_ERROR
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    if !model.optimize_called
        return "MOI.OPTIMIZE_NOT_CALLED"
    end
    st = Clp_status(model)
    if st == 0
        return "0 - optimal"
    elseif st == 1
        return "1 - primal infeasible"
    elseif st == 2
        return "2 - dual infeasible"
    elseif st == 3
        return "3 - stopped on iterations etc"
    else
        @assert st == 4
        return "4 - stopped due to errors"
    end
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    if Clp_primalFeasible(model) != 0
        return 1
    elseif Clp_dualFeasible(model) != 0
        return 1
    elseif Clp_isProvenPrimalInfeasible(model) != 0
        return 1
    elseif Clp_isProvenDualInfeasible(model) != 0
        return 1
    end
    return 0
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif Clp_isProvenDualInfeasible(model) != 0
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp_primalFeasible(model) != 0
        return MOI.FEASIBLE_POINT
    end
    return MOI.UNKNOWN_RESULT_STATUS
end

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif Clp_isProvenPrimalInfeasible(model) != 0
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp_dualFeasible(model) != 0
        return MOI.FEASIBLE_POINT
    end
    return MOI.UNKNOWN_RESULT_STATUS
end

# ===================
#   Primal solution
# ===================

function _unsafe_wrap_clp_array(
    model::Optimizer,
    f::Function,
    n::Integer,
    indices = nothing;
    own::Bool = false,
)
    p = f(model)
    if p == C_NULL
        return map(x -> NaN, indices)
    end
    x = unsafe_wrap(Array, p, (n,); own = own)
    return indices === nothing ? x : x[indices]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_unboundedRay.
        return _unsafe_wrap_clp_array(
            model,
            Clp_unboundedRay,
            Clp_getNumCols(model),
            x.value;
            own = true,
        )
    elseif primal_status == MOI.FEASIBLE_POINT
        return _unsafe_wrap_clp_array(
            model,
            Clp_getColSolution,
            Clp_getNumCols(model),
            x.value,
        )
    end
    return error("Primal solution not available")
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    xs::Vector{MOI.VariableIndex},
)
    MOI.check_result_index_bounds(model, attr)
    col_indices = [idx.value for idx in xs]
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_unboundedRay.
        return _unsafe_wrap_clp_array(
            model,
            Clp_unboundedRay,
            Clp_getNumCols(model),
            col_indices;
            own = true,
        )
    elseif primal_status == MOI.FEASIBLE_POINT
        return _unsafe_wrap_clp_array(
            model,
            Clp_getColSolution,
            Clp_getNumCols(model),
            col_indices,
        )
    end
    return error("Primal solution not available")
end

# TODO: What happens if model is unbounded / infeasible?
function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:SCALAR_SETS},
)
    MOI.check_result_index_bounds(model, attr)
    return _unsafe_wrap_clp_array(
        model,
        Clp_getRowActivity,
        Clp_getNumRows(model),
        c.value,
    )
end

# TODO: What happens if model is unbounded / infeasible?
function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.VariableIndex,<:SCALAR_SETS},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end

# =================
#   Dual solution
# =================

# If sense is maximize, we negate all the duals to follow MOI conventions
# Feasibility problems are treated as a minimization

"""
    _farkas_variable_dual(model::Optimizer, col::Cint)

Return a Farkas dual associated with the variable bounds of `col`.

Compute the Farkas dual as:

    ā * x = λ' * A * x <= λ' * b = -β + sum(āᵢ * Uᵢ | āᵢ < 0) + sum(āᵢ * Lᵢ | āᵢ > 0)

The Farkas dual of the variable is ā, and it applies to the upper bound if ā < 0,
and it applies to the lower bound if ā > 0.
"""
function _farkas_variable_dual(model::Optimizer, col::Integer)
    m, n = Clp_getNumRows(model), Clp_getNumCols(model)
    nnz = Clp_getNumElements(model)
    vbeg = _unsafe_wrap_clp_array(model, Clp_getVectorStarts, n)
    vlen = _unsafe_wrap_clp_array(model, Clp_getVectorLengths, n)
    indices = vbeg[col] .+ (1:vlen[col])
    vind = _unsafe_wrap_clp_array(model, Clp_getIndices, nnz, indices)
    vval = _unsafe_wrap_clp_array(model, Clp_getElements, nnz, indices)
    # We need to claim ownership of the pointer returned by Clp_infeasibilityRay.
    λ = _unsafe_wrap_clp_array(model, Clp_infeasibilityRay, m; own = true)
    return sum(v * λ[i+1] for (i, v) in zip(vind, vval))
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:SCALAR_SETS},
)
    MOI.check_result_index_bounds(model, attr)
    n = Clp_getNumRows(model)
    dual_status = MOI.get(model, MOI.DualStatus())
    sense = Clp_getObjSense(model)
    if dual_status == MOI.FEASIBLE_POINT
        dsol = _unsafe_wrap_clp_array(model, Clp_getRowPrice, n, c.value)
        return sense * dsol
    elseif dual_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_infeasibilityRay.
        return -_unsafe_wrap_clp_array(
            model,
            Clp_infeasibilityRay,
            n,
            c.value;
            own = true,
        )
    end
    return error("Dual solution not available")
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        return min(0.0, _farkas_variable_dual(model, c.value))
    end
    rc = _unsafe_wrap_clp_array(
        model,
        Clp_getReducedCost,
        Clp_getNumCols(model),
        c.value,
    )
    return min(0.0, Clp_getObjSense(model) * rc)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        return max(0.0, _farkas_variable_dual(model, c.value))
    end
    rc = _unsafe_wrap_clp_array(
        model,
        Clp_getReducedCost,
        Clp_getNumCols(model),
        c.value,
    )
    return max(0.0, Clp_getObjSense(model) * rc)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{
        MOI.VariableIndex,
        <:Union{MOI.Interval{Float64},MOI.EqualTo{Float64}},
    },
)
    MOI.check_result_index_bounds(model, attr)
    if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        return _farkas_variable_dual(model, c.value)
    end
    rc = _unsafe_wrap_clp_array(
        model,
        Clp_getReducedCost,
        Clp_getNumCols(model),
        c.value,
    )
    return Clp_getObjSense(model) * rc
end

# Corresponds to the `Status` enum defined in
# https://github.com/coin-or/Clp/blob/8419e63/Clp/src/ClpSimplex.hpp#L114
const _CLP_BASIS_STATUS = Dict(
    # isFree
    0x00 => MOI.BASIC,
    # basic
    0x01 => MOI.BASIC,
    # atUpperBound
    0x02 => MOI.NONBASIC_AT_UPPER,
    # atLowerBound
    0x03 => MOI.NONBASIC_AT_LOWER,
    # superBasic
    0x04 => MOI.SUPER_BASIC,
    # isFixed
    0x05 => MOI.NONBASIC,
)

function _nonbasic_status(status, ::Type{<:MOI.LessThan})
    return status == MOI.NONBASIC_AT_LOWER ? MOI.BASIC : MOI.NONBASIC
end

function _nonbasic_status(status, ::Type{<:MOI.GreaterThan})
    return status == MOI.NONBASIC_AT_LOWER ? MOI.NONBASIC : MOI.BASIC
end

_nonbasic_status(::Any, ::Type{<:MOI.EqualTo}) = MOI.NONBASIC

_nonbasic_status(status, ::Type{<:MOI.Interval}) = status

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintBasisStatus,
    c::MOI.ConstraintIndex{F,S},
) where {F,S}
    code = Clp_getRowStatus(model, c.value - 1)
    status = _CLP_BASIS_STATUS[code]
    if status == MOI.NONBASIC_AT_UPPER || status == MOI.NONBASIC_AT_LOWER
        return _nonbasic_status(status, S)
    end
    return status
end

function MOI.get(
    model::Optimizer,
    ::MOI.VariableBasisStatus,
    vi::MOI.VariableIndex,
)
    code = Clp_getColumnStatus(model, vi.value - 1)
    return _CLP_BASIS_STATUS[code]
end


text:
"""
1D Current Source Density for a set of voltage traces from a linear array of electrodes.

1. data: ch x sample or ch x sample x epoch(volts)
- method:

CSD: finite difference method of second spatial derivatives of potentials, I = - ∇σ⋅∇ϕ,
      σ is conductivities, ϕ is potentials, on the continuous assumption when spacing is around 20-100 um.
     (Nicholson & Freeman, 1975, J Neurophysiol, 38(2): 356-68)

iCSD: solve I -> ϕ forward problem, then get the inverse transformation matrix.
      (Petterson et al., 2006, J Neurosci Methods, 154(1-2):116-33)

      iCSDdelta: δ method where infinitely thin current source discs of radius r.
      iCSDspline: cubic spline method where current source between electrode be modeled by smooth spline

kCSD: kernal CSD
      ()

- h: channel spacing(micrometers)
- c: conductivity of the extracellular medium in siemans/meter
- r: radius(micrometers) of the model current source/sink in reverse method

- return: CSD(amps/meters^3)
"""
function csd(data;method=:iCSDdelta,h=20,c=0.3,r=500)
  nd=ndims(data)
  if nd==3
    nr,nc,n = size(data)
    y = Array{Float64}(undef,nr,nc,n)
    for i in 1:n
      y[:,:,i] = csd(data[:,:,i],method=method,h=h,c=c,r=r)
    end
  else
    nr,nc = size(data)
    𝒉 = h*1e-6
    𝒓 = r*1e-6
    if method==:iCSDdelta
      f = [𝒉^2/(2c) * ( sqrt((j-i)^2 + (𝒓/𝒉)^2) - abs(j-i) ) for j=1:nr,i=1:nr]
      y = f\data
    elseif method==:iCSDspline # have wired problems, doesn't replicate the original code(CSDPlotter) results
      ep = collect(1:nr)*h
      y,zs = csd_cubicspline(ep,data,f_cubicspline(ep,r,c)...)
    elseif method==:kCSD

    else
      y=zeros(nr,nc)
      for i in 2:nr-1
        y[i,:] = -c*(data[i-1,:] + data[i+1,:] - 2*data[i,:]) / 𝒉^2
      end
    end
  end
  return y
end

"""
cubic spline inverse CSD

contact_positions: contact positions
data: measured potentials
Fcs: cubic spline transformation matrix corresponding to the given contact positions
E0,E1,E2,E3: matrixes containing the "recursive rules"
"""
function csd_cubicspline(contact_positions,data,Fcs,E0,E1,E2,E3,num_out_zs=size(data,1))
  N,num_of_timesteps = size(data)
  cs_data = zeros(N+2,num_of_timesteps)
  cs_data[2:N+1,:] = data # ϕ_1 = ϕ_N+2 = 0
  h = mean(diff(contact_positions))
  el_pos_with_ends = vcat(0,contact_positions,contact_positions[end]+h)

  CSD_coeff = Fcs\cs_data
  # The cubic spline polynomial coeffescients
  A0 = E0*CSD_coeff
  A1 = E1*CSD_coeff
  A2 = E2*CSD_coeff
  A3 = E3*CSD_coeff

  out_zs = el_pos_with_ends[1]:(el_pos_with_ends[N+2]-el_pos_with_ends[1])/(num_out_zs-1):el_pos_with_ends[N+2]
  CSD = zeros(length(out_zs),num_of_timesteps)

  i = 1
  for j=1:length(out_zs)
    if out_zs[j]>el_pos_with_ends[i+1]
      i+=1
    end
    CSD[j,:] = A0[i,:] + A1[i,:]*(out_zs[j]-el_pos_with_ends[i]) +
    A2[i,:]*(out_zs[j]-el_pos_with_ends[i])^2 + A3[i,:]*(out_zs[j]-el_pos_with_ends[i])^3
  end
  return CSD,out_zs
end

"""
Computes the E0, E1, E2 and E3 matrixes used in the cubic spline iCSD method.
These matrixes contains the recursive formulas for finding the F matrix
"""
function e_matrixes(el_pos)
  N = length(el_pos)
  h = mean(diff(el_pos))
  z_js = vcat(0,el_pos,el_pos[end]+h) # electrode positions with two imaginary

  C_vec = 1/diff(z_js) # length: N+1
  # Define transformation matrixes
  C_jm1 = zeros(N+2,N+2)
  C_j0 = zeros(N+2,N+2)
  C_mat3 = zeros(N+1,N+1)

  for i=1:N+1
    for j=1:N+1
      if i == j
        C_jm1[i+1,j+1] = C_vec[i]
        C_j0[i,j] = C_jm1[i+1,j+1]
        C_mat3[i,j] = C_vec[i]
      end
    end
  end
  C_jm1[N+2,N+2] = 0
  C_j0[1,1] = 0

  C_jall = copy(C_j0)
  C_jall[1,1] = 1
  C_jall[N+2,N+2] = 1

  Tjp1 = zeros(N+2,N+2) # converting an element k_j to k_j+1
  Tjm1 = zeros(N+2,N+2) # converting an element k_j to k_j-1
  Tj0  = Matrix{Float64}(I,N+2,N+2)
  Tj0[1,1] = 0
  Tj0[N+2,N+2] = 0

  # C to K
  for i=2:N+2
    for j=1:N+2
      if i==j-1
        Tjp1[i,j] = 1
      end
      if i==j+1
        Tjm1[i,j] = 1
      end
    end
  end

  # C to K transformation matrix
  K = ( C_jm1*Tjm1 + 2*C_jm1*Tj0 + 2*C_jall + C_j0*Tjp1 ) \ (3*( C_jm1^2*Tj0 - C_jm1^2*Tjm1 + C_j0^2*Tjp1 - C_j0^2*Tj0 ))

  # Define matrixes for C to A transformation
  Tja  = zeros(N+1,N+2)      # identity matrix except that it cuts off last element
  Tjp1a  = zeros(N+1,N+2)    # converting k_j to k_j+1 and cutting off last element

  # C to A
  for i=1:N+1
    for j=1:N+2
      if i==j-1
        Tjp1a[i,j] = 1
      end
      if i==j
        Tja[i,j] = 1
      end
    end
  end

  # Define spline coeffiscients
  E0  = Tja
  E1  = Tja*K
  E2  = 3*C_mat3^2*(Tjp1a-Tja) - C_mat3*(Tjp1a+2*Tja)*K
  E3  = 2*C_mat3^3*(Tja-Tjp1a) + C_mat3^2*(Tjp1a+Tja)*K
  return E0,E1,E2,E3
end

"F and E matrix of the cubic spline inverse CSD"
function f_cubicspline(el_pos,r,cond=0.3,cond_top=cond)
  N = length(el_pos)
  h = mean(diff(el_pos))
  z_js = vcat(0,el_pos,el_pos[end]+h) # electrode positions with two imaginary

  E0,E1,E2,E3 = e_matrixes(el_pos)
  # Potential functions
  f0(zeta,zj,diam,cond) = 1/(2cond)*(sqrt((diam/2)^2 + (zj-zeta)^2) - abs(zj-zeta))
  f1(zeta,zj,zi,diam,cond) = (zeta-zi)*f0(zeta,zj,diam,cond)
  f2(zeta,zj,zi,diam,cond) = (zeta-zi)^2*f0(zeta,zj,diam,cond)
  f3(zeta,zj,zi,diam,cond) = (zeta-zi)^3*f0(zeta,zj,diam,cond)
  # Define integration matrixes
  F0  = zeros(N,N+1)
  F1  = zeros(N,N+1)
  F2  = zeros(N,N+1)
  F3  = zeros(N,N+1)

  for j = 1:N
    for i = 1:N+1
      F0[j,i] = hquadrature(x->f0(x,z_js[j+1],r,cond),z_js[i],z_js[i+1])[1]
      F1[j,i] = hquadrature(x->f1(x,z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
      F2[j,i] = hquadrature(x->f2(x,z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
      F3[j,i] = hquadrature(x->f3(x,z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
      if cond != cond_top # image technique if conductivity not constant
        F0[j,i] = F0[j,i] + (cond-cond_top)/(cond+cond_top)*hquadrature(x->f0(x,-z_js[j+1],r,cond),z_js[i],z_js[i+1])[1]
        F1[j,i] = F1[j,i] + (cond-cond_top)/(cond+cond_top)*hquadrature(x->f1(x,-z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
        F2[j,i] = F2[j,i] + (cond-cond_top)/(cond+cond_top)*hquadrature(x->f2(x,-z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
        F3[j,i] = F3[j,i] + (cond-cond_top)/(cond+cond_top)*hquadrature(x->f3(x,-z_js[j+1],z_js[i],r,cond),z_js[i],z_js[i+1])[1]
      end
    end
  end

  temp_F = F0*E0+F1*E1+F2*E2+F3*E3  # F matrix(N x N+2)
  # Convert to (N+2 x N+2) matrixes by applying the boundary I_0 = I_N+1 = 0.
  F = zeros(N+2,N+2)
  F[2:N+1,:] = temp_F
  F[1,1] = 1         # implies I_1 = ϕ_1
  F[N+2,N+2] = 1     # implies I_N+2 = ϕ_N+2
  return F,E0,E1,E2,E3
end

ext: jl


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

using Cxx

export TTreeBindings
export TTreeInput
export TTreeOutput
export TChainInput

export ttree_binding_proxy

struct TTreeBindings
    values::Dict{Symbol, Any}
    direct_values::Dict{Symbol, Any}
    proxied_values::Dict{Symbol, Pair{Any, Any}}

    TTreeBindings() = new(
        Dict{Symbol, Any}(),
        Dict{Symbol, Any}(),
        Dict{Symbol, Pair{Any, Any}}()
    )
end

Base.getindex(bindings::TTreeBindings, key::Symbol) = bindings.values[key]

Base.setindex!(bindings::TTreeBindings, value::Any, key::Symbol) = begin
    bindings.values[key] = value

    proxy = ttree_binding_proxy(value)
    if proxy == nothing
        # info("Binding branch name $(key) to value of type $(typeof(value))")
        haskey(bindings.proxied_values, key) && delete!(bindings.proxied_values, key)
        bindings.direct_values[key] = value
    else
        # info("Binding branch name  $(key) to value of type $(typeof(value)) via proxy of type $(typeof(proxy))")
        haskey(bindings.direct_values, key) && delete!(bindings.direct_values, key)
        bindings.proxied_values[key] = Pair{Any, Any}(value, proxy)
    end

    bindings
end

copy_to_proxies!(bindings::TTreeBindings) = begin
    for (_, (value, proxy_value)) in bindings.proxied_values
        copy!(proxy_value, value)
    end
end

copy_from_proxies!(bindings::TTreeBindings) = begin
    for (_, (value, proxy_value)) in bindings.proxied_values
        copy!(value, proxy_value)
    end
end

bind_branches!(ttree::ATTreeInst, bindings::TTreeBindings) = begin
    for (key, value) in bindings.direct_values
        # info("Binding branch $(string(key)) of type $(typeof(value)) in $(ttree)")
        bind_branch!(ttree, string(key), value)
    end
    for (key, (_, proxy)) in bindings.proxied_values
        # info("Binding branch $(string(key)) of type $(typeof(proxy)) in $(ttree)")
        bind_branch!(ttree, string(key), proxy)
    end
end

create_branches!(ttree::ATTreeInst, bindings::TTreeBindings) = begin
    for (key, value) in bindings.direct_values
        # info("Creating branch $(string(key)) of type $(typeof(value)) in $(ttree)")
        create_branch!(ttree, string(key), value)
    end
    for (key, (_, proxy)) in bindings.proxied_values
        # info("Creating branch $(string(key)) of type $(typeof(proxy)) in $(ttree)")
        create_branch!(ttree, string(key), proxy)
    end
end

ttree_binding_proxy(value::Ref{T}) where T = nothing

struct CxxObjWithPtrRef{T,N,CVR}
    x::CppValue{T,N}
    ptrref::Ref{CppPtr{T,CVR}}
end

CxxObjWithPtrRef(x::CppValue) = CxxObjWithPtrRef(x, Ref(icxx"&$x;"))

bind_branch!(ttree::ATTreeInst, name::AbstractString, value::CxxObjWithPtrRef) =
   bind_branch!(ttree, name, value.ptrref)

create_branch!(ttree::ATTreeInst, name::AbstractString, value::CxxObjWithPtrRef; kwargs...) =
   create_branch!(ttree, name, value.ptrref, kwargs...)

ttree_binding_proxy(value::CxxObjWithPtrRef) = nothing

Base.copy!(a, b::CxxObjWithPtrRef) = copy!(a, b.x)
Base.copy!(a::CxxObjWithPtrRef, b) = copy!(a.x, b)
Base.copy!(a::CxxObjWithPtrRef, b::CxxObjWithPtrRef) = copy!(a.x, b.x)

ttree_binding_proxy(value::AbstractVector{T}) where T = CxxObjWithPtrRef(icxx"std::vector<$T>();")
Base.copy!(a::AbstractVector, b::CxxObjWithPtrRef) = (resize!(a, length(b.x)); copy!(a, b.x))
Base.copy!(a::CxxObjWithPtrRef, b::AbstractVector) = (resize!(a.x, length(b)); copy!(a.x, b))

struct TTreeInput <: AbstractVector{Nothing}
    ttree::ATTreeInst
    bindings::TTreeBindings

    TTreeInput(ttree::ATTreeInst, bindings::TTreeBindings) = begin
        input = new(ttree, bindings)
        bind_branches!(ttree, bindings)
        input
    end
end

Base.show(io::IO, x::TTreeInput) = show(io, MIME"text/plain"(), x)
Base.show(io::IO, ::MIME"text/plain", x::TTreeInput) = Base.show_default(io, x)

Base.size(input::TTreeInput) = (length(input.ttree),)

Base.getindex(input::TTreeInput, i::Integer) = begin
    result = getindex(input.ttree, i)
    copy_from_proxies!(input.bindings)
    result
end

mutable struct TChainInput <: AbstractVector{Nothing}
    tchain::TChainPtr
    bindings::TTreeBindings
end

function Base.open(::Type{TChainInput}, bindings::TTreeBindings, treename::AbstractString, filenames::AbstractString...; wildcards::Bool = false)
    tchain = lock(gROOTMutex()) do
        icxx"new TChain($treename);"
    end

    if wildcards
        for filename in filenames
            @cxx tchain->Add(pointer(filename))
        end
    else
        for filename in filenames
            @cxx tchain->AddFile(pointer(filename))
        end
    end

    @cxx tchain->SetCacheSize(-1);
    @cxx tchain->SetBranchStatus(pointer("*"), false);
    bind_branches!(tchain, bindings)

    result = TChainInput(tchain, bindings)
    finalizer(close, result)
    result
end

Base.show(io::IO, x::TChainInput) = show(io, MIME"text/plain"(), x)
Base.show(io::IO, ::MIME"text/plain", x::TChainInput) = Base.show_default(io, x)

Base.open(::Type{TChainInput}, treename::AbstractString, filenames::AbstractString...; kwargs...) =
    open(TChainInput, TTreeBindings(), treename, filenames...; kwargs...)

Base.isopen(input::TChainInput) = (input.tchain != C_NULL)

function Base.close(input::TChainInput)
    if isopen(input)
        tchain = lock(gROOTMutex()) do
            icxx"delete $(input.tchain);"
        end
        input.tchain = icxx"(TChain*)(nullptr);"
    end
    nothing
end

Base.size(input::TChainInput) = (length(input.tchain),)

Base.getindex(input::TChainInput, i::Integer) = begin
    result = getindex(input.tchain, i)
    copy_from_proxies!(input.bindings)
    result
end

struct TTreeOutput
    ttree::TTreePtr
    bindings::TTreeBindings

    TTreeOutput(ttree::TTreePtr, bindings::TTreeBindings) = begin
        output = new(ttree, bindings)
        create_branches!(ttree, bindings)
        output
    end
end

push!(output::TTreeOutput) = begin
    copy_to_proxies!(output.bindings)
    push!(output.ttree)
    nothing
end


text:
using Test
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
import Random

cfg = get_config("test.yaml", nin=4, nout=1)
mutate(i::CGPInd) = uniform_mutate(cfg, i)

function rosenbrock(x::Array{Float64})
    sum([(1.0 - x[i])^2 + 100.0 * (x[i+1] - x[i]^2)^2
         for i in 1:(length(x)-1)])/200
end

function symbolic_evaluate(i::CGPInd; seed::Int=0)
    Random.seed!(seed)
    inputs = rand(i.n_in)
    output = process(i, inputs)
    target = rosenbrock(inputs)
    [-(output[1] - target)^2]
end

@testset "Symbolic Regression Evolution" begin
    e = CGPEvolution(cfg, symbolic_evaluate)

    step!(e)
    @test length(e.population) == cfg.n_population
    best = sort(e.population)[end]
    @test best.fitness[1] <= 0.0
    @test e.gen == 1

    run!(e)
    @test length(e.population) == cfg.n_population
    @test e.gen == cfg.n_gen
    println("Evolution step for symbolic regression")
    @timev step!(e)
    new_best = sort(e.population)[end]
    println("Final fitness: ", new_best.fitness[1])
    @test new_best.fitness[1] <= 0.0
end


text:
module BasnTest
using SimplePNGs
using Test

include("common.jl")
using .TestCommon: load_json

pl(name) = SimplePNGs.load(joinpath("PngSuite", name*".png"))

@testset "Basic format test files (non-interlaced)" begin
    @testset "black & white" begin
        img1 = load_json("basn0g01")
        img2 = pl("basn0g01")
        @test img1 == img2
    end

    @testset "2 bit (4 level) grayscale" begin
        img1 = load_json("basn0g02")
        img2 = pl("basn0g02")
        @test img1 == img2
    end

    @testset "4 bit (16 level) grayscale" begin
        img1 = load_json("basn0g04")
        img2 = pl("basn0g04")
        @test img1 == img2
    end

    @testset "8 bit (256 level) grayscale" begin
        img1 = load_json("basn0g08")
        img2 = pl("basn0g08")
        @test img1 == img2
    end

    @testset "16 bit (64k level) grayscale" begin
        img1 = load_json("basn0g16")
        img2 = pl("basn0g16")
        @test img1 == img2
    end

    @testset "3x8 bits rgb color" begin
        img1 = load_json("basn2c08")
        img2 = pl("basn2c08")
        @test img1 == img2
    end

    @testset "3x16 bits rgb color" begin
        img1 = load_json("basn2c16")
        img2 = pl("basn2c16")
        @test img1 == img2
    end

    @testset "1 bit (2 color) paletted" begin
        img1 = load_json("basn3p01")
        img2 = pl("basn3p01")
        @test img1 == img2
    end

    @testset "2 bit (4 color) paletted" begin
        img1 = load_json("basn3p02")
        img2 = pl("basn3p02")
        @test img1 == img2
    end

    @testset "4 bit (16 color) paletted" begin
        img1 = load_json("basn3p04")
        img2 = pl("basn3p04")
        @test img1 == img2
    end

    @testset "8 bit (256 color) paletted" begin
        img1 = load_json("basn3p08")
        img2 = pl("basn3p08")
        @test img1 == img2
    end

    @testset "8 bit grayscale + 8 bit alpha-channel" begin
        img1 = load_json("basn4a08")
        img2 = pl("basn4a08")
        @test img1 == img2
    end

    @testset "16 bit grayscale + 16 bit alpha-channel" begin
        img1 = load_json("basn4a16")
        img2 = pl("basn4a16")
        @test img1 == img2
    end

    @testset "3x8 bits rgb color + 8 bit alpha-channel" begin
        img1 = load_json("basn6a08")
        img2 = pl("basn6a08")
        @test img1 == img2
    end

    @testset "3x16 bits rgb color + 16 bit alpha-channel" begin
        img1 = load_json("basn6a16")
        img2 = pl("basn6a16")
        @test img1 == img2
    end
end

end # module


text:
"""
# Description

Given a tuple of symbols, return a `Dict{Symbol, Int}`` mapping each symbol to
its index in the input tuple.

# Arguments

1. `column_names::NTuple{N, Symbol} where N`: A tuple of symbols.

# Return Values

1. `index::Dict{Symbol, Int}`: A dictionary mapping symbols to their indices in
  the input tuple.

# Examples

```
julia> index_column_names(((:a, false), (:b, false)))
Dict{ColumnName,Int64} with 2 entries:
  (:a, false) => 1
  (:b, false) => 2
```
"""
function index_column_names(column_names::NTuple{N, ColumnName}) where N
    Dict{ColumnName, Int}(column_names .=> 1:length(column_names))
end


text:
using Test
using SafeTestsets

@safetestset "TightBindingApproximation" begin
    include("TightBindingApproximation.jl")
end


text:
#
# Conjugate gradient
#
# This is an independent implementation of:
#   W. W. Hager and H. Zhang (2006) Algorithm 851: CG_DESCENT, a
#     conjugate gradient method with guaranteed descent. ACM
#     Transactions on Mathematical Software 32: 113–137.
#
# Code comments such as "HZ, stage X" or "HZ, eqs Y" are with
# reference to a particular point in this paper.
#
# Several aspects of the following have also been incorporated:
#   W. W. Hager and H. Zhang (2012) The limited memory conjugate
#     gradient method.
#
# This paper will be denoted HZ2012 below.
#
# There are some modifications and/or extensions from what's in the
# paper (these may or may not be extensions of the cg_descent code
# that can be downloaded from Hager's site; his code has undergone
# numerous revisions since publication of the paper):
#
# cgdescent: the termination condition employs a "unit-correct"
#   expression rather than a condition on gradient
#   components---whether this is a good or bad idea will require
#   additional experience, but preliminary evidence seems to suggest
#   that it makes "reasonable" choices over a wider range of problem
#   types.
#
# linesearch: the Wolfe conditions are checked only after alpha is
#   generated either by quadratic interpolation or secant
#   interpolation, not when alpha is generated by bisection or
#   expansion. This increases the likelihood that alpha will be a
#   good approximation of the minimum.
#
# linesearch: In step I2, we multiply by psi2 only if the convexity
#   test failed, not if the function-value test failed. This
#   prevents one from going uphill further when you already know
#   you're already higher than the point at alpha=0.
#
# both: checks for Inf/NaN function values
#
# both: support maximum value of alpha (equivalently, c). This
#   facilitates using these routines for constrained minimization
#   when you can calculate the distance along the path to the
#   disallowed region. (When you can't easily calculate that
#   distance, it can still be handled by returning Inf/NaN for
#   exterior points. It's just more efficient if you know the
#   maximum, because you don't have to test values that won't
#   work.) The maximum should be specified as the largest value for
#   which a finite value will be returned.  See, e.g., limits_box
#   below.  The default value for alphamax is Inf. See alphamaxfunc
#   for cgdescent and alphamax for linesearch_hz.

immutable ConjugateGradient{T, Tprep<:Union{Function, Void}, L<:Function} <: Optimizer
    eta::Float64
    P::T
    precondprep!::Tprep
    linesearch!::L
end
#= uncomment for v0.8.0
function ConjugateGradient(;
                           linesearch = LineSearches.hagerzhang!,
                           eta::Real = 0.4,
                           P::Any = nothing,
                           precondprep = (P, x) -> nothing)
    ConjugateGradient(Float64(eta),
                                 P, precondprep,
                                 linesearch)
end
=#
function ConjugateGradient(; linesearch! = nothing,
                             linesearch = LineSearches.hagerzhang!,
                             eta::Real = 0.4,
                             P::Any = nothing,
                             precondprep! = nothing,
                             precondprep = (P, x) -> nothing)

    linesearch = get_linesearch(linesearch!, linesearch)
    precondprep = get_precondprep(precondprep!, precondprep)
    ConjugateGradient(Float64(eta),
                      P, precondprep,
                      linesearch)
end

type ConjugateGradientState{T}
    @add_generic_fields()
    x_previous::Array{T}
    g::Array{T}
    g_previous::Array{T}
    f_x_previous::T
    y::Array{T}
    py::Array{T}
    pg::Array{T}
    s::Array{T}
    @add_linesearch_fields()
end

function initial_state{T}(method::ConjugateGradient, options, d, initial_x::Array{T})
    g = similar(initial_x)
    f_x = d.fg!(initial_x, g)
    pg = copy(g)
    @assert typeof(f_x) == T
    # Output messages
    if !isfinite(f_x)
        error("Must have finite starting value")
    end
    if !all(isfinite, g)
        @show g
        @show find(!isfinite.(g))
        error("Gradient must have all finite values at starting point")
    end

    # Determine the intial search direction
    #    if we don't precondition, then this is an extra superfluous copy
    #    TODO: consider allowing a reference for pg instead of a copy
    method.precondprep!(method.P, initial_x)
    A_ldiv_B!(pg, method.P, g)

    ConjugateGradientState("Conjugate Gradient",
                         length(initial_x),
                         copy(initial_x), # Maintain current state in state.x
                         f_x, # Store current f in state.f_x
                         1, # Track f calls in state.f_calls
                         1, # Track g calls in state.g_calls
                         0, # Track h calls in state.h_calls
                         copy(initial_x), # Maintain current state in state.x_previous
                         g, # Store current gradient in state.g
                         copy(g), # Store previous gradient in state.g_previous
                         T(NaN), # Store previous f in state.f_x_previous
                         similar(initial_x), # Intermediate value in CG calculation
                         similar(initial_x), # Preconditioned intermediate value in CG calculation
                         pg, # Maintain the preconditioned gradient in pg
                         -copy(pg), # Maintain current search direction in state.s
                         @initial_linesearch()...) # Maintain a cache for line search results in state.lsr
end

function update_state!{T}(df, state::ConjugateGradientState{T}, method::ConjugateGradient)
        lssuccess = true
        # Reset the search direction if it becomes corrupted
        dphi0 = vecdot(state.g, state.s)
        if dphi0 >= 0
            @simd for i in 1:state.n
                @inbounds state.s[i] = -state.pg[i]
            end
            dphi0 = vecdot(state.g, state.s)
            if dphi0 >= 0
                return true
            end
        end

        # Refresh the line search cache
        LineSearches.clear!(state.lsr)
        @assert typeof(state.f_x) == T
        @assert typeof(dphi0) == T
        push!(state.lsr, zero(T), state.f_x, dphi0)

        # Pick the initial step size (HZ #I1-I2)
        state.alpha, state.mayterminate, f_update, g_update =
          LineSearches.alphatry(state.alpha, df, state.x, state.s, state.x_ls, state.g_ls, state.lsr)
        state.f_calls, state.g_calls = state.f_calls + f_update, state.g_calls + g_update

        # Determine the distance of movement along the search line
        try
            state.alpha, f_update, g_update =
                method.linesearch!(df, state.x, state.s, state.x_ls, state.g_ls, state.lsr,
                                   state.alpha, state.mayterminate)
            state.f_calls, state.g_calls = state.f_calls + f_update, state.g_calls + g_update
        catch ex
            if isa(ex, LineSearches.LineSearchException)
                lssuccess = false
                state.f_calls, state.g_calls = state.f_calls + ex.f_update, state.g_calls + ex.g_update
                state.alpha = ex.alpha
                Base.warn("Linesearch failed, using alpha = $(state.alpha) and exiting optimization.")
            else
                rethrow(ex)
            end
        end

        # Maintain a record of previous position
        copy!(state.x_previous, state.x)

        # Update current position # x = x + alpha * s
        LinAlg.axpy!(state.alpha, state.s, state.x)

        # Maintain a record of the previous gradient
        copy!(state.g_previous, state.g)

        # Update the function value and gradient
        state.f_x_previous, state.f_x = state.f_x, df.fg!(state.x, state.g)
        state.f_calls, state.g_calls = state.f_calls + 1, state.g_calls + 1

        # Check sanity of function and gradient
        if !isfinite(state.f_x)
            error("Function must finite function values")
        end

        # Determine the next search direction using HZ's CG rule
        #  Calculate the beta factor (HZ2012)
        # -----------------
        # Comment on py: one could replace the computation of py with
        #    ydotpgprev = vecdot(y, pg)
        #    vecdot(y, py)  >>>  vecdot(y, pg) - ydotpgprev
        # but I am worried about round-off here, so instead we make an
        # extra copy, which is probably minimal overhead.
        # -----------------
        method.precondprep!(method.P, state.x)
        dPd = dot(state.s, method.P, state.s)
        etak::T = method.eta * vecdot(state.s, state.g_previous) / dPd
        @simd for i in 1:state.n
            @inbounds state.y[i] = state.g[i] - state.g_previous[i]
        end
        ydots = vecdot(state.y, state.s)
        copy!(state.py, state.pg)        # below, store pg - pg_previous in py
        A_ldiv_B!(state.pg, method.P, state.g)
        @simd for i in 1:state.n     # py = pg - py
           @inbounds state.py[i] = state.pg[i] - state.py[i]
        end
        betak = (vecdot(state.y, state.pg) - vecdot(state.y, state.py) * vecdot(state.g, state.s) / ydots) / ydots
        beta = max(betak, etak)
        @simd for i in 1:state.n
            @inbounds state.s[i] = beta * state.s[i] - state.pg[i]
        end
        (lssuccess == false) # break on linesearch error
end

update_g!(d, state, method::ConjugateGradient) = nothing


text:
@testset "Simplex" begin
    s = Simplex(1, 2, 3)
    @test s == Simplex(1, 2, 3)
    @test sort!(values(s)) == [1,2,3]

    s = Simplex([1,2,3])
    @test dim(s) == 2
    @test sort!(values(s)) == [1,2,3]
    allfaces = [Simplex(2, 3), Simplex(1, 3), Simplex(1, 2)]
    @testset "Simplex faces" for o in faces(s)
        idx = findfirst(isequal(o), allfaces)
        @test values(o) == values(allfaces[idx])
    end

    ss = Set{Simplex}()
    push!(ss, s)
    push!(ss, s)
    @test length(ss) == 1

    s = Simplex([1,2], [1,3])
    @test dim(s) == 1

    X = hcat(zeros(3), Matrix(I,3,3))
    S = Simplex(1,2,3,4)
    @test volume(S,X) == 1/6

    @test faces(s)[1] ∪ faces(s)[2] == s
end

@testset "Cube" begin
    c = Cube{Vector{Int}}([0,0,0],Int[]) # 0-cube
    @test dim(c) == 0
    c = Cube{Vector{Int}}([0,0,0],Int[0,0,0]) # 0-cube
    @test dim(c) == 0
    c = Cube{Vector{Int}}([0,0,0],[1]) # 1-cube
    @test dim(c) == 1
    c = Cube{Vector{Int}}([0,0,0],[1,0,1]) # 2-cube
    @test dim(c) == 2
    c = Cube{Vector{Int}}([0,0,0],[1,1,1]) # 3-cube
    @test dim(c) == 3
    @test hash(c) == hash([0,0,0,1,1,1])
    @test c == Cube{Vector{Int}}([0,0,0],[1,1,1])

    c = Cube{Vector{Float64}}([0.,0.,0.],[.5,.0,.5]) # 2-cube
    @test volume(c) == 0.5^2
    c = Cube{Vector{Float64}}([0.,0.,0.],[.5,.5,.5]) # 2-cube
    @test volume(c) == 0.5^3

end

@testset "CW" begin
    a = Cell()
    @test hash(a) == hash(Cell[])
    @test a == Cell()
    @test dim(a) == 0

    b = Cell()
    c = Cell(1, a, b)
    @test dim(c) == 1
    d = Cell(1, a, b, a, b)
    e = c ∪ d
    @test dim(e) == 2

    @testset for (f1, f2) in zip(faces(c), [a, b])
        @test f1 == f2
    end

    @testset for (f1, f2) in zip(faces(e), [c, d])
        @test f1 == f2
    end
end


text:
#import Base.LinAlg.BLAS: nrm2

for (func, elty) in [(:CLBlastSnrm2, Float32), (:CLBlastDnrm2, Float64),
                     (:CLBlastScnrm2, ComplexF32), (:CLBlastDznrm2, ComplexF64)]
    #TODO: (:CLBlastHnrm2, Float16)

    @eval function $func(n::Integer, out_buffer::cl.CL_mem, out_offset::Integer,
                         x_buffer::cl.CL_mem, x_offset::Integer, x_inc::Integer,
                         queue::cl.CmdQueue, event::cl.Event)
        err = ccall(
            ($(string(func)), libCLBlast),
            cl.CL_int,
            (Csize_t, Ptr{Cvoid}, Csize_t, Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}),
            n, out_buffer, out_offset, x_buffer, x_offset, x_inc, Ref(queue), Ref(event)
        )
        if err != cl.CL_SUCCESS
            println(stderr, "Calling function $(string($func)) failed!")
            throw(cl.CLError(err))
        end
        return err
    end

    @eval function nrm2(n::Integer, x::cl.CLArray{$elty}, x_inc::Integer;
                        queue::cl.CmdQueue=cl.queue(x))
        # output buffer and event
        ctx = cl.context(queue)
        out = zeros($elty, 1)
        out_buffer = cl.Buffer($elty, ctx, (:rw, :copy), hostbuf=out)
        event::cl.Event = cl.Event(C_NULL)

        $func(Csize_t(n), pointer(out_buffer), Csize_t(0), pointer(x), Csize_t(0), Csize_t(x_inc),
              queue, event)

        # wait for kernel and read return value
        cl.wait(event)
        cl.enqueue_read_buffer(queue, out_buffer, out, Csize_t(0), nothing, true)

        return real(first(out))
    end

end


text: Dict{String, Any}("utc2" => Dict{String, Any}("value" => "1987-07-05T17:45:56.600000Z", "type" => "datetime"), "utc1" => Dict{String, Any}("value" => "1987-07-05T17:45:56.123456Z", "type" => "datetime"), "wita1" => Dict{String, Any}("value" => "1987-07-05T17:45:56.123456+08:00", "type" => "datetime"), "wita2" => Dict{String, Any}("value" => "1987-07-05T17:45:56.600000+08:00", "type" => "datetime"))


text:
using BinDeps
using Compat

@BinDeps.setup

function validate_libarchive(name, handle)
    try
        p = Libdl.dlsym(handle, :archive_read_free)
        return p != C_NULL
    catch
        return false
    end
end

libarchive = library_dependency("libarchive",
                                aliases=["libarchive", "libarchive-12"],
                                validate=validate_libarchive)

if is_linux()
    provides(Pacman, "libarchive", libarchive)
    provides(AptGet, "libarchive12", libarchive)
    provides(Yum, "libarchive", libarchive)
end

if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "libarchive12", libarchive, os=:Windows)
end

if is_apple()
    using Homebrew
    provides(Homebrew.HB, "libarchive", libarchive, os=:Darwin)
end

@BinDeps.install Dict(:libarchive => :libarchive)


text:
module LinearAlgebra

using Reexport

include("GLA.jl")
include("StructuredGLA.jl")

@reexport using .GraphicalLinearAlgebra
@reexport using .StructuredGraphicalLinearAlgebra

end


text:
## --- Test `MallocString`s

    # Test MallocString constructors
    str = m"Hello, world! 🌍"
    @test isa(str, MallocString)
    @test sizeof(str) == 19
    @test StaticTools.strlen(str) == length(str) == 18
    @test pointer(str) == Base.unsafe_convert(Ptr{UInt8}, str)
    str1 = MallocString(c"Hello, world! 🌍")
    @test isa(str1, MallocString)

    # Test basic string operations
    @test str == str1
    m, p = str*str, str^2
    @test m == p
    @test free(m) == 0
    @test free(p) == 0
    @test codeunit(str) === UInt8
    @test codeunit(str, 5) == UInt8('o')
    @test ncodeunits(str) == length(str)+1
    @test codeunits(str) == codeunits(c"Hello, world! 🌍")
    @test codeunits(c"Hello, world! 🌍") == codeunits(str)

    # Test mutability
    str[8] = 'W'
    @test str[8] == 0x57 # W
    str[:] = c"Hello, world! 🌍"
    @test str[8] == 0x77 # w

    # Test indexing
    @test isa(str[1:end], StringView)
    @test str == str[1:end]
    @test str === str[:]
    @test str[1:2] === str[1:2]
    @test str[1:2] != str[1:3]
    strc = copy(str)
    @test strc == str
    @test free(str) == 0
    @test free(str1) == 0
    @test free(strc) == 0

    # Test ascii escaping
    many_escapes = m"\"\0\a\b\f\n\r\t\v\'\"\\"
    @test isa(many_escapes, MallocString)
    @test length(many_escapes) == 12
    @test codeunits(many_escapes) == codeunits("\"\0\a\b\f\n\r\t\v'\"\\\0")

    # Test unsafe_mallocstring
    s = "Hello there!"
    m = unsafe_mallocstring(pointer(s))
    @test isa(m, MallocString)
    @test length(m) == 12
    @test codeunits(m) == codeunits(c"Hello there!")
    @test free(m) == 0

    # Test constructing from Ptr{Ptr{UInt8}} as in argv
    s1,s2 = "Hello", "there"
    a = [pointer(s1), pointer(s2)]
    argv = pointer(a)
    @test MallocString(argv,1) == c"Hello"
    @test MallocString(argv,2) == c"there"

    # Test consistency with other string types
    abc = m"abc"
    @test abc == "abc"
    @test "abc" == abc
    @test abc == c"abc"
    @test abc == abc[1:3]
    @test abc[1:3] == "abc"
    free(abc)


text:
# Setup for example.
include("KronGraph500NoPerm.jl")
using Random

# Define SCALE and EdgesPerVertex if not already defined
(@isdefined SCALE) || (SCALE = 12)
(@isdefined EdgesPerVertex) || (EdgesPerVertex = 16)

# Calculate number of vertices and edges
Nmax = 2 ^ SCALE
M = EdgesPerVertex * Nmax

# Create directory to write output
data_dir = joinpath(Base.source_dir(),"data/")
isdir(data_dir) || mkdir(data_dir)

# Create graphs and save in files
Nfile = 8
for i = 1:Nfile
    #tic()
    pTime = @elapsed begin
        # Create filename
        fname = data_dir * string(i)
        println("On file: "*string(i))

        #Reset Rand Seed
        Random.seed!(i)

        # Generate data
        v1,v2 = KronGraph500NoPerm(SCALE,round(Int,EdgesPerVertex/Nfile))

        # Convert to strings
        rowStr = join(string.(v1),'\n')
        colStr = join(string.(v2),'\n')
        valStr = join(string.(ones(Int,length(v1))),'\n')

        # Open files, write data, and close files
        fidRow = open(fname * "r.txt","w")
        fidCol = open(fname * "c.txt","w")
        fidVal = open(fname * "v.txt","w")

        write(fidRow,rowStr)
        close(fidRow)

        write(fidCol,colStr)
        close(fidCol)

        write(fidVal,valStr)
        close(fidVal)
    end

    #pTime = toq()
    println("Time: ", pTime)
    println("Edges/sec: ", string(length(v1) / pTime))
end


text:
export buildsif, servertransfer

# todo: as soon as there are option for file generation, take same options here
# could possibly make g take arguments to be interpolated?
function buildsif(;verbose = false, force = true)
    ppath = dirname(Base.active_project(false))

    if force
        cmd = `sudo singularity build --force container/projectcontainer.sif container/Singularity.pack`
    elseif !verbose
        @warn("Image already exists. To overwrite run verbose = true, or force = true")
        return
    else
        cmd = `sudo singularity build container/projectcontainer.sif container/Singularity.pack`
    end

    g() = verbose ? run(cmd) : run(pipeline(cmd, stdout="container/pack.log"))

    try
        println("running $cmd")
        cd(g, ppath)
    catch
        # if !verbose
        @warn("Build failed")
    end

end

# ToDo: add 'prefix' to
function servertransfer(host)
    ppath = dirname(Base.active_project(false))
    pfolder = basename(ppath)

    mkfolder = `ssh greyplover.stats mkdir -p $pfolder/data`
    run(mkfolder)

    transfer = `rsync -azP container/projectcontainer.sif $host:$pfolder`
    g() = run(transfer)
    cd(g, ppath)
end


text:
using LinearAlgebra
import Base: *
import LinearAlgebra: mul!

mutable struct BroydenOperator{R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}}
    H
    theta_bar::R
    function BroydenOperator{R,C,T}(x::T, H, theta_bar) where {R,C,T}
        new(H, theta_bar)
    end
end

BroydenOperator(
    x::T;
    H = I,
    theta_bar = R(0.2),
) where {R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}} =
    BroydenOperator{R,C,T}(x, H, theta_bar)

_sign(x::R) where {R} = x == 0 ? R(1) : sign(x)

function update!(L::BroydenOperator{R,C,T}, s, y) where {R,C,T}
    Hy = L.H * y
    sH = s' * L.H
    delta = dot(Hy, s) / norm(s)^2
    theta = if abs(delta) >= L.theta_bar
        R(1)
    else
        (1 - _sign(delta) * L.theta_bar) / (1 - delta)
    end
    L.H += (s - Hy) / dot(s, (1 / theta - 1) * s + Hy) * sH
end

function reset!(L::BroydenOperator{R,C,T}) where {R,C,T}
    L.H = I
    L.theta_bar = R(0.2)
end

function (*)(L::BroydenOperator, v)
    w = similar(v)
    return mul!(w, L, v)
end

function mul!(d::T, L::BroydenOperator, v::T) where {T}
    mul!(d, L.H, v)
    return d
end

Base.@kwdef struct Broyden{R}
    theta_bar::R = 0.2
end

acceleration_style(::Type{<:Broyden}) = QuasiNewtonStyle()

function initialize(broyden::Broyden, x)
    return BroydenOperator(x, theta_bar=broyden.theta_bar)
end


text:
# define imperial units - from http://en.wikipedia.org/w/index.php?title=Imperial_units&oldid=574336154
# Length
const Inch = DerivedUnit("in", 254//10000*Meter)
const Foot = DerivedUnit("ft", 12*Inch)
const Yard = DerivedUnit("yd", 0.9144*Meter) # wikipedia says this is exact
const Furlong = DerivedUnit("fur", 220*Yard)
const Chain = DerivedUnit("ch", 22*Yard)
const Mile = DerivedUnit("mi", 1760*Yard)
const League = DerivedUnit("lea", 3*Mile)
const Fathom = DerivedUnit("ftm", 1.853184*Meter)

# Area
const Perch = 1//102400*Mile^2
const Rood = 1//2560*Mile^2
const Acre = 1//640*Mile^2

# Mass + Weight
const Pound = DerivedUnit("lb",4.44822162*Newton)
const Slug = DerivedUnit("slug", 14.5939029*KiloGram)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e90001ca-4b3f-11eb-0ccd-7785902a32e3
using DrWatson

# ╔═╡ 6e33e012-4b3d-11eb-30ea-2b5f82c16b8e
begin
	@quickactivate

	using PlaymatSimulator
	using Serialization
	using GameZero
	using PlutoUI
	using Colors

	md"""
	# 🧙 Welcome to Mage.jl
	#### Choose a game variant below to get started ☮
	"""
end

# ╔═╡ 65e7c338-4b42-11eb-014a-017404d2040a
begin
	GAMES = ["EDH"=>"EDH"]
	#[ gn=>gn for gn in readdir(GAME_DIR) if !(occursin(".", gn)) ]

	pushfirst!(GAMES,""=>"")

	md"""
	Select a game to play: $( @bind GAME_NAME Select(GAMES) )
	"""
end

# ╔═╡ e86a6936-4b43-11eb-33b9-3ddefed762c7
if GAME_NAME != ""
	DECKS = [ gn=>gn for gn in readdir("$GAME_NAME/decks") if !(occursin(".", gn)) ]
	pushfirst!(DECKS,""=>"")

	md"""
	Select a deck to play: $( @bind DECK_NAME Select(DECKS) )
	"""
end

# ╔═╡ eb1a3e94-4cca-11eb-2b39-a5794006c375
if GAME_NAME != "" && DECK_NAME != ""

	US = Dict(
		:GAME_NAME => GAME_NAME,
		:DECK_NAME => DECK_NAME,
	)

	serialize("tmp/user_selection.jls", US)
end

# ╔═╡ e11fc3c8-4b44-11eb-0c25-3d787ef9da0a
if GAME_NAME != "" && DECK_NAME != ""
	md"""
	###### Simulating an $GAME_NAME game with a $DECK_NAME deck, please wait...

	*press alt+F4 to quit*
	"""
end

# ╔═╡ a6cc0f5c-4cab-11eb-10d1-b351e17004e3
if (@isdefined DECK_NAME) && !(DECK_NAME == "")
	rungame("$GAME_NAME/notebooks/$GAME_NAME.jl")
end

# ╔═╡ 386f2d76-4ef1-11eb-032b-1715050b4952
md"""
*Want to develop your own PlaymatSimulator games? Check out the docs!
"""

# ╔═╡ 1816db80-4b62-11eb-0a61-91f0662ae287
html"""

"""

# ╔═╡ Cell order:
# ╟─e90001ca-4b3f-11eb-0ccd-7785902a32e3
# ╟─6e33e012-4b3d-11eb-30ea-2b5f82c16b8e
# ╟─65e7c338-4b42-11eb-014a-017404d2040a
# ╟─e86a6936-4b43-11eb-33b9-3ddefed762c7
# ╟─eb1a3e94-4cca-11eb-2b39-a5794006c375
# ╟─e11fc3c8-4b44-11eb-0c25-3d787ef9da0a
# ╟─a6cc0f5c-4cab-11eb-10d1-b351e17004e3
# ╟─386f2d76-4ef1-11eb-032b-1715050b4952
# ╟─1816db80-4b62-11eb-0a61-91f0662ae287


text:
using BenchmarkTools
using Test

import Random
Random.seed!(321)

# Distance function, considering PBC
function r(x,y,side)
  d = 0.
  for i in 1:2
    dx = (x[i]-y[i])%side
    if dx > side/2
      dx = dx - side
    elseif  dx < -side/2
      dx = dx + side
    end
    d = d + dx^2
  end
  return sqrt(d)
end

# Potential energy function
U(r,eps,sig) = 4*eps*((sig/r)^12 - (sig/r)^6)

# Compute total potential energy for all pairs
function UT1(x,eps,sig,side)
  ut = 0.
  for i in 1:length(x)-1, j in i+1:length(x)
    rij = r(x[i],x[j],side)
    ut += U(rij,eps,sig)
  end
  return ut
end

# Compute total potential energy, but do not compute u if rij > 2.
function UT2(x,eps,sig,side)
  ut = 0.
  for i in 1:length(x)-1, j in i+1:length(x)
    rij = r(x[i],x[j],side)
    if rij <= 2.
      ut += U(rij,eps,sig)
    end
  end
  return ut
end

# Generate 100 points
x = [ 10*rand(2) for i in 1:100 ]
# Test if UT1 provides the same result as UT2
eps = 2.0
sig = 0.5
side = 10.
@test UT1(x,eps,sig,side) ≈ UT2(x,eps,sig,side)

# Evaluate code speed
print("UT1: "); @btime UT1($x,$eps,$sig,$side)
print("UT2: "); @btime UT2($x,$eps,$sig,$side)


text:
include("../src/PBCM_IO.jl")

cell_general = cell()
cell_design_op = cylindrical_cell_designer(cell_general)
cost = cost_default()
cell_general, cost = convert_all(cell_general, cost, mult)

using PyPlot

    ###################################################### USD/kWh vs MWH/year  #####################################################################
cell_general.cathode.por = 0.3
cost.cell_costs.anode.AM[1] = 10.57
cost.cell_costs.cathode.AM[1] = 17.0
cell_general.cathode.th       = converter(150.0, mult.pos_th)
cost_calc(cell_general, cost, system="Cell", cost_verbosity = 0)

clf()
data = ["NCA", "NMC622", "LFP"]
AM_rho_data = [4.71, 4.8, 3.6]
rev_sp_cap_data  = [0.180, 0.171, 0.150]
AM_price_data = [26.43, 17.83, 5.07]

for j in range(1, 3, step=1)
    cell_general.cathode.AM_rho = AM_rho_data[j]
    cell_general.cathode.rev_sp_cap = rev_sp_cap_data[j]
    cost.cell_costs.cathode.AM[1] =  AM_price_data[j]
    param_x = []
    param_y = []

    for i in range(10, 120, length=1000)
        # cost.material.pos_AM[1] = i
        cost.general_costs.no_units_mfg = converter(i, mult.units_mfg)
        dollars_per_kWh, mfg_capacity = cost_calc(cell_general, cost, system="Cell", cost_verbosity=0)
        append!(param_y, dollars_per_kWh)
        append!(param_x, mfg_capacity)
    end
    plot(param_x, param_y)
end

xlabel(String("MWh/year"))
ylabel("Dollars per kWh")
title("Cost vs. MFG Capacity")
legend(data)
figure(2)


text:
#!/usr/bin/env julia

##########################################################################
# plot_SNR.jl
#
# Plots SNR curves.
##########################################################################

using CoordinatedPrecoding
using Compat, JLD, ArgParse

##########################################################################
# Load data
s = ArgParseSettings()
@add_arg_table s begin
    "file_names"
        help = "file name with results"
        required = true
        nargs = '+'
end
parsed_args = parse_args(s)

##########################################################################
# Plot parameters
plot_params = @compat Dict(
    "plot_name" => "",

    "objective" => :sum,

    "figure" => Dict(
        :figsize => (8,5),
        :dpi => 125,
        :facecolor => "w",
        :edgecolor => "k",
    ),

    "axes" => Dict(
        :xlabel => "Transmit power [dBm]",
        :ylabel => "Sum rate [bits/s/Hz]",
        :xscale => "linear",
        :yscale => "linear",
    ),

    "legend" => Dict(
        :loc => "best",
        :fontsize => 8,
    ),

    # "confidence_interval_zalpha2" => 1.96,

    "methods" => Dict(
        "Shi2011_WMMSE" => [
            ("logdet_rates", Dict(:color => "b", :linestyle => "-", :label => "WMMSE (logdet)")),
            ("MMSE_rates", Dict(:color => "b", :linestyle => "--",  :label => "WMMSE (MMSE)")),
        ],

        "Gomadam2008_MaxSINR" => [
            ("logdet_rates", Dict(:color => "r", :linestyle => "-", :label => "MaxSINR (logdet)")),
            ("MMSE_rates", Dict(:color => "r", :linestyle => "--",  :label => "MaxSINR (MMSE)")),
        ],

        "Gomadam2008_MinWLI" => [
            ("logdet_rates", Dict(:color => "gray", :linestyle => "-", :label => "MinWLI (logdet)")),
            ("MMSE_rates", Dict(:color => "gray", :linestyle => "--",  :label => "MinWLI (MMSE)")),
        ],

        "Komulainen2013_WMMSE" => [
            ("logdet_rates", Dict(:color => "y", :linestyle => "-", :label => "DiagWMMSE (logdet)")),
            ("MMSE_rates", Dict(:color => "y", :linestyle => "--",  :label => "DiagWMMSE (MMSE)")),
        ],

        "Razaviyayn2013_MinMaxWMMSE" => [
            ("logdet_rates", Dict(:color => "m", :linestyle => "-", :label => "MinMaxWMMSE (logdet)")),
            ("MMSE_rates", Dict(:color => "m", :linestyle => "--",  :label => "MinMaxWMMSE (MMSE)")),
        ],

        "Eigenprecoding" => [
            ("intercell_tdma_logdet_rates", Dict(:color => "c", :linestyle => "-", :label => "TDMA")),
            ("intracell_tdma_logdet_rates", Dict(:color => "c", :linestyle => "-.",  :label => "Intracell TDMA")),
            ("uncoord_logdet_rates", Dict(:color => "k", :linestyle => "-", :label => "Uncoord. transm.")),
        ],
    ),
)

##########################################################################
# Plot it
for file_name in parsed_args["file_names"]
    data = load(file_name)
    processed_results = postprocess(data["raw_results"], data["simulation_params"], plot_params)
    plot(processed_results, data["simulation_params"], plot_params)
end


text:
function basis_rels(b::Array{nf_elem, 1}, c; bd::fmpz = fmpz(10^35), no_b::Int = 250, no_rel::Int = 10000, no_coeff::Int = 4, no_id::fmpz = fmpz(0) )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)
  for i=1:no_rel
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Nemo.sub!(a, a, b[p])
      end
    end
    iszero(a) && continue
    n = norm_div(a, one, no_b)
    if cmpabs(numerator(n), bd) <= 0
      if no_id != 0
        g = gcd(no_id, numerator(n))
        if g==1 || gcd(div(numerator(n), g), g) == 1
          if class_group_add_relation(c, a, n, one)
            a = b[1].parent()
          end
        end
      else
        if class_group_add_relation(c, a, n, one)
          a = b[1].parent()
        end
      end
    end
  end
end

function basis_rels_2(b::Array{nf_elem, 1}, bd::fmpz = fmpz(10^35), no_b::Int = 250, no_rel::Int = 10000, no_coeff::Int = 4, no_id::fmpz = fmpz(0), smooth = 0 )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)
  rels = Array{Hecke.nf_elem}(no_rel)
  i = 1
  l = 0
  while i < no_rel + 1
    l = l + 1
    #println(l)
    if l % 1000 == 0
      println("so far $l tries")
    end
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Hecke.sub!(a, a, b[p])
      end
    end
    iszero(a) && continue
    if no_id != 0
      n = norm_div(a, no_id, no_b)
    else
      n = norm_div(a, one, no_b)
    end
    if smooth != 0
      !issmooth(smooth, numerator(n))[1] && continue
    end
    if cmpabs(numerator(n), bd) <= 0
      if no_id != 0
        g = gcd(no_id, numerator(n))
        if g==1 || gcd(div(numerator(n), g), g) == 1
          rels[i] = deepcopy(a)
          i = i +1
          println(i)
        end
      else
        rels[i] = deepcopy(a)
        i = i + 1
        println(i)
      end
    end
  end
  return rels
end

function basis_rels_3(b::Array{nf_elem, 1}, no_b::Int = 250, no_rel::Int = 10000, no_coeff::Int = 5, no_id::fmpz = fmpz(1), smooth = 0 )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  rels = Dict{fmpz, nf_elem}()
  i = 1
  l = 0
  while i < no_rel + 1
    l = l + 1
    #println(l)
    if l % 1000 == 0
      println("so far $l tries")
    end
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      cf = 1
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Hecke.sub!(a, a, b[p])
      end
    end
    iszero(a) && continue
    n = norm_div(a, no_id, no_b)
    if smooth != 0
      !issmooth(smooth, numerator(n))[1] && continue
    end
    nn = abs(numerator(n))
    if !haskey(rels, nn)
      rels[nn] = deepcopy(a)
      i = i + 1
      println(i)
    end
  end
  return rels
end

function local_norm!(n::fmpz, ap::Array{fq_nmod, 1}, me::Hecke.modular_env)
  nn = UInt(1)
  np = UInt(1)
  for j=1:length(ap)
    ccall((:fq_nmod_norm, :libflint), Nothing, (Ref{fmpz}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), n, ap[j], ap[j].parent)
    nn = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt), nn, UInt(n), me.up, me.upinv)
  end
  ccall((:fmpz_set_ui, :libflint), Nothing, (Ref{fmpz}, UInt), n, nn)
  return n
end

function local_norm(a::nf_elem, me::Hecke.modular_env)
  Fpx = me.Fpx
  aa = Fpx(a)
  ff = Fpx(parent(a).pol)
  np = resultant(aa, ff, !true)
  return fmpz(np)
end

function basis_rels_4(b::Array{nf_elem, 1}, no_b::Int = 250, no_rel::Int = 10000, no_coeff::Int = 5, smooth = 0 )
  a = b[1].parent()
  t = b[1].parent()

  K = b[1].parent

  n = degree(K)

  if Int==Int32
    p = 2^30
  else
    p = 2^60
  end
  sp = PrimesSet(p, 2*p, 1024, 1)
  st = start(sp)
  p, st = next(sp, st)
  lp = [fmpz(p)]
  i = no_b
  while i>0
    p, st = next(sp, st)
    push!(lp, fmpz(p))
    i -= nbits(p[end])
  end

  crt_env = Hecke.crt_env(lp)
  np = Array{fmpz}(length(lp))
  for i=1:length(lp)
    np[i] = fmpz(0)
  end

  lpx = [modular_init(K, p) for p=lp]
  bp = Array{fq_nmod}(length(lpx), n, n)
  for i=1:length(lpx)
    for k=1:n
      ap = Hecke.modular_proj(b[k], lpx[i])
      for l=1:n
        bp[i, k, l] = ap[l]
      end
    end
  end
#  bp = [[deepcopy(Hecke.modular_proj(x, me)) for x = b] for me = lpx]

  tmp = Array{fq_nmod}(length(lpx), n)
  for i=1:length(lpx)
    for j=1:n
      tmp[i,j] = zero(parent(bp[i, 1, 1]))
    end
  end

#  tmp = [[zero(parent(x)) for x=bp[i][1]] for i=1:length(lpx)]

  lc = Array{Int, 1}()
  for i=1:no_coeff
    push!(lc, 0)
  end

  nb = length(b)

  rels = Dict{fmpz, Array{Int, 1}}()
  i = 1
  ll = 0
  sum_nb = 0
  while i < no_rel + 1
    ll = ll + 1
    if ll % 1000 == 0
      println("so far $ll tries, avg nbits of norm ", sum_nb/ll)
    end
    for j=1:no_coeff
      p  = rand(1:nb)
      @inbounds lc[j] = p
    end
#    println("lc: $lc")
    for j=1:length(lpx)
      for k=1:n
        @inbounds ccall((:fq_nmod_set, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}), tmp[j, k], bp[j, lc[1], k])
        for l = 2:no_coeff
#          tmp[j][k] +=  bp[j][lc[l]][k]
          @inbounds add!(tmp[j, k], tmp[j, k], bp[j, lc[l], k])
        end
      end
#      @assert tmp[j] == Hecke.modular_proj(sum(b[lc]), lpx[j])
      local_norm!(np[j], tmp[j,:], lpx[j])
#      println(lpx[j])
#      println(np[j])
#      println(local_norm(sum(b[lc]), lpx[j]))
#      @assert local_norm(sum(b[lc]), lpx[j]) == np[j]
    end

    no = Hecke.crt_signed(np, crt_env)
    sum_nb += nbits(no)
#    println("testing $no")

    if smooth != 0
      !issmooth(smooth, no)[1] && continue
    end
    nn = abs(no)
    if !haskey(rels, nn)
      rels[nn] = deepcopy(lc)
      i = i + 1
      println(i)
    else
      println("again $nn")
    end
  end
  return rels
end

function local_norm!(n::fmpz, ap::nmod_mat, me::Hecke.modular_env)
  nn = UInt(1)
  np = UInt(1)
  for j=1:nrows(ap)
    np = Nemo.getindex_raw(ap, j, 1)
    nn = ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt), nn, np, me.up, me.upinv)
  end
  ccall((:fmpz_set_ui, :libflint), Nothing, (Ref{fmpz}, UInt), n, nn)
  return n
end

function basis_rels_5(b::Array{nf_elem, 1}, no_b::Int = 250, no_rel::Int = 10000, no_coeff::Int = 5, smooth = 0 )
  a = b[1].parent()
  t = b[1].parent()

  K = b[1].parent

  n = degree(K)

  if Int==Int32
    p = 2^30
  else
    p = 2^60
  end
  sp = PrimesSet(p, 2*p, 1024, 1)
  st = start(sp)
  p, st = next(sp, st)
  lp = [fmpz(p)]
  i = no_b
  while i>0
    p, st = next(sp, st)
    push!(lp, fmpz(p))
    i -= nbits(p[end])
  end
  println("using primes ", lp, " for modular norm")

  crt_env = Hecke.crt_env(lp)
  np = Array{fmpz}(length(lp))
  for i=1:length(lp)
    np[i] = fmpz(0)
  end

  lpx = [modular_init(K, p) for p=lp]
  bp = Array{gfp_mat}(length(lpx))
  for i=1:length(lpx)
    bp[i] = zero_matrix(GF(lpx[i].p, cached=false), n, n)
    for k=1:n
      ap = modular_proj(b[k], lpx[i])
      for l=1:n
        bp[i][l, k] = coeff(ap[l], 0)
      end
    end
  end
#  bp = [[deepcopy(modular_proj(x, me)) for x = b] for me = lpx]

  tmp = Array{gfp_mat}(length(lpx))
  lcp = Array{gfp_mat}(length(lpx))
  for i=1:length(lpx)
    tmp[i] = zero_matrix(GF(lpx[i].p, cached=false), n, 1)
    lcp[i] = zero_matrix(GF(lpx[i].p, cached=false), n, 1)
  end

  lc = Array{Int, 1}()
  for i=1:no_coeff
    push!(lc, 0)
  end

  nb = length(b)

  rels = Dict{fmpz, Array{Int, 1}}()
  i = 1
  ll = 0
  sum_nb = 0
  @inbounds while i < no_rel + 1
    ll = ll + 1
    if ll % 1000 == 0
      println("so far $ll tries, avg nbits of norm ", sum_nb/ll)
    end
    for j=1:no_coeff
      lc[j]=0
    end
    for j=1:no_coeff
      p  = rand(1:nb)
      while p in lc
        p = rand(1:nb)
      end
      lc[j] = p
    end
    zero = false;
    for j=1:length(lpx)
      zero!(lcp[j])
      for k=lc
        Nemo.setindex_raw!!(lcp[j], Nemo.getindex_raw(lcp[j], k, 1) + UInt(1), k, 1)
        #should be unlikely - but seems to happen a lot:
        # duplication!
      end
      mul!(tmp[j], bp[j], lcp[j])
      local_norm!(np[j], tmp[j], lpx[j])
      if np[j] == UInt(0) zero = true; break; end;
#      if local_norm(sum(b[lc]), lpx[j]) != np[j]
#        return j, lpx, np[j], tmp[j], lc, bp, lcp
#      end
    end

    if zero continue; end

    no = crt_signed(np, crt_env)
#    @assert no == norm(sum(b[lc]))
    sum_nb += nbits(no)
    if iszero(no) continue; end #cannot (should not) happen. Tested above
#    println("testing $no")

    if smooth != 0
      !issmooth(smooth, no)[1] && continue
    end
    nn = abs(no)
    if !haskey(rels, nn)
#      z = sum(b[lc])
#      nz = norm_div(z, fmpz(1), no_b)
#      @assert norm_div(sum(b[lc]), fmpz(1), no_b) == no
      rels[nn] = deepcopy(lc)
      i = i + 1
      println(i)
    else
      println("again $nn")
    end
  end
  return rels
end

#=

Qx,x = PolynomialRing(FlintQQ, "a")
K, a = MaximalRealSubfield(1024, "a");
@time fb_int = Hecke.int_fb_max_real(1024, 2^20);
h = Hecke.auto_of_maximal_real(K, 3);
b = [K(1), a]
while length(b) < 256 push!(b, h(b[end])); end
fb_int = FactorBase(fmpz[x for x = vcat(fb_int[1], fb_int[2], fb_int[3])]);
@time Hecke.basis_rels_4(b, 600, 10, 5, fb_int)
@time Hecke.basis_rels_5(b, 600, 10, 5, fb_int)

Qx,x = PolynomialRing(FlintQQ, "a")
K, a = MaximalRealSubfield(512, "a");
@time fb_int = Hecke.int_fb_max_real(512, 2^18);
h = Hecke.auto_of_maximal_real(K, 3);
b = [K(1), a]
while length(b) < 128 push!(b, h(b[end])); end
fb_int = FactorBase(fmpz[x for x = vcat(fb_int[1], fb_int[2], fb_int[3])]);
@time basis_rels_4(b, 300, 100, 5, fb_int)
@time basis_rels_5(b, 300, 100, 5, fb_int)

=#

function improve(c::Hecke.ClassGrpCtx)
  H = sub(c.M, 1:nrows(c.M), 1:ncols(c.M))
  Hecke.upper_triangular(H, mod = 17)
  p = setdiff(Set(1:ncols(H)), Set([x.pos[1] for x=H.rows]))
  p = maximum(p)
  b = Hecke.bkz_basis(c.FB.ideals[p]);
#  b = rels_stat(b, ...)
  for x=b
    class_group_add_relation(c, b, n, one)
  end
end

function rels_stat(b::Array{Hecke.nf_elem, 1}; no_b = 250, no_rel::Int = 10000, no_coeff::Int = 4, fixed = 0, smooth=0 )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)

  stat = Dict{Int, Int}()
  if fixed != 0
    stat[-1] = 0
  end
  if smooth != 0
    stat[-2] = 0
  end
  all_g = Array{Any, 1}()
  for i=1:no_rel
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Nemo.sub!(a, a, b[p])
      end
    end
    if iszero(a) continue; end;
    if fixed!=0 && fixed(a) == a
      stat[-1] += 1
    end
    n = norm_div(a, one, no_b)
    k = div(ndigits(numerator(n), 2), 5)
    if haskey(stat, k)
      stat[k] += 1
    else
      stat[k] = 1
    end
    if smooth != 0 && issmooth(smooth, numerator(n))[1]
      stat[-2] += 1
      push!(all_g, a)
      a =   a = b[1].parent()
    end
  end
  return stat, all_g
end

function int_fb_max_real(f::Int, B::Int)
  p = 2
  l1 = []
  lr = []
  lu = []
  sB = sqrt(B)

  while p <=  B
    if p % f == 1 || p % f == f-1
      push!(l1, p)
    elseif f % p == 0
      push!(lr, p)
    elseif p <= sB
      pp = fmpz(p)^2
      while pp <= B
        if pp % f == 1 || pp % f == -1
          push!(lu, p)
          break
        else
          pp *= p
        end
      end
    end
    p = next_prime(p)
  end
  return l1, lr, lu
end

function statistic(st::Dict{Int,T}) where T
  s = T(0)
  s2 = T(0)
  n = 0
  for k in keys(st)
    if k<1
      continue
    end
    t = st[k]
    s += k*t
    s2 += k^2*t
    n += t
  end
  av = s/n
  si = s2/n-av^2

  return av, si
end

function statistic(st::Array{T, 1}) where T
  s = T(0)
  s2 = T(0)
  n = length(st)
  for t in st
    s += t
    s2 += t^2
  end
  av = s/n
  si = s2/n-av^2

  return av, si
end

# If M is upper-triangular with more columns then rows,
# this function returns the non-pivot column indices.
function _find_missing_pivot(M::SMat)
  return setdiff(Set(1:ncols(M)), Set([y.pos[1] for y = M.rows ]))
end

function res_degree_in_max_real(p::Int, n::Int)
  r = remove(n, 2)[2]
  r != 1 && error("Not yet implemented")
  f = Hecke.modord(p, n)
  return (powermod(p, div(f,2), n) == n-1) ? (return div(f, 2)) : (return f)
end


text:
using CellMLModelRepository
using Test

@testset "CellMLModelRepository.jl" begin
    @testset "data" begin include("data.jl") end
end


text:
@testset "Other" begin
    @test LearnBase.Minimizable <: Any
    @test LearnBase.Transformation <: Any
    @test LearnBase.StochasticTransformation <: LearnBase.Transformation

    @test typeof(LearnBase.transform) <: Function
    @test typeof(LearnBase.transform!) <: Function
    @test typeof(LearnBase.learn) <: Function
    @test typeof(LearnBase.learn!) <: Function
    @test typeof(LearnBase.update) <: Function
    @test typeof(LearnBase.update!) <: Function

    @test typeof(LearnBase.grad) <: Function
    @test typeof(LearnBase.grad!) <: Function
    @test typeof(LearnBase.prox) <: Function
    @test typeof(LearnBase.prox!) <: Function

    # IntervalSet
    let s = LearnBase.IntervalSet(-1,1)
        @test typeof(s) == LearnBase.IntervalSet{Int}
        @test typeof(s) <: AbstractSet
        for x in (-1,0,0.5,1,1.0)
            @test x in s
        end
        for x in (-1-1e-10, 1+1e-10, -Inf, Inf, 2, NaN)
            @test !(x in s)
        end
        for i=1:10
            x = rand(s)
            @test typeof(x) == Float64
            @test x in s
        end
        xs = rand(s, 10)
        @test typeof(xs) == Vector{Float64}
        for x in xs
            @test typeof(x) == Float64
            @test x in s
        end
        @test LearnBase.randtype(s) == Float64
        # @show s LearnBase.randtype(s)
    end
    let s = LearnBase.IntervalSet(-1,1.0)
        @test typeof(s) == LearnBase.IntervalSet{Float64}
        @test typeof(s) <: AbstractSet
        @test 1 in s
        # @show s LearnBase.randtype(s)
        @test length(s) == 1
    end

    # IntervalSet{Vector}
    let s = LearnBase.IntervalSet([-1.,0.], [1.,1.])
        @test typeof(s) == LearnBase.IntervalSet{Vector{Float64}}
        @test typeof(s) <: AbstractSet
        @test LearnBase.randtype(s) == Vector{Float64}
        @test typeof(rand(s)) == Vector{Float64}
        @test rand(s) in s
        @test [-1, 0] in s
        @test !([-1.5,0] in s)
        @test !([0,2] in s)
        @test length(s) == 2
    end

    # DiscreteSet
    let s = LearnBase.DiscreteSet([-1,1])
        @test typeof(s) == LearnBase.DiscreteSet{Vector{Int}}
        @test typeof(s) <: AbstractSet
        for x in (-1, 1, -1.0, 1.0)
            @test x in s
        end
        for x in (0, Inf, -Inf, NaN)
            @test !(x in s)
        end
        for i=1:10
            x = rand(s)
            @test typeof(x) == Int
            @test x in s
        end
        xs = rand(s, 10)
        @test typeof(xs) == Vector{Int}
        for x in xs
            @test typeof(x) == Int
            @test x in s
        end
        @test LearnBase.randtype(s) == Int
        @test length(s) == 2
        @test s[1] == -1
    end
    let s = LearnBase.DiscreteSet([-1,1.0])
        @test typeof(s) == LearnBase.DiscreteSet{Vector{Float64}}
        @test typeof(s) <: AbstractSet
        @test typeof(rand(s)) == Float64
        @test typeof(rand(s, 2)) == Vector{Float64}
    end

    # TupleSet
    let s = LearnBase.TupleSet(LearnBase.IntervalSet(0,1), LearnBase.DiscreteSet([0,1]))
        @test typeof(s) == LearnBase.TupleSet{Tuple{LearnBase.IntervalSet{Int}, LearnBase.DiscreteSet{Vector{Int}}}}
        @test typeof(s) <: AbstractSet
        for x in ([0,0], [0.0,0.0], [0.5,1.0])
            @test x in s
        end
        for x in ([0,0.5], [-1,0])
            @test !(x in s)
        end
        @test typeof(rand(s)) == Vector{Float64}
        @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
        @test typeof(rand(s, Tuple)) == Tuple{Float64,Int}
        @test typeof(rand(s, Tuple, 2)) == Vector{Tuple{Float64,Int}}
        @test LearnBase.randtype(s) == Vector{Float64}

        tot = 0
        for (i,x) in enumerate(s)
            @test x == s.sets[i]
            tot += length(x)
        end
        @test length(s) == tot
    end

    # arrays of sets
    let s = [LearnBase.IntervalSet(0,1), LearnBase.DiscreteSet([0,1])]
        @test typeof(s) == Vector{AbstractSet}
        for x in ([0,0], [0.0,0.0], [0.5,1.0])
            @test x in s
        end
        for x in ([0,0.5], [-1,0])
            @test !(x in s)
        end
        @test typeof(rand(s)) == Vector{Float64}
        @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
    end
end


text:
#VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

using Base.Dates, TimeSeries, KernelDensity, Colors

module Watchmaker

using Base.Dates, Requires, TimeSeries, Requests, UnicodePlots, Winston, KernelDensity, Colors

import Base: hist, first, last

export yahoo, fred,
       plot, kplot,
       uplot, uhist,
       colors, palette, default_colors

###### include ##################

include("timearray.jl")
include("fred.jl")
include("yahoo.jl")
include("plot.jl")
include("kplot.jl")
include("uplot.jl")
include("uhist.jl")
include("colors.jl")
include("const.jl")
include("first.jl")
include("last.jl")

end


text:
"""
Multiply a matrix with its own adjoint, obtaining a symmetric/hermitian result where
relevant.

The only exported symbol is [`SELF`](@ref), use eg `SELF' * A`.
"""
module SymmetricProducts

using LinearAlgebra: Diagonal, Hermitian, Symmetric

import Base: show
import LinearAlgebra: adjoint, *

export SELF

struct Self end

"""
A placeholder for the other argument in a binary operation.

Only a few methods are supported, typically for the adjoint `SELF'`, where the
result is special.

# Examples

`SELF'*A` is equivalent to `A'*A`, except for producing a symmetric matrix.
"""
const SELF = Self()

show(io::IO, ::Self) = print(io, "SELF")

struct SelfAdjoint end

show(io::IO, ::SelfAdjoint) = print(io, "SELF'")

adjoint(::Self) = SelfAdjoint()

# matrices (general)

*(A::AbstractMatrix{<: Real}, ::SelfAdjoint) = Symmetric(A*A')

*(::SelfAdjoint, A::AbstractMatrix{<: Real}) = Symmetric(A'*A)

*(A::AbstractMatrix{<: Complex}, ::SelfAdjoint) = Hermitian(A*A')

*(::SelfAdjoint, A::AbstractMatrix{<: Complex}) = Hermitian(A'*A)

# matrices (special)

*(D::Diagonal, ::SelfAdjoint) = Diagonal(abs2.(parent(D)))

*(::SelfAdjoint, D::Diagonal) = D*SELF'

# vectors

*(v::AbstractVector{<: Real}, ::SelfAdjoint) = Symmetric(v*v')

*(v::AbstractVector{<: Complex}, ::SelfAdjoint) = Hermitian(v*v')

*(::SelfAdjoint, v::AbstractVector) = sum(abs2, v)

end # module

max_stars_repo_name: tpapp/SymmetricProducts

max_issues_repo_name: tpapp/SymmetricProducts

max_forks_repo_name: tpapp/SymmetricProducts


text:
"""
solve for the parameters in the similarity function in the rep-FAM model
"""
function solve_rep_param(; d=[1.0, 10.0], p=[.01, .5])
  alpha = log(log(1-p[1]) / log(1-p[2])) / log(abs(d[1] / d[2]))
  phi = - abs(d[1]) ^ alpha / log(1-p[1])
  return (alpha=alpha, phi=phi)
end

#= Test
using Distributions
using RCall
include("repFAM.jl")
d = [7, 10]
p = [.01, .99]
alpha, phi = solve_rep_param(d=d, p=p)
@assert abs(1-exp(-d[2]^alpha/phi) - p[2]) < 1E-6

R"d = 0:32; plot(d, 1-exp(-d^$alpha/$phi), type='o', pch=20); abline(v=$d[2], h=$p[2], lty=2)"
=#

"""
Simulate random bit vectors

For a bit vector of length K, estimate the probability of any two random bit vectors having
at most k ∈ 1:K repeats.
"""
function simulate_prob_differing_bits(bit_vec_length::Int; nsims=100)
  base_bitvec = fill(false, bit_vec_length)
  sim_vec() = [rand() > .5 for j in 1:bit_vec_length]
  num_diffs() = sum(sim_vec() .!= base_bitvec)
  est_prob_num_diffs(k_repeats::Int) = mean(num_diffs() <= k_repeats for i in 1:nsims)

  return [est_prob_num_diffs(k) for k in 1:bit_vec_length]
end

"""
Find the tolerable number of repeated elements in a bit vec of length `bit_vec_length`,
based on a maximum probability of repeats being `thresh_prob`.
"""
function tolerable_num_diffs(bit_vec_length::Int; thresh_prob::Float64=.01, nsims::Int=100)
  sims = simulate_prob_differing_bits(bit_vec_length, nsims=nsims)
  return minimum(findall(sims .> thresh_prob))
end

"""
default parameters psi=(alpha, phi) for repFAM similarity function
"""
function default_psi(bit_vec_length::Int; thresh_probs=(.01, .5), nsims::Int=100)
  tol_num_diff_bits = tolerable_num_diffs(bit_vec_length, thresh_prob=thresh_probs[1], nsims=nsims)
  return solve_rep_param(d=[tol_num_diff_bits, tol_num_diff_bits+1], p=thresh_probs)
end

"""
return default similarity function for repFAM
"""
function gen_similarity_fn(bit_vec_length::Int, thresh_probs=(.01, .5), nsims::Int=100)
  psi = default_psi(bit_vec_length, thresh_probs=thresh_probs, nsims=nsims)
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> exp(-sum(abs.(z1 - z2)) ^ psi.alpha / psi.phi)
end

"""
return default similarity function for repFAM, given probability of repeats
"""
function gen_similarity_fn(; repeats=(1,3), thresh_probs=(.01, .8))
  psi = solve_rep_param(d=repeats, p=thresh_probs)
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> exp(-sum(abs.(z1 - z2)) ^ psi.alpha / psi.phi)
end

"""
Generate a function of threshold, which returns similarity of one if distance
between two vectors < threshold, and zero otherwise.
"""
function sim_fn_abs(threshold::Int)::Function
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> sum(abs.(z1 - z2)) <= threshold ? 1.0 : 0.0
end

function sim_fn_exp_decay(z1::Vector{Bool}, z2::Vector{Bool}, phi::Float64)
  return exp(-sum(abs.(z1 - z2)) / phi)
end

function penalty_exp_decay(z1::Vector{Bool}, z2::Vector{Bool}, phi::Float64)
  return 1 - sim_fn_exp_decay(z1, z2, phi)
end

function log_penalty_exp_decay(z1::Vector{Bool}, z2::Vector{Bool}, phi::Float64)
  return MCMC.log1m(sim_fn_exp_decay(z1, z2, phi))
end

function sim_fn_exp_decay_generator(phi::Float64)::Function
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> begin
    sim_fn_exp_decay(z1, z2, phi)
  end
end

function log_repulsive_fn_generator(phi::Float64)
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> let
    d = sum(abs.(z1 - z2))
    phi * MCMC.log1m(exp(-d))
  end
end

#= Test
using RCall, Distributions
include("repFAM.jl")
@time x = simulate_prob_differing_bits(25, nsims=1000)
R"plot(1:length($x), $x, xlab='k = number of diffs', ylab='probability of at most k differences', type='o')"
R"abline(h=.01)"

tolerable_num_diffs(25, thresh_prob=.01, nsims=1000)
=#

#=
using RCall, Distributions

J = 32
sim_fn = gen_similarity_fn()
sim_fn = sim_fn_abs(2)
z1 = [rand() > .5 ? 1 : 0 for i in 1:J]
z2 = z1 .+ 0

p = []
for j in 1:J
  z2[j] = z2[j] == 1 ? 0 : 1
  append!(p, sim_fn(z1, z2))
end

R"plot(1:$J, $(1 .- p), type='o', xlab='differences', ylab='prior prob')"
=#

include("update_Z_repFAM.jl")