x = x + beta * A(:,j), where x is a dense vector and A(:,j) is sparse.
494226
private static int scatter(SparseMatrix A, int j, double beta, int[] w, double[] x, int mark, SparseMatrix C, int nz) {
        int[] Ap = A.colIndex;
        int[] Ai = A.rowIndex;
        double[] Ax = A.x;

        int[] Ci = C.rowIndex;
        for (int p = Ap[j]; p < Ap[j + 1]; p++) {
            int i = Ai[p];                // A(i,j) is nonzero
            if (w[i] < mark) {
                w[i] = mark;              // i is new entry in column j
                Ci[nz++] = i;             // add i to pattern of C(:,j)
                x[i] = beta * Ax[p];      // x(i) = beta*A(i,j)
            } else {
                x[i] += beta * Ax[p];     // i exists in C(:,j) already
            }
        }

        return nz;
    }
1230368
public static void jacobian_Control3( DMatrixRMaj L_full ,
										  double beta[] ,  DMatrixRMaj A)
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2];

		final double ld[] = L_full.data;

		for( int i = 0; i < 3; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+3]*b1 +   ld[li+4]*b2;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+4]*b1 + 2*ld[li+5]*b2;
		}
	}
1352783
def _zero_many(self, i, j):
    
    i, j, M, N = self._prepare_indices(i, j)

    n_samples = len(i)
    offsets = np.empty(n_samples, dtype=self.indices.dtype)
    ret = _sparsetools.csr_sample_offsets(M, N, self.indptr, self.indices,
                                          n_samples, i, j, offsets)
    if ret == 1:
        # rinse and repeat
        self.sum_duplicates()
        _sparsetools.csr_sample_offsets(M, N, self.indptr,
                                        self.indices, n_samples, i, j,
                                        offsets)

    # only assign zeros to the existing sparsity structure
    self.data[list(offsets[offsets > -1])] = 0
243161
def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """RHS for condensed system"""
    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))
    dxnu = mysolve(LU, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
201410
public static void multAddOuter( double alpha , DMatrix4x4 A , double beta , DMatrix4 u , DMatrix4 v , DMatrix4x4 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
    }
713083
private void computeA( DenseMatrix64F A , DenseMatrix64F H , double lambda )
    {
        final int numParam = param.getNumElements();

        A.set(H);
        for( int i = 0; i < numParam; i++ ) {
            A.set(i,i, A.get(i,i) + lambda);
        }
    }
1257331
def _rand_sparse(m, n, density, format='csr'):
    """"""
    nnz = max(min(int(m*n*density), m*n), 0)

    row = np.random.randint(low=0, high=m-1, size=nnz)
    col = np.random.randint(low=0, high=n-1, size=nnz)
    data = np.ones(nnz, dtype=float)

    # duplicate (i,j) entries will be summed together
    return sp.sparse.csr_matrix((data, (row, col)), shape=(m, n))
668901
def _sparse_fit(self, X, strategy, missing_values, fill_value):
        """"""
        mask_data = _get_mask(X.data, missing_values)
        n_implicit_zeros = X.shape[0] - np.diff(X.indptr)

        statistics = np.empty(X.shape[1])

        if strategy == "constant":
            # for constant strategy, self.statistcs_ is used to store
            # fill_value in each column
            statistics.fill(fill_value)

        else:
            for i in range(X.shape[1]):
                column = X.data[X.indptr[i]:X.indptr[i + 1]]
                mask_column = mask_data[X.indptr[i]:X.indptr[i + 1]]
                column = column[~mask_column]

                # combine explicit and implicit zeros
                mask_zeros = _get_mask(column, 0)
                column = column[~mask_zeros]
                n_explicit_zeros = mask_zeros.sum()
                n_zeros = n_implicit_zeros[i] + n_explicit_zeros

                if strategy == "mean":
                    s = column.size + n_zeros
                    statistics[i] = np.nan if s == 0 else column.sum() / s

                elif strategy == "median":
                    statistics[i] = _get_median(column,
                                                n_zeros)

                elif strategy == "most_frequent":
                    statistics[i] = _most_frequent(column,
                                                   0,
                                                   n_zeros)
        return statistics
523198
public static void jacobian_Control4( DMatrixRMaj L_full ,
										  double beta[] , DMatrixRMaj A )
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2]; double b3 = beta[3];

		final double ld[] = L_full.data;

		for( int i = 0; i < 6; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2 +   ld[li+3]*b3;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+4]*b1 +   ld[li+5]*b2 +   ld[li+6]*b3;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+5]*b1 + 2*ld[li+7]*b2 +   ld[li+8]*b3;
			A.data[indexA++] =   ld[li+3]*b0 +   ld[li+6]*b1 +   ld[li+8]*b2 + 2*ld[li+9]*b3;
		}
	}
