text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# src/solvers.jl
# Handles conversion of the JuMP Model into a format that can be passed
# through the MathProgBase interface to solvers, and ongoing updating of
# that representation if supported by the solver.
#############################################################################

# Analyze a JuMP Model to determine its traits, and thus what solvers can
# be used to solve the problem
type ProblemTraits
    int::Bool  # has integer variables
    lin::Bool  # has only linear objectives and constraints
    qp ::Bool  # has a quadratic objective function
    qc ::Bool  # has a quadratic constraint
    nlp::Bool  # has general nonlinear objective or constraints
    soc::Bool  # has a second-order cone constraint
    sdp::Bool  # has an SDP constraint (or SDP variable bounds)
    sos::Bool  # has an SOS constraint
    conic::Bool  # has an SDP or SOC constraint
end
function ProblemTraits(m::Model; relaxation=false)
    int = !relaxation && any(c-> !(c == :Cont || c == :Fixed), m.colCat)
    qp = !isempty(m.obj.qvars1)
    qc = !isempty(m.quadconstr)
    nlp = m.nlpdata !== nothing
    soc = !isempty(m.socconstr)
    # will need to change this when we add support for arbitrary variable cones
    sdp = !isempty(m.sdpconstr) || !isempty(m.varCones)
    sos = !isempty(m.sosconstr)
    ProblemTraits(int, !(qp|qc|nlp|soc|sdp|sos), qp, qc, nlp, soc, sdp, sos, soc|sdp)
end
function default_solver(traits::ProblemTraits)
    if traits.nlp
        MathProgBase.defaultNLPsolver
    elseif traits.int || traits.sos
        MathProgBase.defaultMIPsolver
    elseif traits.sdp
        MathProgBase.defaultSDPsolver
    elseif traits.conic
        MathProgBase.defaultConicsolver
    elseif traits.qp || traits.qc
        MathProgBase.defaultQPsolver
    else
        MathProgBase.defaultLPsolver
    end
end

function fillConicRedCosts(m::Model)
    bndidx = 0
    numlinconstr = length(m.linconstr)
    vardual = MathProgBase.getvardual(m.internalModel)
    if m.objSense == :Min
        scale!(vardual, -1)
    end
    for i in 1:m.numCols
        lower = false
        upper = false
        lb, ub = m.colLower[i], m.colUpper[i]

        if lb != -Inf && lb != 0.0
            lower = true
            bndidx += 1
        end
        if ub != Inf && ub != 0.0
            upper = true
            bndidx += 1
        end

        if lower && !upper
            m.redCosts[i] = m.conicconstrDuals[numlinconstr + bndidx]
        elseif !lower && upper
            m.redCosts[i] = m.conicconstrDuals[numlinconstr + bndidx]
        elseif lower && upper
            m.redCosts[i] = m.conicconstrDuals[numlinconstr + bndidx]+m.conicconstrDuals[numlinconstr + bndidx-1]
        end
        m.redCosts[i] += vardual[i]
    end
end

function fillConicDuals(m::Model)

    numRows, numCols = length(m.linconstr), m.numCols

    numBndRows = getNumBndRows(m)
    numSOCRows = getNumSOCRows(m)
    m.conicconstrDuals = try
        MathProgBase.getdual(m.internalModel)
    catch
        fill(NaN, numRows+numBndRows+numSOCRows)
    end
    if isfinite(m.conicconstrDuals[1]) # NaN could mean unavailable
        if m.objSense == :Min
            scale!(m.conicconstrDuals, -1)
        end
        m.linconstrDuals = m.conicconstrDuals[1:length(m.linconstr)]
        m.redCosts = zeros(numCols)
        fillConicRedCosts(m)
    end

end

function solve(m::Model; suppress_warnings=false,
                ignore_solve_hook=(m.solvehook===nothing),
                relaxation=false,
                kwargs...)
    # If the user or an extension has provided a solve hook, call
    # that instead of solving the model ourselves
    if !ignore_solve_hook
        return m.solvehook(m; suppress_warnings=suppress_warnings, kwargs...)::Symbol
    end

    isempty(kwargs) || error("Unrecognized keyword arguments: $(join([k[1] for k in kwargs], ", "))")

    # Clear warning counters
    m.getvalue_counter = 0
    m.operator_counter = 0

    # Remember if the solver was initially unset so we can restore
    # it to be unset later
    unset = m.solver == UnsetSolver()

    # Analyze the problems traits to determine what solvers we can use
    traits = ProblemTraits(m, relaxation=relaxation)

    # Build the MathProgBase model from the JuMP model
    build(m, traits=traits, suppress_warnings=suppress_warnings, relaxation=relaxation)

    # If the model is a general nonlinear, use different logic in
    # nlp.jl to solve the problem
    traits.nlp && return solvenlp(m, traits, suppress_warnings=suppress_warnings)

    # Solve the problem
    MathProgBase.optimize!(m.internalModel)
    stat::Symbol = MathProgBase.status(m.internalModel)

    # Extract solution from the solver
    numRows, numCols = length(m.linconstr), m.numCols
    m.objVal = NaN
    m.colVal = fill(NaN, numCols)
    m.linconstrDuals = Array(Float64, 0)

    discrete = !relaxation && (traits.int || traits.sos)
    if stat == :Optimal
        # If we think dual information might be available, try to get it
        # If not, return an array of the correct length
        if !discrete && !traits.conic
            m.redCosts = try
                MathProgBase.getreducedcosts(m.internalModel)[1:numCols]
            catch
                fill(NaN, numCols)
            end

            m.linconstrDuals = try
                MathProgBase.getconstrduals(m.internalModel)[1:numRows]
            catch
                fill(NaN, numRows)
            end
        end
        # conic duals (currently, SOC only)
        if !discrete && traits.soc && !traits.qp && !traits.qc && !traits.sdp
            fillConicDuals(m)
        end
    else
        # Problem was not solved to optimality, attempt to extract useful
        # information anyway
        suppress_warnings || warn("Not solved to optimality, status: $stat")
        # Some solvers provide infeasibility rays (dual) or unbounded
        # rays (primal) for linear problems. Store these as the solution
        # if the exist.
        if traits.lin
            if stat == :Infeasible
                m.linconstrDuals = try
                    infray = MathProgBase.getinfeasibilityray(m.internalModel)
                    @assert length(infray) == numRows
                    infray
                catch
                    suppress_warnings || warn("Infeasibility ray (Farkas proof) not available")
                    fill(NaN, numRows)
                end
            elseif stat == :Unbounded
                m.colVal = try
                    unbdray = MathProgBase.getunboundedray(m.internalModel)
                    @assert length(unbdray) == numCols
                    unbdray
                catch
                    suppress_warnings || warn("Unbounded ray not available")
                    fill(NaN, numCols)
                end
            end
        end
        # conic duals (currently, SOC only)
        if !discrete && traits.soc && !traits.qp && !traits.qc && !traits.sdp
            if stat == :Infeasible
                fillConicDuals(m)
            end
        end
    end

    # If the problem was solved, or if it terminated prematurely, try
    # to extract a solution anyway. This commonly occurs when a time
    # limit or tolerance is set (:UserLimit)
    if !(stat == :Infeasible || stat == :Unbounded)
        try
            objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
            colVal = MathProgBase.getsolution(m.internalModel)[1:numCols]
            # Rescale off-diagonal terms of SDP variables
            if traits.sdp
                offdiagvars = offdiagsdpvars(m)
                colVal[offdiagvars] /= sqrt(2)
            end
            # Don't corrupt the answers if one of the above two calls fails
            m.objVal = objVal
            m.colVal = colVal
        end
    end

    # The MathProgBase interface defines a conic problem to always be
    # a minimization problem, so we need to flip the objective before
    # reporting it to the user
    if traits.conic && m.objSense == :Max
        m.objVal *= -1
    end

    # If the solver was initially not set, we will restore this status
    # and drop the internal MPB model. This is important for the case
    # where the solver used changes between solves because the user
    # has changed the problem class (e.g. LP to MILP)
    if unset
        m.solver = UnsetSolver()
        if traits.int
            m.internalModelLoaded = false
        end
    end

    # don't keep relaxed model in memory
    relaxation && (m.internalModelLoaded = false)

    # Return the solve status
    stat
end

function isquadsoc(m::Model)
    # check if all quadratic constraints are actually conic
    all_conic = true
    for qconstr in m.quadconstr
        q = copy(qconstr.terms)
        if qconstr.sense == :(>=)
            q *= -1
        end
        if !(all(t->t==0, q.aff.coeffs) && q.aff.constant == 0)
            all_conic = false
            break
        end
        n_pos_on_diag = 0
        off_diag_idx  = 0
        neg_diag_idx  = 0
        n = length(q.qvars1)
        nz = 0
        for i in 1:n
            q.qcoeffs[i] == 0 && continue
            nz += 1
            if q.qvars1[i].col == q.qvars2[i].col
                if q.qcoeffs[i] == 1
                    n_pos_on_diag += 1
                elseif q.qcoeffs[i] == -1
                    if !(neg_diag_idx == off_diag_idx == 0)
                        all_conic = false; break
                    end
                    neg_diag_idx = i
                else
                    all_conic = false; break
                end
            else
                if q.qcoeffs[i] == -1
                    if !(neg_diag_idx == off_diag_idx == 0)
                        all_conic = false; break
                    end
                    off_diag_idx = i
                else
                    all_conic = false; break
                end
            end
        end
        if n_pos_on_diag == nz-1 && neg_diag_idx > 0
            # Plain SOC
        elseif n_pos_on_diag == nz-1 && off_diag_idx > 0
            # Rotated SOC
        else
            all_conic = false
        end
    end
    return all_conic && length(m.quadconstr) > 0

end

# Converts the JuMP Model into a MathProgBase model based on the
# traits of the model
function build(m::Model; suppress_warnings=false, relaxation=false, traits=ProblemTraits(m,relaxation=relaxation))
    # Set solver based on the model's traits if it hasn't provided
    if isa(m.solver, UnsetSolver)
        m.solver = default_solver(traits)
    end

    # If the model is nonlinear, use different logic in nlp.jl
    # to build the problem
    traits.nlp && return _buildInternalModel_nlp(m, traits)

    ## Temporary hack for Mosek, see https://github.com/JuliaOpt/Mosek.jl/pull/67
    if contains("$(typeof(m.solver))","MosekSolver") && isquadsoc(m)
        traits.conic = true
    end

    if traits.conic
        # If there are semicontinuous/semi-integer variables, we will have to
        # adjust the b vector below to construct a valid relaxation. This seems
        # like a pretty marginal case, so let's punt on it for now.
        if relaxation && any(x -> (x == :SemiCont || x == :SemiInt), m.colCat)
            error("Relaxations of conic problem with semi-integer/semicontinuous variables are not currently supported.")
        end

        traits.qp && error("JuMP does not support quadratic objectives for conic problems")

        # Obtain a fresh MPB model for the solver
        # If the problem is conic, we rebuild the problem from
        # scratch every time
        m.internalModel = MathProgBase.ConicModel(m.solver)

        # Build up the objective, LHS, RHS and cones from the JuMP Model...
        f, A, b, var_cones, con_cones = conicdata(m)
        # ... and pass to the solver
        MathProgBase.loadproblem!(m.internalModel, f, A, b, con_cones, var_cones)
    else
        # Extract objective coefficients and linear constraint bounds
        f, rowlb, rowub = prepProblemBounds(m)
        # If we already have an MPB model for the solver...
        if m.internalModelLoaded
            # ... and if the solver supports updating bounds/objective
            if applicable(MathProgBase.setvarLB!, m.internalModel, m.colLower) &&
               applicable(MathProgBase.setvarUB!, m.internalModel, m.colUpper) &&
               applicable(MathProgBase.setconstrLB!, m.internalModel, rowlb) &&
               applicable(MathProgBase.setconstrUB!, m.internalModel, rowub) &&
               applicable(MathProgBase.setobj!, m.internalModel, f) &&
               applicable(MathProgBase.setsense!, m.internalModel, m.objSense)
                MathProgBase.setvarLB!(m.internalModel, copy(m.colLower))
                MathProgBase.setvarUB!(m.internalModel, copy(m.colUpper))
                MathProgBase.setconstrLB!(m.internalModel, rowlb)
                MathProgBase.setconstrUB!(m.internalModel, rowub)
                MathProgBase.setobj!(m.internalModel, f)
                MathProgBase.setsense!(m.internalModel, m.objSense)
            else
                # The solver doesn't support changing bounds/objective
                # We need to build the model from scratch
                if !suppress_warnings
                    Base.warn_once("Solver does not appear to support hot-starts. Model will be built from scratch.")
                end
                m.internalModelLoaded = false
            end
        end
        # If we don't already have a MPB model
        if !m.internalModelLoaded
            # Obtain a fresh MPB model for the solver
            m.internalModel = MathProgBase.LinearQuadraticModel(m.solver)
            # Construct a LHS matrix from the linear constraints
            A = prepConstrMatrix(m)

            # Load the problem data into the model...
            collb = copy(m.colLower)
            colub = copy(m.colUpper)
            if relaxation
                for i in 1:m.numCols
                    if m.colCat[i] in (:SemiCont,:SemiInt)
                        collb[i] = min(0.0, collb[i])
                        colub[i] = max(0.0, colub[i])
                    end
                end
            end
            MathProgBase.loadproblem!(m.internalModel, A, collb, colub, f, rowlb, rowub, m.objSense)
            # ... and add quadratic and SOS constraints separately
            addQuadratics(m)
            if !relaxation
                addSOS(m)
            end
        end
        # Update solver callbacks, if any
        if !relaxation
            registercallbacks(m)
        end
    end

    # Update the type of each variable
    if applicable(MathProgBase.setvartype!, m.internalModel, Symbol[])
        if relaxation
            MathProgBase.setvartype!(m.internalModel, fill(:Cont, m.numCols))
        else
            colCats = vartypes_without_fixed(m)
            MathProgBase.setvartype!(m.internalModel, colCats)
        end
    elseif traits.int
        # Solver that do not implement anything other than continuous
        # variables do not need to implement this method, so throw an
        # error if the model has anything but continuous
        error("Solver does not support discrete variables")
    end

    # Provide a primal solution to the solve, if the problem is integer
    # and the user has provided one.
    # TODO: change this so you can warm start continuous problems?
    if !relaxation && traits.int && !all(isnan(m.colVal))
        if applicable(MathProgBase.setwarmstart!, m.internalModel, m.colVal)
            MathProgBase.setwarmstart!(m.internalModel, m.colVal)
        else
            suppress_warnings || Base.warn_once("Solver does not appear to support providing initial feasible solutions.")
        end
    end

    # Record that we have a MPB model constructed
    m.internalModelLoaded = true
    nothing
