<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _M2_const(Xvar, mask_X, xvarsum, xconst, Yvar, mask_Y, yvarsum, yconst, weights=None):
 r""" Computes the unnormalized covariance matrix between X and Y, exploiting constant input columns Computes the unnormalized covariance matrix :math:`C = X^\top Y` (for symmetric=False) or :math:`C = \frac{1}{2} (X^\top Y + Y^\top X)` (for symmetric=True). Suppose the data matrices can be column-permuted to have the form .. math: X &=& (X_{\mathrm{var}}, X_{\mathrm{const}}) Y &=& (Y_{\mathrm{var}}, Y_{\mathrm{const}}) with rows: .. math: x_t &=& (x_{\mathrm{var},t}, x_{\mathrm{const}}) y_t &=& (y_{\mathrm{var},t}, y_{\mathrm{const}}) where :math:`x_{\mathrm{const}},\:y_{\mathrm{const}}` are constant vectors. The resulting matrix has the general form: .. math: C &=& [X_{\mathrm{var}}^\top Y_{\mathrm{var}} x_{sum} y_{\mathrm{const}}^\top ] & & [x_{\mathrm{const}}^\top y_{sum}^\top x_{sum} x_{sum}^\top ] where :math:`x_{sum} = \sum_t x_{\mathrm{var},t}` and :math:`y_{sum} = \sum_t y_{\mathrm{var},t}`. Parameters Xvar : ndarray (T, m) Part of the data matrix X with :math:`m \le M` variable columns. mask_X : ndarray (M) Boolean array of size M of the full columns. False for constant column, True for variable column in X. xvarsum : ndarray (m) Column sum of variable part of data matrix X xconst : ndarray (M-m) Values of the constant part of data matrix X Yvar : ndarray (T, n) Part of the data matrix Y with :math:`n \le N` variable columns. mask_Y : ndarray (N) Boolean array of size N of the full columns. False for constant column, True for variable column in Y. yvarsum : ndarray (n) Column sum of variable part of data matrix Y yconst : ndarray (N-n) Values of the constant part of data matrix Y weights : None or ndarray (N) weights for all time steps. Returns ------- C : ndarray (M, N) Unnormalized covariance matrix. """

    C = np.zeros((len(mask_X), len(mask_Y)))
    # Block 11
    C[np.ix_(mask_X, mask_Y)] = _M2_dense(Xvar, Yvar, weights=weights)
    # other blocks
    xsum_is_0 = _is_zero(xvarsum)
    ysum_is_0 = _is_zero(yvarsum)
    xconst_is_0 = _is_zero(xconst)
    yconst_is_0 = _is_zero(yconst)
    # TODO: maybe we don't need the checking here, if we do the decision in the higher-level function M2
    # TODO: if not zero, we could still exploit the zeros in const and compute (and write!) this outer product
    # TODO: only to a sub-matrix
    # Block 12 and 21
    if weights is not None:
        wsum = np.sum(weights)
        xvarsum = np.sum(weights[:, None] * Xvar, axis=0)
        yvarsum = np.sum(weights[:, None] * Yvar, axis=0)
    else:
        wsum = Xvar.shape[0]
    if not (xsum_is_0 or yconst_is_0) or not (ysum_is_0 or xconst_is_0):
        C[np.ix_(mask_X, ~mask_Y)] = np.outer(xvarsum, yconst)
        C[np.ix_(~mask_X, mask_Y)] = np.outer(xconst, yvarsum)
    # Block 22
    if not (xconst_is_0 or yconst_is_0):
        C[np.ix_(~mask_X, ~mask_Y)] = np.outer(wsum*xconst, yconst)
    return C