725689
private static void scatter(SparseVector v, double[] z) {
        int[] index = v.getIndex();
        int used = v.getUsed();
        double[] data = v.getData();
        Arrays.fill(z, 0);
        for (int i = 0; i < used; ++i)
            z[index[i]] = data[i];
    }
469883
def leastsq_generalized(x, y, x2, y2, res_std=None, res2_std=None, fit_offset=False, fit_offset2=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
        x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()

    n_obs, n_var = x.shape
    offset, offset_ss = np.zeros(n_var, dtype="float32"), np.zeros(n_var, dtype="float32")
    gamma = np.ones(n_var, dtype="float32")

    if (res_std is None) or (res2_std is None): res_std, res2_std = np.ones(n_var), np.ones(n_var)
    ones, zeros = np.ones(n_obs), np.zeros(n_obs)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        x, y = np.vstack((make_dense(x)/res_std, x2/res2_std)), np.vstack((make_dense(y)/res_std, y2/res2_std))

    if fit_offset and fit_offset2:
        for i in range(n_var):
            A = np.c_[np.vstack((np.c_[ones/res_std[i], zeros], np.c_[zeros, ones/res2_std[i]])), x[:, i]]
            offset[i], offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset:
        for i in range(n_var):
            A = np.c_[np.hstack((ones/res_std[i], zeros)), x[:, i]]
            offset[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset2:
        for i in range(n_var):
            A = np.c_[np.hstack((zeros, ones/res2_std[i])), x[:, i]]
            offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    else:
        for i in range(n_var):
            A = np.c_[x[:, i]]
            gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))

    offset[np.isnan(offset)] = 0
    offset_ss[np.isnan(offset_ss)] = 0
    gamma[np.isnan(gamma)] = 0

    return offset, offset_ss, gamma
393746
public double dot(SparseVector b) {
		SparseVector a = this;
		if (a.N != b.N) throw new IllegalArgumentException("Vector lengths disagree. " + a.N + " != " + b.N);
		double sum = 0.0;

		// iterate over the vector with the fewest nonzeros
		if (a.symbolTable.size() <= b.symbolTable.size()) {
			for (int i : a.symbolTable)
				if (b.symbolTable.contains(i)) sum += a.get(i) * b.get(i);
		}
		else  {
			for (int i : b.symbolTable)
				if (a.symbolTable.contains(i)) sum += a.get(i) * b.get(i);
		}
		return sum;
	}
621622
def _M2(Xvar, Yvar, mask_X=None, mask_Y=None, xsum=0, xconst=0, ysum=0, yconst=0, weights=None, diag_only=False):
    
    if mask_X is None and mask_Y is None:
        return _M2_dense(Xvar, Yvar, weights=weights, diag_only=diag_only)
    else:
        # Check if one of the masks is not None, modify it and also adjust the constant columns:
        if mask_X is None:
            mask_X = np.ones(Xvar.shape[1], dtype=np.bool)
            xconst = np.ones(0, dtype=float)
        if mask_Y is None:
            mask_Y = np.ones(Yvar.shape[1], dtype=np.bool)
            yconst = np.ones(0, dtype=float)
    if _is_zero(xsum) and _is_zero(ysum) or _is_zero(xconst) and _is_zero(yconst):
        return _M2_sparse(Xvar, mask_X, Yvar, mask_Y, weights=weights)
    else:
        return _M2_const(Xvar, mask_X, xsum[mask_X], xconst, Yvar, mask_Y, ysum[mask_Y], yconst, weights=weights)
932371
double refined_vector_value(SparseVector composite, SparseVector vec, int sign)
    {
        double sum = 0.0;
        for (Map.Entry<Integer, Double> entry : vec.entrySet())
        {
            sum += Math.pow(entry.getValue(), 2) + sign * 2 * composite.get(entry.getKey()) * entry.getValue();
        }
        return sum;
    }
1238289
def array2d(X, dtype=None, order=None, copy=False):
    """"""
    if sp.issparse(X):
        raise TypeError('A sparse matrix was passed, but dense data '
                        'is required. Use X.toarray() to convert to dense.')
    X_2d = np.asarray(np.atleast_2d(X), dtype=dtype, order=order)
    _assert_all_finite(X_2d)
    if X is X_2d and copy:
        X_2d = safe_copy(X_2d)
    return X_2d
552652
def MarginalBeta(self, i):
        
        alpha0 = self.params.sum()
        alpha = self.params[i]
        return Beta(alpha, alpha0 - alpha)
549333
private SparseDoubleVector updateRow(int row) {
        SparseDoubleVector rowVec = rowToColumns.get(row);
        if (rowVec == null) {
            rowVec = emptyVector.instanceCopy();
            rowToColumns.put(row, rowVec);
        }

        return rowVec;
    }