end

# Add the quadratic part of the objective and all quadratic constraints
# to the internal MPB model
function addQuadratics(m::Model)
    # The objective function is always a quadratic expression, but
    # may have no quadratic terms (i.e. be just affine)
    if length(m.obj.qvars1) != 0
        # Check that no coefficients are NaN/Inf
        assert_isfinite(m.obj)
        # Check that quadratic term variables belong to this model
        # Affine portion is checked in prepProblemBounds
        if !(verify_ownership(m, m.obj.qvars1) &&
                verify_ownership(m, m.obj.qvars2))
            error("Variable not owned by model present in objective")
        end
        # Check for solver support for quadratic objectives happens in MPB
        MathProgBase.setquadobjterms!(m.internalModel,
            Cint[v.col for v in m.obj.qvars1],
            Cint[v.col for v in m.obj.qvars2], m.obj.qcoeffs)
    end

    # Add quadratic constraint to solver
    const sensemap = Dict(:(<=) => '<', :(>=) => '>', :(==) => '=')
    for k in 1:length(m.quadconstr)
        qconstr = m.quadconstr[k]::QuadConstraint
        if !haskey(sensemap, qconstr.sense)
            error("Invalid quadratic constraint sense $(qconstr.sense)")
        end
        s = sensemap[qconstr.sense]

        terms::QuadExpr = qconstr.terms
        # Check that no coefficients are NaN/Inf
        assert_isfinite(terms)
        # Check that quadratic and affine term variables belong to this model
        if !(verify_ownership(m, terms.qvars1) &&
                verify_ownership(m, terms.qvars2) &&
                verify_ownership(m, terms.aff.vars))
            error("Variable not owned by model present in quadratic constraint")
        end
        # Extract indices for MPB, and add the constraint (if we can)
        affidx  = Cint[v.col for v in terms.aff.vars]
        var1idx = Cint[v.col for v in terms.qvars1]
        var2idx = Cint[v.col for v in terms.qvars2]
        if applicable(MathProgBase.addquadconstr!, m.internalModel, affidx, terms.aff.coeffs, var1idx, var2idx, terms.qcoeffs, s, -terms.aff.constant)
            MathProgBase.addquadconstr!(m.internalModel,
                affidx, terms.aff.coeffs,           # aᵀx +
                var1idx, var2idx, terms.qcoeffs,    # xᵀQx
                s, -terms.aff.constant)             # ≤/≥ b
        else
            error("Solver does not support quadratic constraints")
        end
    end
    nothing
end

function addSOS(m::Model)
    for i in 1:length(m.sosconstr)
        sos = m.sosconstr[i]
        indices = Int[v.col for v in sos.terms]
        if sos.sostype == :SOS1
            if applicable(MathProgBase.addsos1!, m.internalModel, indices, sos.weights)
                MathProgBase.addsos1!(m.internalModel, indices, sos.weights)
            else
                error("Solver does not support SOS constraints")
            end
        elseif sos.sostype == :SOS2
            if applicable(MathProgBase.addsos2!, m.internalModel, indices, sos.weights)
                MathProgBase.addsos2!(m.internalModel, indices, sos.weights)
            else
                error("Solver does not support SOS constraints")
            end
        end
    end
end

# Returns coefficients for the affine part of the objective and the
# affine constraint lower and upper bounds, all as dense vectors
function prepProblemBounds(m::Model)

    # Create dense objective vector
    objaff::AffExpr = m.obj.aff
    # Check that no coefficients are NaN/Inf
    assert_isfinite(objaff)
    if !verify_ownership(m, objaff.vars)
        error("Variable not owned by model present in objective")
    end
    f = zeros(m.numCols)
    @inbounds for ind in 1:length(objaff.vars)
        f[objaff.vars[ind].col] += objaff.coeffs[ind]
    end

    # Create dense affine constraint bound vectors
    linconstr = m.linconstr::Vector{LinearConstraint}
    numRows = length(linconstr)
    # -Inf means no lower bound, +Inf means no upper bound
    rowlb = fill(-Inf, numRows)
    rowub = fill(+Inf, numRows)
    @inbounds for ind in 1:numRows
        rowlb[ind] = linconstr[ind].lb
        rowub[ind] = linconstr[ind].ub
    end

    return f, rowlb, rowub
end

# Convert all the affine constraints into a sparse column-wise
# matrix of coefficients.
function prepConstrMatrix(m::Model)

    linconstr = m.linconstr::Vector{LinearConstraint}
    numRows = length(linconstr)
    # Calculate the maximum number of nonzeros
    # The actual number may be less because of cancelling or
    # zero-coefficient terms
    nnz = 0
    for c in 1:numRows
        nnz += length(linconstr[c].terms.coeffs)
    end
    # Non-zero row indices
    I = Array(Int,nnz)
    # Non-zero column indices
    J = Array(Int,nnz)
    # Non-zero values
    V = Array(Float64,nnz)

    # Fill it up!
    # Number of nonzeros seen so far
    nnz = 0
    for c in 1:numRows
        # Check that no coefficients are NaN/Inf
        assert_isfinite(linconstr[c].terms)
        coeffs = linconstr[c].terms.coeffs
        vars   = linconstr[c].terms.vars
        # Check that variables belong to this model
        if !verify_ownership(m, vars)
            error("Variable not owned by model present in a constraint")
        end
        # Record all (i,j,v) triplets
        @inbounds for ind in 1:length(coeffs)
            nnz += 1
            I[nnz] = c
            J[nnz] = vars[ind].col
            V[nnz] = coeffs[ind]
        end
    end

    # sparse() handles merging duplicate terms and removing zeros
    A = sparse(I,J,V,numRows,m.numCols)
end

function vartypes_without_fixed(m::Model)
    colCats = copy(m.colCat)
    for i in 1:length(colCats)
        if colCats[i] == :Fixed
            @assert m.colLower[i] == m.colUpper[i]
            colCats[i] = :Cont
        end
    end
    return colCats
end

function collect_expr!(m, tmprow, terms::AffExpr)
    empty!(tmprow)
    assert_isfinite(terms)
    coeffs = terms.coeffs
    vars = terms.vars
    # collect duplicates
    for ind in 1:length(coeffs)
        if !is(vars[ind].m, m)
            error("Variable not owned by model present in constraints")
        end
        addelt!(tmprow,vars[ind].col, coeffs[ind])
    end
    tmprow
end

# Returns a boolean vector indicating if variable in the model
# is an off-diagonal element of an SDP matrix.
# This is needed because we have to rescale coefficients that
# touch these variables.
function offdiagsdpvars(m::Model)
    offdiagvars = falses(m.numCols)
    for (name,idx) in m.varCones
        if name == :SDP
            conelen = length(idx)
            n = round(Int,sqrt(1/4+2*conelen)-1/2)
            @assert n*(n+1)/2 == conelen
            r = 1
            for i in 1:n
                for j in i:n
                    if i != j
                        offdiagvars[idx[r]] = true
                    end
                    r += 1
                end
            end
        end
    end
    return offdiagvars
end

function conicdata(m::Model)
    var_cones = Any[cone for cone in m.varCones]
    con_cones = Any[]
    nnz = 0

    # find starting column indices for sdp matrices
    numSDPRows = 0
    numSymRows = 0
    for c in m.sdpconstr
        n = size(c.terms,1)
        @assert n == size(c.terms,2)
        @assert ndims(c.terms) == 2
        numSDPRows += convert(Int, n*(n+1)/2)
        for i in 1:n, j in i:n
            nnz += length(c.terms[i,j].coeffs)
        end
        if !Compat.issymmetric(c.terms)
            # symmetry constraints
            numSymRows += convert(Int, n*(n-1)/2)
        end
    end

    soc_cones  = Any[]
    rsoc_cones = Any[]
    numQuadRows = 0
    for qconstr in m.quadconstr
        q = copy(qconstr.terms)
        if qconstr.sense == :(>=)
            q *= -1
        end
        if !(all(t->t==0, q.aff.coeffs) && q.aff.constant == 0)
            error("Quadratic constraint $qconstr must be in second-order cone form")
        end
        n_pos_on_diag = 0
        off_diag_idx  = 0
        neg_diag_idx  = 0
        n = length(q.qvars1)
        nz = 0
        for i in 1:n
            q.qcoeffs[i] == 0 && continue
            nz += 1
            if q.qvars1[i].col == q.qvars2[i].col
                if q.qcoeffs[i] == 1
                    n_pos_on_diag += 1
                elseif q.qcoeffs[i] == -1
                    neg_diag_idx == off_diag_idx == 0 || error("Invalid SOC constraint $qconstr")
                    neg_diag_idx = i
                end
            else
                if q.qcoeffs[i] == -1
                    neg_diag_idx == off_diag_idx == 0 || error("Invalid rotated SOC constraint $qconstr")
                    off_diag_idx = i
                    nz += 1
                end
            end
        end
        cone = Array(Int, nz)
        if n_pos_on_diag == nz-1 && neg_diag_idx > 0
            cone[1] = q.qvars1[neg_diag_idx].col
            r = 1
            for i in 1:n
                (q.qcoeffs[i] == 0 || i == neg_diag_idx) && continue
                r += 1
                cone[r] = q.qvars1[i].col
            end
            push!(soc_cones, cone)
        elseif n_pos_on_diag == nz-2 && off_diag_idx > 0
            cone[1] = q.qvars1[off_diag_idx].col
            cone[2] = q.qvars2[off_diag_idx].col
            r = 2
            for i in 1:n
                (q.qcoeffs[i] == 0 || i == off_diag_idx) && continue
                r += 1
                cone[r] = q.qvars1[i].col
            end
            push!(rsoc_cones, cone)
        else
            error("Quadratic constraint $qconstr is not conic representable")
        end
        numQuadRows += length(cone)
    end

    linconstr = m.linconstr::Vector{LinearConstraint}
    numLinRows = length(linconstr)
    numBounds = 0
    nonNeg  = Int[]
    nonPos  = Int[]
    free    = Int[]
    zeroVar = Int[]
    for i in 1:m.numCols
        seen = false
        lb, ub = m.colLower[i], m.colUpper[i]
        for (_,cone) in m.varCones
            if i in cone
                seen = true
                @assert lb == -Inf && ub == Inf
                break
            end
        end

        if !seen
            if lb != -Inf && lb != 0
                numBounds += 1
            end
            if ub != Inf && ub != 0
                numBounds += 1
            end
            if lb == 0 && ub == 0
                push!(zeroVar, i)
            elseif lb == 0
                push!(nonNeg, i)
            elseif ub == 0
                push!(nonPos, i)
            else
                push!(free, i)
            end
        end
    end

    if !isempty(zeroVar)
        push!(var_cones, (:Zero,zeroVar))
    end
    if !isempty(nonNeg)
        push!(var_cones, (:NonNeg,nonNeg))
    end
    if !isempty(nonPos)
        push!(var_cones, (:NonPos,nonPos))
    end
    if !isempty(free)
        push!(var_cones, (:Free,free))
    end

    nnz += numBounds
    for c in 1:numLinRows
        nnz += length(linconstr[c].terms.coeffs)
    end

    numSOCRows = 0
    numNormRows = 0
    for con in m.socconstr
        numNormRows += 1
        numSOCRows += length(con.normexpr.norm.terms) + 1
    end
    numRows = numLinRows + numBounds + numQuadRows + numSOCRows + numSDPRows + numSymRows

    # should maintain the order of constraints in the above form
    # throughout the code c is the conic constraint index
    # TODO: only added linear+bound+soc support, extend to all
    constr_dual_map = Array(Vector{Int}, numLinRows + numBounds + numNormRows)

    b = Array(Float64, numRows)

    I = Int[]
    J = Int[]
    V = Float64[]
    sizehint!(I, nnz)
    sizehint!(J, nnz)
    sizehint!(V, nnz)

    # Fill it up
    nnz = 0
    tmprow = IndexedVector(Float64,m.numCols)
    tmpelts = tmprow.elts
    tmpnzidx = tmprow.nzidx
    nonneg_rows = Int[]
    nonpos_rows = Int[]
    eq_rows     = Int[]
    for c in 1:numLinRows
        if linconstr[c].lb == -Inf
            b[c] = linconstr[c].ub
            push!(nonneg_rows, c)
        elseif linconstr[c].ub == Inf
            b[c] = linconstr[c].lb
            push!(nonpos_rows, c)
        elseif linconstr[c].lb == linconstr[c].ub
            b[c] = linconstr[c].lb
            push!(eq_rows, c)
        else
            error("We currently do not support ranged constraints with conic solvers")
        end

        assert_isfinite(linconstr[c].terms)
        coeffs = linconstr[c].terms.coeffs
        vars = linconstr[c].terms.vars
        # collect duplicates
        for ind in 1:length(coeffs)
            if !is(vars[ind].m, m)
                error("Variable not owned by model present in constraints")
            end
            addelt!(tmprow,vars[ind].col, coeffs[ind])
        end
        nnz = tmprow.nnz
        append!(I, fill(c, nnz))
        indices = tmpnzidx[1:nnz]
        append!(J, indices)
        append!(V, tmpelts[indices])
        empty!(tmprow)
        constr_dual_map[c] = vec(collect(c))
    end

    c = numLinRows
    bndidx = 0
    for idx in 1:m.numCols
        lb = m.colLower[idx]
        if lb != -Inf && lb != 0
            bndidx += 1
            nnz += 1
            c   += 1
            push!(I, c)
            push!(J, idx)
            push!(V, 1.0)
            b[c] = lb
            push!(nonpos_rows, c)
            constr_dual_map[numLinRows + bndidx] = vec(collect(c))
        end
        ub = m.colUpper[idx]
        if ub != Inf && ub != 0
            bndidx += 1
            c   += 1
            push!(I, c)
            push!(J, idx)
            push!(V, 1.0)
            b[c] = ub
            push!(nonneg_rows, c)
            constr_dual_map[numLinRows + bndidx] = vec(collect(c))
        end
    end

    if !isempty(nonneg_rows)
        push!(con_cones, (:NonNeg,nonneg_rows))
    end
    if !isempty(nonpos_rows)
        push!(con_cones, (:NonPos,nonpos_rows))
    end
    if !isempty(eq_rows)
        push!(con_cones, (:Zero,eq_rows))
    end
    @assert c == numLinRows + numBounds

    for cone in soc_cones
        n = length(cone)
        rng = (c+1):(c+n)
        append!(I, rng)
        append!(J, copy(cone))
        append!(V, [-1.0; ones(n-1)])
        push!(con_cones, (:SOC,rng))
        b[rng] = 0
        c += n
    end
    for cone in rsoc_cones
        n = length(cone)
        rng = (c+1):(c+n)
        append!(I, rng)
        append!(J, copy(cone))
        append!(V, [-1/sqrt(2); -1/sqrt(2); ones(n-2)])
        push!(con_cones, (:SOCRotated,rng))
        b[rng] = 0
        c += n
    end
    @assert c == numLinRows + numBounds + numQuadRows

    tmpelts = tmprow.elts
    tmpnzidx = tmprow.nzidx
    socidx = 0
    for con in m.socconstr
        socidx += 1
        expr = con.normexpr
        c += 1
        soc_start = c
        collect_expr!(m, tmprow, expr.aff)
        nnz = tmprow.nnz
        indices = tmpnzidx[1:nnz]
        append!(I, fill(c, nnz))
        append!(J, indices)
        append!(V, tmpelts[indices])
        b[c] = -expr.aff.constant
        for term in expr.norm.terms
            c += 1
            collect_expr!(m, tmprow, term)
            nnz = tmprow.nnz
            indices = tmpnzidx[1:nnz]
            append!(I, fill(c, nnz))
            append!(J, indices)
            append!(V, -expr.coeff*tmpelts[indices])
            b[c] = expr.coeff*term.constant
        end
        push!(con_cones, (:SOC, soc_start:c))
        constr_dual_map[numLinRows + numBounds + socidx] = collect(soc_start:c)
    end
    @assert c == numLinRows + numBounds + numQuadRows + numSOCRows

    numDroppedSym = 0
    for con in m.sdpconstr
        sdp_start = c + 1
        n = size(con.terms,1)
        for i in 1:n, j in i:n
            c += 1
            terms::AffExpr = con.terms[i,j]
            collect_expr!(m, tmprow, terms)
            nnz = tmprow.nnz
            indices = tmpnzidx[1:nnz]
            append!(I, fill(c, nnz))
            append!(J, indices)
            # scale to svec form
            scale = (i == j) ? 1.0 : sqrt(2)
            append!(V, -scale*tmpelts[indices])
            b[c] = scale*terms.constant
        end
        push!(con_cones, (:SDP, sdp_start:c))
        if !Compat.issymmetric(con.terms)
            sym_start = c + 1
            # add linear symmetry constraints
            for i in 1:n, j in 1:(i-1)
                collect_expr!(m, tmprow, con.terms[i,j] - con.terms[j,i])
                nnz = tmprow.nnz
                # if the symmetry-enforcing row is empty, just drop it
                if nnz == 0
                    numDroppedSym += 1
                    continue
                end
                c += 1
                indices = tmpnzidx[1:nnz]
                append!(I, fill(c, nnz))
                append!(J, indices)
                append!(V, tmpelts[indices])
                b[c] = 0
            end
            push!(con_cones, (:Zero, sym_start:c))
        end
    end
    numRows -= numDroppedSym
    resize!(b, numRows)
    @assert c == numRows

    m.constrDualMap = constr_dual_map

    # Use only the objective coefficients from prepProblemBounds
    f,_,_ = prepProblemBounds(m)

    # The conic MPB interface defines conic problems as
    # always being minimization problems, so flip if needed
    m.objSense == :Max && scale!(f, -1.0)

    # Objective coefficients and columns of A matrix are
    # rescaled for SDP variables
    offdiagvars = offdiagsdpvars(m)
    f[offdiagvars] /= sqrt(2)
    for k in 1:length(J)
        if offdiagvars[J[k]]
            V[k] /= sqrt(2)
        end
    end

    A = sparse(I, J, V, numRows, m.numCols)
    # @show full(A), b
    # @show var_cones, con_cones

    # TODO: uncomment these lines when they work with Mosek
    # supported = MathProgBase.supportedcones(m.internalModel)
    # @assert (:NonNeg in supported) && (:NonPos in supported) && (:Free in supported) && (:SDP in supported)
    f, A, b, var_cones, con_cones
end

constraintbounds(m::Model) = constraintbounds(m, ProblemTraits(m))

function constraintbounds(m::Model,traits::ProblemTraits)

    if traits.conic
        error("Not implemented for conic problems")
    elseif traits.sos
        error("Not implemented for SOS constraints")
    end

    linobj, linrowlb, linrowub = prepProblemBounds(m)

    quadrowlb = Float64[]
    quadrowub = Float64[]
    for c::QuadConstraint in m.quadconstr
        if c.sense == :(<=)
            push!(quadrowlb, -Inf)
            push!(quadrowub, 0.0)
        elseif c.sense == :(>=)
            push!(quadrowlb, 0.0)
            push!(quadrowub, Inf)
        elseif c.sense == :(==)
            push!(quadrowlb, 0.0)
            push!(quadrowub, 0.0)
        else
            error("Unrecognized quadratic constraint sense $(c.sense)")
        end
    end

    nlrowlb = Float64[]
    nlrowub = Float64[]

    if traits.nlp
        nldata::NLPData = m.nlpdata
        for c in nldata.nlconstr
            push!(nlrowlb, c.lb)
            push!(nlrowub, c.ub)
        end
    end

    lb = [linrowlb;quadrowlb;nlrowlb]
    ub = [linrowub;quadrowub;nlrowub]

    return lb, ub

end

# returns (unsorted) column indices and coefficient terms for merged vector
# assume that v is zero'd
function merge_duplicates{CoefType,IntType<:Integer}(::Type{IntType},aff::GenericAffExpr{CoefType,Variable}, v::IndexedVector{CoefType}, m::Model)
    resize!(v, m.numCols)
    for ind in 1:length(aff.coeffs)
        var = aff.vars[ind]
        is(var.m, m) || error("Variable does not belong to this model")
        addelt!(v, aff.vars[ind].col, aff.coeffs[ind])
    end
    indices = Array(IntType,v.nnz)
    coeffs = Array(CoefType,v.nnz)
    for i in 1:v.nnz
        idx = v.nzidx[i]
        indices[i] = idx
        coeffs[i] = v.elts[idx]
    end
    empty!(v)

    return indices, coeffs

end


text:
"""
    initialize!(rbm, [data]; ϵ = 1e-6)

Initializes the RBM and returns it.
If provided, matches average visible unit activities from `data`.

    initialize!(layer, [data]; ϵ = 1e-6)

Initializes a layer and returns it.
If provided, matches average unit activities from `data`.
"""
function initialize! end

function initialize!(rbm::RBM; ϵ::Real = 1e-6)
    initialize!(visible(rbm))
    initialize!(hidden(rbm))
    initialize_w!(rbm; ϵ)
    zerosum!(rbm)
    return rbm
end

function initialize!(rbm::RBM, data::AbstractArray; ϵ::Real = 1e-6, wts = nothing)
    @assert 0 < ϵ < 1/2
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    initialize!(visible(rbm), data; ϵ, wts)
    initialize!(hidden(rbm))
    initialize_w!(rbm, data; ϵ, wts)
    zerosum!(rbm)
    return rbm
end

function initialize!(layer::Binary, data::AbstractArray; ϵ::Real=1e-6, wts=nothing)
    @assert size(layer) == size(data)[1:ndims(layer)]
    @assert 0 < ϵ < 1/2
    μ = batchmean(layer, data; wts)
    μϵ = clamp.(μ, ϵ, 1 - ϵ)
    @. layer.θ = -log(1/μϵ - 1)
    return layer
end

function initialize!(layer::Spin, data::AbstractArray; ϵ::Real=1e-6, wts=nothing)
    @assert size(layer) == size(data)[1:ndims(layer)]
    @assert 0 < ϵ < 1/2
    μ = batchmean(layer, data; wts)
    μϵ = clamp.(μ, ϵ - 1, 1 - ϵ)
    layer.θ .= atanh.(μϵ)
    return layer
end

function initialize!(layer::Potts, data::AbstractArray; ϵ::Real=1e-6, wts=nothing)
    @assert size(layer) == size(data)[1:ndims(layer)]
    @assert 0 < ϵ < 1/2
    μ = batchmean(layer, data; wts)
    μϵ = clamp.(μ, ϵ, 1 - ϵ)
    layer.θ .= log.(μϵ)
    return layer # does not do zerosum!
end

function initialize!(layer::Gaussian, data::AbstractArray; ϵ::Real=1e-6, wts=nothing)
    @assert size(layer) == size(data)[1:ndims(layer)]
    @assert 0 < ϵ < 1/2
    μ = batchmean(layer, data; wts)
    ν = batchmean(layer, (data .- μ).^2; wts)
    layer.γ .= inv.(ν .+ ϵ)
    layer.θ .= μ .* layer.γ
    return layer
end

function initialize!(layer::Union{Potts, Binary, Spin})
    layer.θ .= 0
    return layer
end

function initialize!(layer::Union{Gaussian, ReLU})
    layer.θ .= 0
    layer.γ .= 1
    return layer
end

function initialize!(layer::dReLU)
    layer.θp .= layer.θn .= 0
    layer.γp .= layer.γn .= 1
    return layer
end

function initialize!(layer::pReLU)
    layer.θ .= layer.Δ .= layer.η .= 0
    layer.γ .= 1
    return layer
end

function initialize!(layer::xReLU)
    layer.θ .= layer.Δ .= layer.ξ .= 0
    layer.γ .= 1
    return layer
end

"""
    initialize_w!(rbm, data; λ = 0.1)

Initializes `rbm.w` such that typical inputs to hidden units are λ.
"""
function initialize_w!(
    rbm::RBM, data::AbstractArray;
    λ::Real = 0.1, ϵ::Real = 1e-6, wts::Union{Nothing,AbstractVector} = nothing
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    if isnothing(wts)
        d = dot(data, data / size(data, ndims(data)))
    else
        @assert length(wts) == size(data)[end]
        x = reshape(data, length(visible(rbm)), size(data)[end])
        d = dot(x * Diagonal(wts), x / sum(wts))
    end
    randn!(weights(rbm))
    weights(rbm) .*= λ / √(d + ϵ)
    return rbm # does not impose zerosum
end

function initialize_w!(rbm::RBM; λ::Real = 0.1, ϵ::Real = 1e-6)
    d = sum(var_from_inputs(visible(rbm)) .+ mean_from_inputs(visible(rbm)).^2)
    randn!(weights(rbm))
    weights(rbm) .*= λ / √(d + ϵ)
    return rbm
end

ext: jl


text:
# Import Turing and Distributions.
using Turing, Distributions
# GLM for comparison
using GLM
# Import RDatasets.
using RDatasets
# Import MCMCChains, Plots, and StatPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots
# Functionality for splitting and normalizing the data.
using MLDataUtils: shuffleobs, splitobs, rescale!
# Functionality for evaluating the model predictions.
using Distances
# set random seeds
using Random
Random.seed!(12212020)

# import RDatasets
df = RDatasets.dataset("datasets","ToothGrowth");
# Dummy code the Supp column
df.Supp_num = convert.(Float64, df.Supp .== "VC")

# split train-test
train, test = splitobs(shuffleobs(df), at=.7)

# put DV and IV into separate vars
function split_iv_dv(df)
    y = Array(df.Len)
    X = Array(df[:,["Supp_num","Dose"]])
    return y, X
end
y_train, X_train = split_iv_dv(train)
y_test, X_test = split_iv_dv(test)

# OLS for comparison
ols = lm(@formula(Len ~ Supp * Dose), train)

# define the generative function
@model function linear_regression(y, X)
    # construct interaction term
    X = [X X[:,1] .* X[:,2]]
    # prior for intercept
    b0 ~ Normal(0,3)
    # prior for betas, mvnormal with 0 means
    bs ~ MvNormal(size(X,2), 3)
    # error variance prior
    e ~ InverseGamma(2,3)

    ## condition on obs
    mu = b0 .+ X * bs
    y ~ MvNormal(mu, e)
end

# prior check, plot sample from prior
prior_sample = sample(
    linear_regression(y_train, X_train),
    Prior(),
    1000
)
plot(prior_sample)
describe(prior_sample)

# NUTS sampling posterior
ch = sample(
    linear_regression(y_train, X_train),
    NUTS(.65),
    1000
)

# visualize the samples
plot(ch)
# summarize the samples
describe(ch)

## check against hold-out data

max_stars_repo_name: zidongzhao/julia-learning

max_issues_repo_name: zidongzhao/julia-learning

max_forks_repo_name: zidongzhao/julia-learning


text:
using ITensors
using Printf

function ITensors.op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
  h =
    1 / 2 * op("S+", s1) * op("S-", s2) +
    1 / 2 * op("S-", s1) * op("S+", s2) +
    op("Sz", s1) * op("Sz", s2)
  return exp(τ * h)
end

"""
Given a Vector of numbers, returns
the average and the standard error
(= the width of distribution of the numbers)
"""
function avg_err(v::Vector)
  N = length(v)
  avg = v[1] / N
  avg2 = v[1]^2 / N
  for j in 2:N
    avg += v[j] / N
    avg2 += v[j]^2 / N
  end
  return avg, √((avg2 - avg^2) / N)
end

function main(; N=10, cutoff=1E-8, δτ=0.1, beta=2.0, NMETTS=3000, Nwarm=10)

  # Make an array of 'site' indices
  s = siteinds("S=1/2", N)

  # Make gates (1,2),(2,3),(3,4),...
  gates = ops([("expτSS", (n, n + 1), (τ=-δτ / 2,)) for n in 1:(N - 1)], s)
  # Include gates in reverse order to complete Trotter formula
  append!(gates, reverse(gates))

  # Make y-rotation gates to use in METTS collapses
  Ry_gates = ops([("Ry", n, (θ=π / 2,)) for n in 1:N], s)

  # Arbitrary initial state
  psi = randomMPS(s)

  # Make H for measuring the energy
  terms = OpSum()
  for j in 1:(N - 1)
    terms += 1 / 2, "S+", j, "S-", j + 1
    terms += 1 / 2, "S-", j, "S+", j + 1
    terms += "Sz", j, "Sz", j + 1
  end
  H = MPO(terms, s)

  # Make τ_range and check δτ is commensurate
  τ_range = δτ:δτ:(beta / 2)
  if norm(length(τ_range) * δτ - beta / 2) > 1E-10
    error("Time step δτ=$δτ not commensurate with beta/2=$(beta/2)")
  end

  energies = Float64[]

  for step in 1:(Nwarm + NMETTS)
    if step <= Nwarm
      println("Making warmup METTS number $step")
    else
      println("Making METTS number $(step-Nwarm)")
    end

    # Do the time evolution by applying the gates
    for τ in τ_range
      psi = apply(gates, psi; cutoff)
      normalize!(psi)
    end

    # Measure properties after >= Nwarm
    # METTS have been made
    if step > Nwarm
      energy = inner(psi', H, psi)
      push!(energies, energy)
      @printf("  Energy of METTS %d = %.4f\n", step - Nwarm, energy)
      a_E, err_E = avg_err(energies)
      @printf(
        "  Estimated Energy = %.4f +- %.4f  [%.4f,%.4f]\n",
        a_E,
        err_E,
        a_E - err_E,
        a_E + err_E
      )
    end

    # Measure in X or Z basis on alternating steps
    if step % 2 == 1
      psi = apply(Ry_gates, psi)
      samp = sample!(psi)
      new_state = [samp[j] == 1 ? "X+" : "X-" for j in 1:N]
    else
      samp = sample!(psi)
      new_state = [samp[j] == 1 ? "Z+" : "Z-" for j in 1:N]
    end
    psi = productMPS(s, new_state)
  end

  return nothing
end


text:
#
# Joe Fowler and Galen O'Neil
# NIST Boulder Laboratories
# December 2014 - February 2015
#

module DynamicTimeWarp

include("WindowedMatrix.jl")

# Dynamic Time Warping with a user-specified distance function

function dtw(seq1::Vector, seq2::Vector, distance::Function=Distance.square)

    # Build the cost matrix
    const m=length(seq2)
    const n=length(seq1)
    cost11 = distance(seq1[1], seq2[1])
    cost = zeros(typeof(cost11), m, n)

    # Initialize first column and first row
    cost[1,1] = cost11
    for r=2:m
        cost[r,1] = cost[r-1,1] + distance(seq1[1], seq2[r])
    end
    for c=2:n
        cost[1,c] = cost[1,c-1] + distance(seq1[c], seq2[1])
    end

    # Complete the cost matrix
    for c=2:n
        for r=2:m
            best_neighbor_cost = min(cost[r-1,c], cost[r-1,c-1], cost[r,c-1])
            cost[r,c] = best_neighbor_cost + distance(seq1[c], seq2[r])
        end
    end

    trackcols, trackrows = trackback(cost)
    cost[end,end], trackcols, trackrows
end

# Compute the optimal track backwards through the cost matrix from end to beginning.
# Return (columns, rows) of the optimal track.

function trackback(D)
    r,c = size(D)
    rows,cols = [r],[c]
    while r > 1 && c > 1
        tb = indmin([D[r-1,c-1], D[r-1,c], D[r,c-1]])
        tb in [1,2] && (r-=1)
        tb in [1,3] && (c-=1)
        push!(rows,r)
        push!(cols,c)
    end
    # Possibly either r>1 or c>1 at this point (but not both).
    # Add the unfinished part of the track to reach [1,1]
    for r=r-1:-1:1
        push!(rows,r)
        push!(cols,1)
    end
    for c=c-1:-1:1
        push!(rows,1)
        push!(cols,c)
    end
    reverse(cols), reverse(rows)
end

# The FastDTW approximation to the DTW, described in "FastDTW: Toward Accurate
# Dynamic Time Warping in Linear Time and Space", S Salvador & P Chan, __Intelligent
# Data Analysis__ (2007).

function fastdtw(seq1::Vector, seq2::Vector, radius::Integer,
                 distance::Function=Distance.square)
    const MinSize = max(radius + 2, 10)
    const N1 = length(seq1)
    const N2 = length(seq2)
    if N1 <= MinSize || N2 <= MinSize
        return (dtw(seq1, seq2, distance))
    end

    # Call recursively on a pair of sequences half this length
    compressed1 = compress(seq1)
    compressed2 = compress(seq2)
    _cost, lowrescol, lowresrow = fastdtw(compressed1, compressed2, radius, distance)

    # Now resample that path to the finer resolution, find the correct
    # window around it, and get the DTW given that window.
    hirescol, hiresrow = expandpath(lowrescol, lowresrow, N1, N2)
    idx2min, idx2max = computewindow(hirescol, hiresrow, radius)
    cost1, newcol, newrow = dtwwindowed(seq1, seq2, idx2min, idx2max, distance)
end

# Given a path through low-res space, generate an approximate path
# through high-res space. It should have dimension Ncol x Nrow

function expandpath(lowrescol, lowresrow, Ncol, Nrow)
    @assert div(Ncol+1,2) == lowrescol[end]
    @assert div(Nrow+1,2) == lowresrow[end]
    const Np = length(lowrescol)
    @assert Np == length(lowresrow)

    hirescol = zeros(eltype(lowrescol), 2*Np)
    hiresrow = zeros(eltype(lowresrow), 2*Np)
    hirescol[1] = hiresrow[1] = c = r = 1
    for i=1:Np-1
        # Select plan according to the next move in lowres path.
        if lowrescol[i+1] == lowrescol[i]  # Next move is up
            r += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            r += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r

        elseif lowresrow[i+1] == lowresrow[i] # Next move is sideways
            c += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            c += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r

        else  # Next move is diagonal.
            c += 1; r += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            c += 1; r += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r
        end
    end
    hirescol[end] = Ncol
    hiresrow[end] = Nrow
    # When expanding to an odd numbered size, it's possible to repeat
    # the last step.  Fix that:
    if hirescol[end]==hirescol[end-1] && hiresrow[end]==hiresrow[end-1]
        hirescol = hirescol[1:end-1]
        hiresrow = hiresrow[1:end-1]
    end
    hirescol, hiresrow
end

# Given the lists of (col,row) indices for the optimal path, compute a "window"
# around that path of the given radius.
# Returns (rowmin, rowmax), each a vector of length pathcols[end], representing
# for each column, the minimum and maximum row numbers used in that column.

function computewindow(pathcols, pathrows, radius)
    const Np = length(pathcols)
    @assert Np == length(pathrows)
    const Ncol = pathcols[end]
    const Nrow = pathrows[end]

    # Find the min/max row at each column in the path.
    pathmin = zeros(Int, Ncol)
    pathmax = zeros(Int, Ncol)
    for i=1:Np
        c,r = pathcols[i], pathrows[i]
        pathmax[c] = r
        if pathmin[c] == 0
            pathmin[c] = r
        end
    end

    # The window in each column for "radius" r starts at the pathmin
    # of the rth-previous column and ends at the pathmax of the
    # rth-next column, plus (in each case) the radius.
    if radius < Ncol-1 && radius < Nrow-1
        rowmin = vcat(fill(1,radius), pathmin[1:end-radius]-radius)
        rowmax = vcat(pathmax[radius+1:end]+radius, fill(Nrow,radius))

        # Window values must be in the range [1:Nrow].
        for c=1:Ncol
            if rowmin[c]<1; rowmin[c]=1; end
            if rowmax[c]>Nrow; rowmax[c]=Nrow; end
        end
    else
        rowmin = fill(1,Ncol)
        rowmax = fill(Nrow,Ncol)
    end
    rowmin, rowmax
end

# Do DTW in a subset of the full space, the subset specified by
# a "window". Arguments [idx2min,idx2max] give the inclusive lowest
# and highest index in the seq2 direction, one element for each index along
# the seq1 direction. Thus seq1, idx2min, and idx2max should all be of
# equal length

function dtwwindowed(seq1::Vector, seq2::Vector,
                     idx2min::Vector, idx2max::Vector,
                     distance::Function=Distance.square)

    const m=length(seq2) # of rows  in cost matrix
    const n=length(seq1) # of columns in cost matrix
    @assert n==length(idx2min)
    @assert n==length(idx2max)
    @assert 1==minimum(idx2min)
    @assert m==maximum(idx2max)

    # Build the (n x m) cost matrix into a WindowedMatrix, because it's ragged.
    # That type gives efficient storage with convenient [r,c] indexing and returns
    # Inf when accessed outside the window.
    cost = WindowedMatrix(idx2min, idx2max, Inf)

    # First column first
    cost[1,1] = distance(seq1[1], seq2[1])
    for r=2:idx2max[1]
        cost[r,1] = cost[r-1,1]  + distance(seq1[1], seq2[r])
    end

    # Complete the cost matrix from columns 2 to m.
    for c=2:n
        for r=idx2min[c]:idx2max[c]
            best_neighbor_cost = min(cost[r-1,c], cost[r-1,c-1], cost[r,c-1])
            cost[r,c] = best_neighbor_cost + distance(seq1[c], seq2[r])
        end
    end
    trackcols, trackrows = trackback(cost)
    cost[end,end], trackcols, trackrows
end

function compress(seq::Vector)
    Navg = div(length(seq), 2)
    evenseq = 0.5*(seq[1:2:end-1]+seq[2:2:end])
    if length(seq)%2 == 1
        return vcat(evenseq, [seq[end]])
    end
    evenseq
end

using PyPlot

function dtwspeedtest(N::Integer, functype=1)
    @assert N>10
    if functype == 1
        stretch = div(N,10)
        stretch = max(stretch, 10)
        t = linspace(0,1,N-stretch)
        x = vcat(fill(0.0, stretch), sin((t.^1.3).*5))
        y = vcat(sin((t.^0.7).*5), fill(sin(5), stretch))
    else
        # Random path
        rows=Array(Int, N)
        cols=Array(Int, N)
        r=c=1
        const PSTRETCH=0.8
        for i=1:N
            rnum = rand()
            if rnum > PSTRETCH # go diag
                r+=1; c+=1
            elseif rnum > 0.5*PSTRETCH
                r+=1
            else
                c+=1
            end
            rows[i] = r
            cols[i] = c
        end
        t = linspace(0, N, max(r[end],c[end]))
        f = sin((N-t).^0.7)+ (t.*(3./N))
        x = f[rows]
        y = f[cols]
    end

    clf()
    plt.subplot(211)
    dstd = std(x)+std(y)
    plot(x-mean(x), "r")
    plot(y-mean(y)+dstd, "b")

    plt.subplot(224)
    radii=[40,20,10,5,2]
    colors=["red","orange","gold","green","cyan"]
    for (rad,color) in zip(radii, colors)
        print("Testing size $N radius $rad: ")
        @time _,r1,c1=fastdtw(x,y,rad)
        plot(r1,c1,color=color)
    end

    if N<=2500
        @time _,r,c=dtw(x,y)
        plot(r,c,"k")
    end
    nothing
end

function plotdtw(seq1::Vector, seq2::Vector, offset=0.0)
    cost, match1, match2 = dtw(seq1, seq2)
    if offset == 0.0
        offset = 2*(std(seq1) + std(seq2)) + mean(seq1) - mean(seq2)
    end
    clf()
    plot(seq1, "-r", seq2+offset, "-b")
    for i=1:length(match1)
        plot([match1[i],match2[i]]-1, [seq1[match1[i]], seq2[match2[i]]+offset],
             color="gray")
    end
end

#function dtwbaryavg_iteration(dbavg::Vector, sequences::Array{Array{1},1})
function dtwbaryavg_iteration(dbavg::Vector, sequences::Array)
    const Nseq = length(sequences)
    count = zeros(Int, length(dbavg))
    sumcoords = zeros(Float64, length(dbavg))

    for i=1:Nseq
        cost, match1, match2 = dtw(dbavg, sequences[i])
        for j=1:length(match2)
            count[match1[j]] += 1
            sumcoords[match1[j]] += sequences[i][match2[j]]
        end
        println("Compared $i to the standard")
    end
    sumcoords = sumcoords ./ count
end

#function dtwbaryavg{T<:Real}(sequences::Array{Array{T,N},1})
function dtwbaryavg(sequences)
    dbavg = copy(sequences[1])
    clf()
    for i = 1:5
        dbavg = dtwbaryavg_iteration(dbavg, sequences)
        plot(dbavg)
    end
    dbavg
end

# This sub-module contains various functions returning pointwise distances
# between elements from two sequences.

module Distance
square(x,y) = (x-y)^2
absval(x,y) = abs(x-y)

# And now a Poisson-sensitive distance.
# Call this with the 2 sequences, and it returns a closure, a function that
# measures Poisson-distance for elements of those two particular sequences.
# Alternately, call it with two integers, representing the number of elements
# in each sequence.
function poissonpenalty(seq1::Vector, seq2::Vector)
    n1=sum(seq1)
    n2=sum(seq2)
    poissonpenalty(n1, n2)
end

function poissonpenalty(n1::Integer, n2::Integer)
    function pdist(x,y)
        if x<=0 && y<=0; return 0; end
        mu = float(x+y)/float(n1+n2)
        return ((x-mu*n1)^2 + (y-mu*n2)^2) / (x+y)
    end
end

end # module Distance

export
dtw,
dtwwindowed,
fastdtw,
dtwbaryavg

end # module


text:
export FacElem, FacElemMon, simplify, factor_coprime

export transform

################################################################################
#
#  Insert in a dictionary
#
################################################################################

function add_to_key!(D::Dict{S, T}, k::S, v; remove_zero::Bool = true) where S where T <: Union{fmpz, Integer}
  add_to_key!(D, k, T(v), remove_zero = remove_zero)
  return nothing
end

function add_to_key!(D::Dict{S, T}, k::S, v::T; remove_zero::Bool = true) where S where T <: Union{fmpz, Integer}
  hash_k = Base.ht_keyindex2!(D, k)
  if hash_k > 0
    #The key is in the dictionary, we only need to add
    w = D.vals[hash_k]
    if remove_zero
      if w != v && iszero(cmpabs(w, v))
        Base._delete!(D, hash_k)
      else
        @inbounds D.vals[hash_k] = w + v
      end
    else
      @inbounds D.vals[hash_k] = w + v
    end
  else
    pos = -hash_k
    @inbounds D.slots[pos] = 0x1
    @inbounds D.keys[pos] = k
    @inbounds D.vals[pos] = v
    D.count += 1
    if pos < D.idxfloor
      D.idxfloor = pos
    end
  end
  D.age += 1
  return nothing
end

function cmpabs(a::Int, b::Int)
  a = abs(a)
  b = abs(b)
  if a > b
    return 1
  elseif a == b
    return 0
  else
    return -1
  end
end
################################################################################
#
#  Multiplicative representation
#
################################################################################

# abstract nonsense

const FacElemMonDict = IdDict()

function (x::FacElemMon{S})() where S
  z = FacElem{elem_type(S), S}()
  z.fac = Dict{elem_type(S), fmpz}()
  z.parent = x
  return z
end

@doc Markdown.doc"""
    FacElem{B}(R, base::Vector{B}, exp::Vector{fmpz}) -> FacElem{B}

Returns the element $\prod b_i^{e_i}$, un-expanded.
"""
function FacElem(R, base::Vector{B}, exp::Vector{fmpz}) where {B}
  if elem_type(R) !== B
    throw(ArgumentError("Parent of elements wrong."))
  end

  length(base) != length(exp) && throw(error("not the same length"))

  z = FacElem{B, typeof(R)}()

  for i in 1:length(base)
    if iszero(exp[i])
      continue
    end
    add_to_key!(z.fac, base[i], exp[i], remove_zero = true)
  end

  z.parent = FacElemMon(R)
  return z
end

@doc Markdown.doc"""
    FacElem{B}(base::Vector{B}, exp::Vector{fmpz}) -> FacElem{B}

Returns the element $\prod b_i^{e_i}$, un-expanded.
"""
function FacElem(base::Vector{B}, exp::Vector{fmpz}) where B

  length(base) == 0 && error("Array must not be empty")

  R = parent(base[1])

  return FacElem(R, base, exp)
end

@doc Markdown.doc"""
    FacElem{B}(R, d::Dict{B, fmpz}) -> FacElem{B}
    FacElem{B}(R, d::Dict{B, Integer}) -> FacElem{B}

Returns the element $\prod b^{d[p]}$, un-expanded.
"""
function FacElem(R, d::Dict{B, fmpz}) where B
  if elem_type(R) !== B
    throw(ArgumentError("Parent of elements wrong."))
  end

  z = FacElem{B, typeof(R)}()
  z.fac = d
  z.parent = FacElemMon(R)
  return z
end

@doc Markdown.doc"""
    FacElem{B}(d::Dict{B, fmpz}) -> FacElem{B}
    FacElem{B}(d::Dict{B, Integer}) -> FacElem{B}

Returns the element $\prod b^{d[p]}$, un-expanded.
"""
function FacElem(d::Dict{B, fmpz}) where B

  length(d) == 0 && error("Dictionary must not be empty")

  R = parent(first(keys(d)))

  return FacElem(R, d)
end

function FacElem(R::Ring)
  z = FacElem{elem_type(R), typeof(R)}()
  z.fac = Dict{elem_type(R), fmpz}()
  z.parent = FacElemMon(R)
  return z
end

function FacElem(R, d::Dict{B, T}) where {B, T <: Integer}

  z = FacElem{B, typeof(R)}()
  z.fac = Dict{B, fmpz}((k,fmpz(v)) for (k,v) = d)

  z.parent = FacElemMon(R)
  return z
end

function FacElem(d::Dict{B, T}) where {B, T <: Integer}

  length(d) == 0 && error("Dictionary must not be empty")

  R = parent(first(keys(d)))

  return FacElem(R, d)
end

parent(x::FacElem) = x.parent

base_ring(x::FacElemMon) = x.base_ring

base_ring(x::FacElem) = base_ring(parent(x))

base(x::FacElem) = keys(x.fac)

function Base.deepcopy_internal(x::FacElem{B, S}, dict::IdDict) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.deepcopy_internal(x.fac, dict)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
#  z.hash = x.hash # this needs to be deleted in ! operations
  return z
end

function Base.copy(x::FacElem{B, S}) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.copy(x.fac)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
#  z.hash = x.hash
  return z
end

Base.iterate(a::FacElem) = Base.iterate(a.fac)

Base.iterate(a::FacElem, state) = Base.iterate(a.fac, state)

Base.length(a::FacElem) = Base.length(a.fac)

check_parent(x::FacElem{S, T}, y::FacElem{S, T}) where { S, T } = base_ring(x) == base_ring(y)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FacElemMon)
  print(io, "Factored elements over $(x.base_ring)")
end

function show(io::IO, x::FacElem)
  print(io, "Factored element with data\n$(x.fac)")
end

################################################################################
#
#  Inverse
#
################################################################################
function inv!(x::FacElem)
  for i = x.fac.idxfloor:length(x.fac.vals)
    if isassigned(x.fac.vals, i)
      x.fac.vals[i] = -x.fac.vals[i]
    end
  end
  return x
end

function inv(x::FacElem)
  y = copy(x)
  inv!(y)
  return y
end

################################################################################
#
#  Exponentiation
#
################################################################################

function pow!(z::FacElem, x::FacElem, y::T) where T <: Union{fmpz, Integer}
  z.fac = copy(x.fac)
  for i = z.fac.idxfloor:length(z.fac.vals)
    if isassigned(z.fac.vals, i)
      z.fac.vals[i] = z.fac.vals[i]*y
    end
  end
  return nothing
end

function pow!(z::FacElem, y::T) where T <: Union{fmpz, Integer}
  for i = z.fac.idxfloor:length(z.fac.vals)
    if isassigned(z.fac.vals, i)
      z.fac.vals[i] = z.fac.vals[i]*y
    end
  end
end

# ^(x::FacElem, y::Union{fmpz, Integer}) is ambiguous
for T in [:Integer, fmpz]
  @eval begin
    function ^(x::FacElem, y::($T))
      z = parent(x)()
      if y == 0
        return z
      end
      if y == 1
        return copy(x)
      else
        z.fac = copy(x.fac)
        for i = z.fac.idxfloor:length(z.fac.vals)
          if isassigned(z.fac.vals, i)
            z.fac.vals[i] = z.fac.vals[i]*y
          end
        end
        return z
      end
    end
  end
end

################################################################################
#
#  Multiplication
#
################################################################################

function mul!(z::FacElem{B, S}, x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  @assert check_parent(x, y) "Elements must have the same parent"
  @assert check_parent(x, z) "Elements must have the same parent"
  z.fac = copy(x.fac)
  for (a, v) in y
    add_to_key!(z.fac, a, v, remove_zero = true)
  end
  return z
end

function *(x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  @assert check_parent(x, y) "Elements must have the same parent"
  if length(x.fac) == 0
    return copy(y)
  end

  if length(y.fac) == 0
    return copy(x)
  end

  z = copy(x)
  for (a, v) in y
    add_to_key!(z.fac, a, v, remove_zero = true)
  end

  return z
end

function *(x::FacElem{B}, y::B) where B
  @assert base_ring(x) == parent(y)
  z = copy(x)
  add_to_key!(z.fac, y, 1, remove_zero = true)
  return z
end

function *(y::B, x::FacElem{B}) where B
  @assert base_ring(x) == parent(y)
  z = copy(x)
  add_to_key!(z.fac, y, 1, remove_zero = true)
  return z
end

function div(x::FacElem{B}, y::FacElem{B}) where B
  @assert check_parent(x, y) "Elements must have the same parent"
  z = copy(x)
  for (a, v) in y
    add_to_key!(z.fac, a, -v, remove_zero = true)
  end
  return z
end

################################################################################
#
#  Apply transformation
#
################################################################################

# return (x1,...,xr)*y
function _transform(x::Vector{FacElem{T, S}}, y::fmpz_mat) where {T, S}
  length(x) != nrows(y) &&
              error("Length of array must be number of rows of matrix")

  z = Array{FacElem{T, S}}(undef, ncols(y))

  t = parent(x[1])()

  for i in 1:ncols(y)
    z[i] = x[1]^y[1,i]
    for j in 2:nrows(y)
       if iszero(y[j, i])
        continue
      end
      if isone(y[j, i])
        mul!(z[i], z[i], x[j])
      else
        pow!(t, x[j], y[j, i])
        mul!(z[i], z[i], t)
      end
    end
  end
  return z
end

function transform(x::Vector{FacElem{S, T}}, y::fmpz_mat) where {S, T}
  return _transform(x, y)
end

function transform_left!(x::Vector{FacElem{S, T}}, y::TrafoSwap{fmpz}) where {S, T}
  x[y.i], x[y.j] = x[y.j], x[y.i]
  nothing
end

function transform_left!(x::Vector{FacElem{S, T}}, y::TrafoAddScaled{fmpz}) where {S, T}
  x[y.j] = x[y.j] * x[y.i]^y.s
  nothing
end

function transform_left!(x::Vector{FacElem{S, T}}, y::TrafoPartialDense{R}) where {S, T, R}
  z = view(deepcopy(x), y.rows)
  xx = view(x, y.rows)
  for i in 1:nrows(y.U)
    xx[i] = z[1]^Int(y.U[i, 1])
    for j in 2:ncols(y.U)
      xx[i] *= z[j]^Int(y.U[i, j])
    end
  end
end

function transform_left!(x::Vector{T}, t::TrafoDeleteZero{S}) where {S, T}
  # move ith position to the back
  for j in t.i:length(x)-1
    r = x[j]
    x[j] = x[j + 1]
    x[j + 1] = r
  end
end

function transform_left!(x::Vector{FacElem{S, T}}, y::TrafoParaAddScaled{R}) where {S, T, R}
  # [ Ri; Rj] <- [a b; c d]* [ Ri; Rj ]
  ri = deepcopy(x[y.i])
  rj = deepcopy(x[y.j])
  x[y.i] = ri^y.a * rj^y.b
  x[y.j] = ri^y.c * rj^y.d
end

################################################################################
#
#  Evaluate
#
################################################################################

function evaluate(x::FacElem{NfOrdIdl, NfOrdIdlSet}; coprime::Bool = false)
  O = order(base_ring(x))
  if !coprime
    x = simplify(x) # the other method won't work due to one()
  end
  if length(x.fac)==0
    return fractional_ideal(O, O(1))
  end
  # still doesn't work
  D = collect(x.fac)
  A = (D[1][1]//1)^Int(D[1][2])
  for i in 2:length(D)
    A = A * (D[i][1]//1)^Int(D[i][2])
  end
  return A
end

function _ev(d::Dict{nf_elem, fmpz}, oe::nf_elem)
  z = deepcopy(oe)
  if length(d)==0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^x[2]
  end
  b = empty(d)
  for (k, v) in d
    if iszero(v)
      continue
    end
    if v>-10 && v<10
      mul!(z, z, k^Int(v))
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv != 0
        b[k] = vv
      end
      if r != 0
        mul!(z, z, k)
      end
    end
  end
  if isempty(b)
    return z
  end
  res = _ev(b, oe)
  mul!(res, res, res)
  mul!(z, z, res)
  return z
end

function _ev(d::Dict{fq_nmod, fmpz}, z::fq_nmod)
  Fq = parent(z)
  if length(d) == 0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^Int(x[2])
  end
  b = empty(d)
  for (k, v) in d
    if iszero(v)
      continue
    end
    if abs(v) < 10
      if v >0
        kv = Fq()
        ccall((:fq_nmod_pow, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), kv, k, v, Fq)
        mul!(z, z, kv)
      else
        kv = Fq()
        ccall((:fq_nmod_inv, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), kv, k, Fq)
        ccall((:fq_nmod_pow, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), kv, kv, -v, Fq)
        mul!(z, z, kv)
      end
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv != 0
        b[k] = vv
      end
      if r != 0
        mul!(z, z, k)
      end
    end
  end
  if isempty(b)
    return z
  end
  res = _ev(b, one(Fq))
  mul!(res, res, res)
  mul!(z, z, res)
  return z
end

function _ev(d::Dict{T, fmpz}, oe::T) where T
  z = copy(oe)
  if length(d)==0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^x[2]
  end
  b = empty(d)
  for (k,v) in d
    if v>-10 && v<10
      z *= k^Int(v)
    else
      r = isodd(v) ? 1 : 0
      vv = div(v-r, 2)
      if vv!=0
        b[k] = vv
      end
      if r!=0
        z*= k
      end
    end
  end
  return _ev(b, oe)^2*z
end

@doc Markdown.doc"""
    evaluate{T}(x::FacElem{T}) -> T

Expands or evaluates the factored element, i.e. actually computes the
value.
Does "square-and-multiply" on the exponent vectors.
"""
function evaluate(x::FacElem{T}) where T
  return _ev(x.fac, one(base_ring(x)))
end

@doc Markdown.doc"""
    evaluate_naive{T}(x::FacElem{T}) -> T

Expands or evaluates the factored element, i.e. actually computes the
value. Uses the obvious naive algorithm. Faster for input in finite rings.
"""
function evaluate_naive(x::FacElem{T}) where T
  z = one(base_ring(x))
  d = x.fac
  if length(d)==0
    return z
  end
  for (k,v) in d
    z *= k^v
  end
  return z
end

#################################################################################
@doc Markdown.doc"""
    max_exp(a::FacElem)

Finds the largest exponent in the factored element $a$.
"""
function max_exp(a::FacElem)
  return maximum(values(a.fac))
end

@doc Markdown.doc"""
    min_exp(a::FacElem)

Finds the smallest exponent in the factored element $a$.
"""
function min_exp(a::FacElem)
  return minimum(values(a.fac))
end

@doc Markdown.doc"""
    maxabs_exp(a::FacElem)

Finds the largest exponent by absolute value in the factored element $a$.
"""
function maxabs_exp(a::FacElem)
  return maximum(abs, values(a.fac))
end

function Base.hash(a::FacElem, u::UInt)
  if a.hash == UInt(0)
    h = hash(u, UInt(3127))
    for (k,v) = a.fac
      h = hash(k, hash(v, h))
    end
    a.hash = h
  end
  return a.hash
end

#################################################################################
##
##  Auxillary deep copy functions
##
#################################################################################
#
#function _deepcopy{K, V}(x::Dict{K, V})
#  z = typeof(x)()
#  for a in keys(x)
#    z[a] = deepcopy(x[a])
#  end
#  return z
#end

################################################################################
#
#  Parent object overloading
#
################################################################################

function (F::FacElemMon)(a::T) where T
  z = F()
  z.fac[a] = fmpz(1)
  return z
end

#function order(A::FacElemMon{IdealSet})
#  return order(A.base_ring)
#end

function order(A::FacElemMon)
  return order(A.base_ring)
end

(::Type{FacElem{T}})(a::FacElem{T}) where {T} = a

(::Type{FacElem})(a::FacElem{T}) where {T} = a


text:
module WaveFunctionCollapse

include("overlapping_model.jl")

end


text:
@with_kw struct DGaussianPulse <: AbstractWaveform
    fc::Float64 = 0.0
    half_bandwidth::Float64
    half_width::Float64 = 1 / (2 * sqrt(log(2))) * 1 / (π * half_bandwidth)
    time_delay::Float64 = 6 * half_width
end
(s::DGaussianPulse)(t) = -(t - s.time_delay) / s.half_width * exp(-(t - s.time_delay)^2 / s.half_width^2)
@with_kw struct GaussianPulse <: AbstractWaveform
    fc::Float64 = 0.0
    half_bandwidth::Float64
    half_width::Float64 = 1 / (2 * sqrt(log(2))) * 1 / (π * half_bandwidth)
    time_delay::Float64 = 6 * half_width
end

(s::GaussianPulse)(t) = cos(2π * s.fc * (t - s.time_delay)) * exp(-(t - s.time_delay)^2 / s.half_width^2)

const GaussianWaveForm = Union{GaussianPulse,DGaussianPulse}

@recipe function f(s::GaussianWaveForm; nsamples = 20, spectrum = true)
    if spectrum
        xlabel := "frequency"
        xrange := [-s.fc - 10 * s.half_bandwidth, s.fc + 10 * s.half_bandwidth]
    else
        xrange := [s.time_delay - 3s.half_width, s.time_delay + 3s.half_width]
    end
    δt = 1 / (s.fc + s.half_bandwidth) / nsamples
    t_samples = range(s.time_delay - 20s.half_width, s.time_delay + 20s.half_width, step = δt)
    if spectrum
        return fourier(t_samples, s.(t_samples))
    else
        return (t_samples, s.(t_samples))
    end
end


text:
"""
Test minimize function with no constraints
"""

#push!(LOAD_PATH,"../src/")
using joptimise

function rosenbrock(x)
    f = (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
    return f
end

function rosenbrock!(g, x)
    # compute objective
    f = rosenbrock(x)
    g[1] = 0.0
    return f
end

# initial guess
x0 = [4.0; 4.0]
# bounds on variables
lx = [-5.0; -5.0]
ux = [5.0; 5.0]
# bounds on constriants
lg = [0.0]
ug = [0.0]
# number of constraints
ng = 1

## run minimizer with IPOPT
ip_options = Dict(
    "max_iter" => 2500,   # 1500 ~ 2500
    "tol" => 1e-6
)

xopt, fopt, info = minimize(rosenbrock!, x0, ng; lx=lx, ux=ux, lg=lg, ug=ug, solver="ipopt", options=ip_options, derivatives=ForwardAD());

println("Done with IPOPT!")
println(info)
println(xopt)

## run minimizer with SNOPT
sn_options = Dict(
    "Major feasibility tolerance" => 1.e-6,
    "Major optimality tolerance"  => 1.e-6,
    "Minor feasibility tolerance" => 1.e-6,
    "Major iterations limit" => 1000,
    "Major print level" => 1,
)

xopt, fopt, info = minimize(rosenbrock!, x0, ng; lx=lx, ux=ux, lg=lg, ug=ug, solver="snopt", options=sn_options, derivatives=ForwardAD());

println("Done with SNOPT!")
println(info)
println(xopt)

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
using Pkg

Pkg.add("JLD")
Pkg.add("JLD2")
Pkg.add("FFTW")
Pkg.add("Trapz")
Pkg.add("Printf")
Pkg.add("FileIO")
Pkg.add("PyPlot")
Pkg.add("Random")
Pkg.add("Statistics")
Pkg.add("Parameters")
Pkg.add("SparseArrays")
Pkg.add("ProgressMeter")
Pkg.add("Distributions")
Pkg.add("LinearAlgebra")
Pkg.add("DelimitedFiles")

max_stars_repo_name: raminkhajeh/sparse-balance

max_issues_repo_name: raminkhajeh/sparse-balance

max_forks_repo_name: raminkhajeh/sparse-balance


text:
module ToggleMenuItems

Base.Experimental.@optlevel 2

using Redux
using CImGui
using ..MenuItems
import ..MenuItems: AbstractMenuItemAction, AbstractGenericMenuItem, AbstractMenuItem,
                    MenuItem, get_label

# actions
const Rename = MenuItems.Rename
const SetTriggeredTo = MenuItems.SetTriggeredTo
const ChangeShortcut = MenuItems.ChangeShortcut
const Enable = MenuItems.Enable
const Disable = MenuItems.Disable

"""
    SetSelectedTo(label, is_selected)
Set widget selected state to `is_selected`.
"""
struct SetSelectedTo <: AbstractMenuItemAction
    label::String
    is_selected::Bool
end

"""
    Toggle(label)
Toggle the menu item `is_selected` state.
"""
struct Toggle <: AbstractMenuItemAction
    label::String
end

"""
    AddMenuItem(label::AbstractString, shortcut::AbstractString = "", is_selected = false, is_enabled = true)
Add a menu item.
"""
struct AddMenuItem <: AbstractMenuItemAction
    menu_item::MenuItem
    is_selected::Bool
end
AddMenuItem(label::AbstractString, shortcut::AbstractString = "", is_selected = false, is_enabled = true) =
    AddMenuItem(MenuItem(label, shortcut, is_enabled, false), is_selected)

# any subtypes of AbstractMenuItemAction should have a `label` field, so we fake it.
function Base.getproperty(obj::AddMenuItem, sym::Symbol)
    if sym === :label
        return obj.menu_item.label
    else
        return getfield(obj, sym)
    end
end

"""
    DeleteMenuItem(label::AbstractString)
Delete the menu item.
"""
struct DeleteMenuItem <: AbstractMenuItemAction
    label::String
end

# state
"""
    ToggleMenuItem(label::AbstractString, shortcut::AbstractString = "", is_selected = false, is_enabled = true)
A [`MenuItem`](@ref) with a toggle state.
"""
struct ToggleMenuItem <: AbstractMenuItem
    menu_item::MenuItem
    is_selected::Bool
end
ToggleMenuItem(label::AbstractString, shortcut::AbstractString = "", is_selected = false, is_enabled = true) =
    ToggleMenuItem(MenuItem(label, shortcut, is_enabled, false), is_selected)

# any subtypes of AbstractMenuItem should have a `label` field, so we fake it.
function Base.getproperty(obj::ToggleMenuItem, sym::Symbol)
    if sym === :label
        return obj.menu_item.label
    else
        return getfield(obj, sym)
    end
end

get_shortcut(s::ToggleMenuItem) = s.menu_item.shortcut
is_triggered(s::ToggleMenuItem) = s.menu_item.is_triggered
is_enabled(s::ToggleMenuItem) = s.menu_item.is_enabled
is_selected(s::ToggleMenuItem) = s.is_selected

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::ToggleMenuItem, a::AbstractMenuItemAction) = ToggleMenuItem(MenuItems.reducer(s.menu_item, a), s.is_selected)
reducer(s::ToggleMenuItem, a::Toggle) = ToggleMenuItem(s.menu_item, !s.is_selected)
reducer(s::ToggleMenuItem, a::SetSelectedTo) = ToggleMenuItem(s.menu_item, a.is_selected)

reducer(s::Dict{String,<:AbstractImmutableState}, a::AbstractMenuItemAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractImmutableState}, a::AbstractMenuItemAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end
reducer(s::Vector{<:AbstractImmutableState}, a::AddMenuItem) =
    [s..., ToggleMenuItem(a.menu_item, a.is_selected)]
reducer(s::Vector{<:AbstractImmutableState}, a::DeleteMenuItem) =
    filter(s -> get_label(s) !== get_label(a), s)

# UI
"""
    (::ToggleMenuItem)(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered/activated.
`get_state` is a router function that tells how to find the target state from `store`.
`chain_action` is for chaining upstream actions.
"""
function (::ToggleMenuItem)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity)
    s = get_state(store)
    is_activated = CImGui.MenuItem(get_label(s), get_shortcut(s), s.is_selected, is_enabled(s))
    dispatch!(
        store,
        SetTriggeredTo(get_label(s), is_activated) |> chain_action,
    )
    is_activated && dispatch!(
        store,
        Toggle(get_label(s)) |> chain_action,
    )
    return is_activated
end

end # module


text:
##
# Sparse BroadcastStyle
##

for Typ in (:Diagonal, :SymTridiagonal, :Tridiagonal)
    @eval begin
        BroadcastStyle(::StructuredMatrixStyle{<:$Typ}, ::BandedStyle) =
            BandedStyle()
        BroadcastStyle(::BandedStyle, ::StructuredMatrixStyle{<:$Typ}) =
            BandedStyle()
    end
end

# Here we implement the banded matrix interface for some key examples
isbanded(::Zeros) = true
bandwidths(::Zeros) = (0,0)
inbands_getindex(::Zeros{T}, k::Integer, j::Integer) where T = zero(T)

isbanded(::Eye) = true
bandwidths(::Eye) = (0,0)
inbands_getindex(::Eye{T}, k::Integer, j::Integer) where T = one(T)

isbanded(::Diagonal) = true
bandwidths(::Diagonal) = (0,0)
inbands_getindex(D::Diagonal, k::Integer, j::Integer) = D.diag[k]
inbands_setindex!(D::Diagonal, v, k::Integer, j::Integer) = (D.diag[k] = v)

isbanded(::SymTridiagonal) = true
bandwidths(::SymTridiagonal) = (1,1)
inbands_getindex(J::SymTridiagonal, k::Integer, j::Integer) =
    k == j ? J.dv[k] : J.ev[min(k,j)]
inbands_setindex!(J::SymTridiagonal, v, k::Integer, j::Integer) =
    k == j ? (J.dv[k] = v) : (J.ev[min(k,j)] = v)

isbanded(::Tridiagonal) = true
bandwidths(::Tridiagonal) = (1,1)

isbanded(K::Kron{<:Any,2}) = all(isbanded, K.arrays)
function bandwidths(K::Kron{<:Any,2})
    A,B = K.arrays
    (size(B,1)*bandwidth(A,1) + max(0,size(B,1)-size(B,2))*size(A,1)   + bandwidth(B,1),
        size(B,2)*bandwidth(A,2) + max(0,size(B,2)-size(B,1))*size(A,2) + bandwidth(B,2))
end

const BandedMatrixTypes = (:AbstractBandedMatrix, :(AdjOrTrans{<:Any,<:AbstractBandedMatrix}),
                                    :(AbstractTriangular{<:Any, <:AbstractBandedMatrix}),
                                    :(Symmetric{<:Any, <:AbstractBandedMatrix}))

const OtherBandedMatrixTypes = (:Zeros, :Eye, :Diagonal, :SymTridiagonal)

for T1 in BandedMatrixTypes, T2 in BandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

for T1 in BandedMatrixTypes, T2 in OtherBandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

for T1 in OtherBandedMatrixTypes, T2 in BandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end

###
# Specialised multiplication for arrays padded for zeros
# needed for ∞-dimensional banded linear algebra
###

function _copyto!(::VcatLayout{<:Tuple{<:Any,ZerosLayout}}, y::AbstractVector,
                 M::MatMulVec{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}})
    A,x = M.args
    length(y) == size(A,1) || throw(DimensionMismatch())
    length(x) == size(A,2) || throw(DimensionMismatch())

    ỹ,_ = y.arrays
    x̃,_ = x.arrays

    length(ỹ) ≥ min(length(M),length(x̃)+bandwidth(A,1)) ||
        throw(InexactError("Cannot assign non-zero entries to Zero"))

    ỹ .= Mul(view(A, axes(ỹ,1), axes(x̃,1)) , x̃)
    y
end

function similar(M::MatMulVec{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}}, ::Type{T}) where T
    A,x = M.args
    xf,_ = x.arrays
    n = max(0,min(length(xf) + bandwidth(A,1),length(M)))
    Vcat(Vector{T}(undef, n), Zeros{T}(size(A,1)-n))
end

###
# MulMatrix
###

bandwidths(M::MulMatrix) = bandwidths(M.applied)
isbanded(M::MulMatrix) = all(isbanded, M.applied.args)

const MulBandedMatrix{T} = MulMatrix{T, <:Mul{<:LayoutApplyStyle{<:Tuple{Vararg{<:AbstractBandedLayout}}}}}

BroadcastStyle(::Type{<:MulBandedMatrix}) = BandedStyle()

Base.replace_in_print_matrix(A::MulBandedMatrix, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

function _banded_mul_getindex(::Type{T}, (A, B), k::Integer, j::Integer) where T
    Al, Au = bandwidths(A)
    Bl, Bu = bandwidths(B)
    n = size(A,2)
    ret = zero(T)
    for ν = max(1,j-Bu,k-Al):min(n,j+Bl,k+Au)
        ret += A[k,ν] * B[ν,j]
    end
    ret
end

getindex(M::MatMulMat{<:AbstractBandedLayout,<:AbstractBandedLayout}, k::Integer, j::Integer) =
    _banded_mul_getindex(eltype(M), M.args, k, j)

getindex(M::Mul{<:LayoutApplyStyle{<:Tuple{Vararg{<:AbstractBandedLayout}}}}, k::Integer, j::Integer) =
    _banded_mul_getindex(eltype(M), (first(M.args), Mul(tail(M.args)...)), k, j)

@inline _sub_materialize(::MulLayout{<:Tuple{Vararg{<:AbstractBandedLayout}}}, V) = BandedMatrix(V)

MemoryLayout(V::SubArray{T,2,<:MulBandedMatrix,I}) where {T,I<:Tuple{Vararg{AbstractUnitRange}}} =
    MemoryLayout(parent(V))

@inline getindex(A::MulBandedMatrix, kr::Colon, jr::Colon) = _lazy_getindex(A, kr, jr)
@inline getindex(A::MulBandedMatrix, kr::Colon, jr::AbstractUnitRange) = _lazy_getindex(A, kr, jr)
@inline getindex(A::MulBandedMatrix, kr::AbstractUnitRange, jr::Colon) = _lazy_getindex(A, kr, jr)
@inline getindex(A::MulBandedMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = _lazy_getindex(A, kr, jr)


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

module TimeVectors

using Base.Dates, TimeSeries, FunctionalCollections

export TimeVector,
       timevector

include("timevector.jl")

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Advisory reentrant lock
"""
    ReentrantLock()

Creates a reentrant lock for synchronizing Tasks.
The same task can acquire the lock as many times as required.
Each `lock` must be matched with an `unlock`.

This lock is NOT threadsafe. See `Threads.Mutex` for a threadsafe lock.
"""
type ReentrantLock
    locked_by::Nullable{Task}
    cond_wait::Condition
    reentrancy_cnt::Int

    ReentrantLock() = new(nothing, Condition(), 0)
end

"""
    islocked(the_lock) -> Status (Boolean)

Check whether the lock is held by any task/thread.
This should not be used for synchronization (see instead `trylock`).
"""
function islocked(rl::ReentrantLock)
    return rl.reentrancy_cnt != 0
end

"""
    trylock(the_lock) -> Success (Boolean)

Acquires the lock if it is available,
returning `true` if successful.
If the lock is already locked by a different task/thread,
returns `false`.

Each successful `trylock` must be matched by an `unlock`.
"""
function trylock(rl::ReentrantLock)
    t = current_task()
    if rl.reentrancy_cnt == 0
        rl.locked_by = t
        rl.reentrancy_cnt = 1
        return true
    elseif t == get(rl.locked_by)
        rl.reentrancy_cnt += 1
        return true
    end
    return false
end

"""
    lock(the_lock)

Acquires the lock when it becomes available.
If the lock is already locked by a different task/thread,
it waits for it to become available.

Each `lock` must be matched by an `unlock`.
"""
function lock(rl::ReentrantLock)
    t = current_task()
    while true
        if rl.reentrancy_cnt == 0
            rl.locked_by = t
            rl.reentrancy_cnt = 1
            return
        elseif t == get(rl.locked_by)
            rl.reentrancy_cnt += 1
            return
        end
        wait(rl.cond_wait)
    end
end

"""
    unlock(the_lock)

Releases ownership of the lock.

If this is a recursive lock which has been acquired before, it
just decrements an internal counter and returns immediately.
"""
function unlock(rl::ReentrantLock)
    if rl.reentrancy_cnt == 0
        error("unlock count must match lock count")
    end
    rl.reentrancy_cnt -= 1
    if rl.reentrancy_cnt == 0
        rl.locked_by = nothing
        notify(rl.cond_wait)
    end
    return
end

function lock(f, l)
    lock(l)
    try
        return f()
    finally
        unlock(l)
    end
end

function trylock(f, l)
    if trylock(l)
        try
            return f()
        finally
            unlock(l)
        end
    end
    return false
end

"""
    Semaphore(sem_size)

Creates a counting semaphore that allows at most `sem_size`
acquires to be in use at any time.
Each acquire must be mached with a release.

This construct is NOT threadsafe.
"""
type Semaphore
    sem_size::Int
    curr_cnt::Int
    cond_wait::Condition
    Semaphore(sem_size) = sem_size > 0 ? new(sem_size, 0, Condition()) : throw(ArgumentError("Semaphore size must be > 0"))
end

"""
    acquire(s::Semaphore)

Wait for one of the `sem_size` permits to be available,
blocking until one can be acquired.
"""
function acquire(s::Semaphore)
    while true
        if s.curr_cnt < s.sem_size
            s.curr_cnt = s.curr_cnt + 1
            return
        else
            wait(s.cond_wait)
        end
    end
end

"""
    release(s::Semaphore)

Return one permit to the pool,
possibly allowing another task to acquire it
and resume execution.
"""
function release(s::Semaphore)
    @assert s.curr_cnt > 0 "release count must match acquire count"
    s.curr_cnt -= 1
    notify(s.cond_wait; all=false)
end


text:
using JSON

for exercise in readdir(joinpath("exercises", "concept"))
    exercise_path = joinpath("exercises", "concept", exercise)
    isdir(exercise_path) || continue
    cfg_path = joinpath(exercise_path, ".meta", "config.json")

    cfg = JSON.parsefile(cfg_path, use_mmap=false) # use_mmap=false makes it work on Windows
    if haskey(cfg, "editor")
        @info "⏩ $exercise already contains editor config settings"
        continue
    end

    cfg["editor"] = Dict(
        "test_files" => ["runtests.jl"],
        "solution_files" => ["$exercise.jl"],
    )

    open(cfg_path, "w") do io
        JSON.print(io, cfg, 4)
    end
    @info "✔ Added editor settings for $exercise" cfg["editor"]
end


text:
import PAPI

function slow_code()
    tmp = 1.1
    for i=1:200_000_000
        tmp = (tmp+100)/i
    end
    return tmp
end

function main()
    precompile(slow_code, ())

    flpops, mflops, real_time, proc_time  = @PAPI.flops slow_code()

    @printf("Real_time: %f Proc_time: %f Total flpops: %lld MFLOPS: %f\n",
             real_time, proc_time, flpops, mflops)
end

main()


text:
# =========================================================================== #
# TimeInterval
# =========================================================================== #

struct TimeInterval
    t_start::Float64
    t_end::Float64
end

function TimeInterval(t_end)
    TimeInterval(0.0, t_end)
end

function TimeInterval(x::Vector{T}) where {T}
    TimeInterval(x[1], x[2])
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export write_fileheader

function write_fileheader(s::IO, fh::FileHeader)

    # Check text header length
    if sizeof(fh.th) > 3200 @error "Text Header longer than 3200 bytes" end

    # Check dsf
    if fh.bfh.DataSampleFormat != 5
        @warn "DataSampleFormat not supported for writing. Attempting to convert to IEEE Float32"
        fh.bfh.DataSampleFormat = 5
    end

    ##0000
    # Write Text header #DEVNOTE# Big Endian?
    th_length = write(s, fh.th)

    # If too short, pad with blanks
    if th_length < 3200
        pad = 3200 - th_length
        write(s, " "^pad)
    end

    ##3200
    # Write first section of assigned values
    for field in fieldnames(typeof(fh.bfh))[1:27]
        write(s, bswap(getfield(fh.bfh, field)))
    end

    ##3260
    # Skip to next block
    write(s, Array{UInt8,1}(undef,240))

    ##3500
    # Write second section of assigned values
    for field in fieldnames(typeof(fh.bfh))[28:end]
        write(s, bswap(getfield(fh.bfh, field)))
    end

    ##3506
    # Skip to end
    write(s, Array{UInt8,1}(undef,94))

    ##3600

end


text:
#=
    The File contains `CSR_Sparsik<T>` and related funcs implementation

    `CSR_Sparsik` implements `AbstractSparseMatrix` interface
=#

#------------------------------------------------------------------------------

# It is convenient to have 0-based indexing in this file

ZeroBased(arr::Array) = OffsetArray(arr, Origin(0))
ZeroBased(T::Type) = OffsetArray(T[], Origin(0))

#------------------------------------------------------------------------------

# T
mutable struct CSR_Sparsik{T <: AbstractAlgebra.Field} <: AbstractSparseMatrix{T}
    # The struct implements the `Compressed Sparse Row` storage strategy
    # No amortized algorithms here!

    m::Int
    n::Int

    field::T

    #
    #
    row_ptr::OffsetVector{Int}
    col_idx::OffsetVector{Int}
    val::OffsetVector{<:AbstractAlgebra.FieldElem}
end

#------------------------------------------------------------------------------

nnz(A::CSR_Sparsik) = length(A.val)

Base.size(A::CSR_Sparsik) = (A.m, A.n)
Base.size(A::CSR_Sparsik, i::Int) = (A.m, A.n)[i]

issquare(A::CSR_Sparsik) = A.n == A.m

order(A::CSR_Sparsik) = A.n

base_ring(v::CSR_Sparsik) = v.field

#------------------------------------------------------------------------------

# constructs a CSR_Sparsik object from the given dense matrix `A`
# while coercing each element into the `field`
function dense_to_CSR(A::AbstractMatrix, field)
    nnz_coords = [
        (idx[1], idx[2], A[idx])
        for idx in findall(!iszero, A)
    ]
    return COO_to_CSR(size(A)..., nnz_coords, field)
end

#------------------------------------------------------------------------------

# constructs a CSR_Sparsik object from the given list of entries of form:
#   nnz_coords = { (x, y, val) }
# where (x, y) stands for the coordinates and val - for the value
# of an entry
function COO_to_CSR(m, n, nnz_coords, field)
    row_ptr = ZeroBased(Int[0])
    col_idx = ZeroBased(Int)
    data = ZeroBased(elem_type(field))

    sort!(nnz_coords)

    j = 1

    for i in 1 : m
        while j <= length(nnz_coords) && nnz_coords[j][1] == i
            push!(col_idx, nnz_coords[j][2] - 1)
            push!(data, field(nnz_coords[j][3]))
            j += 1
        end
        push!(row_ptr, j - 1)
    end

    return CSR_Sparsik(m, n, field, row_ptr, col_idx, data)
end

#------------------------------------------------------------------------------

function __apply_vector_seq(A::CSR_Sparsik, x)
    field = base_ring(A)
    nnz = Int[]
    data = Dict{Int, elem_type(field)}()

    for i in 1 : A.m
        v = zero(field)
        for j in A.row_ptr[i] : A.row_ptr[i + 1] - 1
            v += A.val[j] * x[A.col_idx[j] + 1]
        end
        if !iszero(v)
            push!(nnz, i)
            data[i] = v
        end
    end

    return Sparsik(A.n, field, nnz, data)
end

function __apply_vector_par(A::CSR_Sparsik, x)
    field = base_ring(A)
    m, n = size(A)
    y = zeros(field, n)

    @threads for i in 1 : m
        for j in A.row_ptr[i] : A.row_ptr[i + 1] - 1
            y[i] += A.val[j] * x[A.col_idx[j] + 1]
        end
    end

    return from_dense(y, field)
end

#------------------------------------------------------------------------------

# The resulting vectors are dense almost surely
# ...
# And that is not true
function Base.prod(A::CSR_Sparsik, x)
    field = base_ring(A)
    m, n = size(A)
    y = zeros(field, n)

    for i in 1 : m
        for j in A.row_ptr[i - 1] : A.row_ptr[i] - 1
            y[i] += A.val[j] * x[A.col_idx[j] + 1]
        end
    end

    return from_dense(y, field)
end

# TODO
function Base.prod(A::CSR_Sparsik, B::CSR_Sparsik)
    field = base_ring(A)
    m, n = size(A)
    y = zeros(field, n)

    for i in 1 : m
        for j in A.row_ptr[i - 1] : A.row_ptr[i] - 1
            y[i] += A.val[j] * x[A.col_idx[j] + 1]
        end
    end

    return from_dense(y, field)
end

#------------------------------------------------------------------------------

Base.iszero(A::CSR_Sparsik{T}) where {T} = isempty(A.val)

function first_nonzero(A::CSR_Sparsik)
    if iszero(A)
        return -1
    end

    i = 0
    while A.row_ptr[i] == 0
        i += 1
    end

    return to_plain(size(A), i, A.col_idx[A.row_ptr[i - 1]] + 1)
end

#------------------------------------------------------------------------------

function Base.zero(A::CSR_Sparsik)
    row_ptr = OffsetArray(zeros(Int, size(A, 1) + 1), Origin(0))
    return CSR_Sparsik(
        size(A)...,
        base_ring(A),
        row_ptr,
        ZeroBased(Int),
        ZeroBased(elem_type(base_ring(A)))
    )
end

#------------------------------------------------------------------------------

function scale!(A::CSR_Sparsik, c)
    if iszero(c)
        return A = zero(A)
    end

    for i in LinearIndices(A.val)
        A.val[i] *= c
    end

    return A
end

#------------------------------------------------------------------------------

function Base.reduce(A::CSR_Sparsik{T}, B::CSR_Sparsik{T}, c) where {T}
    field = base_ring(A)
    m, n = size(A)

    val = ZeroBased(elem_type(field))
    col_idx = ZeroBased(Int)
    row_ptr = ZeroBased(Int)

    push!(row_ptr, 0)

    for ridx in LinearIndices(A.row_ptr)[1:end]
        i, end_i = A.row_ptr[ridx - 1], A.row_ptr[ridx] - 1
        j, end_j = B.row_ptr[ridx - 1], B.row_ptr[ridx] - 1
        added = 0

        while i <= end_i || j <= end_j
            x = zero(field)

            if i > end_i
                new_idx = B.col_idx[j]
                x = c * B.val[j]
                j += 1
            elseif j > end_j
                new_idx = A.col_idx[i]
                x = A.val[i]
                i += 1
            else
                new_idx = min(A.col_idx[i], B.col_idx[j])
                if A.col_idx[i] > B.col_idx[j]
                    x = c * B.val[j]
                    j += 1
                elseif A.col_idx[i] < B.col_idx[j]
                    x = A.val[i]
                    i += 1
                else
                    x = A.val[i] + c * B.val[j]
                    i += 1
                    j += 1
                end
            end

            if !iszero(x)
                push!(col_idx, new_idx)
                push!(val, x)
                added += 1
            end

        end

        push!(row_ptr, last(row_ptr) + added)

    end

    return CSR_Sparsik(m, n, field, row_ptr, col_idx, val)
end

#------------------------------------------------------------------------------

# Do we really need this..

#------------------------------------------------------------------------------

==(A::CSR_Sparsik{T}, B::CSR_Sparsik{T}) where {T} = (size(A) == size(B) &&
                    A.col_idx == B.col_idx;
                    A.val == B.val;)

#------------------------------------------------------------------------------

function to_dense(A::CSR_Sparsik)
    ans = zeros(base_ring(A), size(A)...)
    for ridx in LinearIndices(A.row_ptr)[1:end]
        for i in A.row_ptr[ridx - 1] : A.row_ptr[ridx] - 1
            ans[ridx, A.col_idx[i] + 1] = A.val[i]
        end
    end
    return ans
end