code: 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)
docs: direct (nonsymmetric) second moment matrix. Decide if we need dense, sparse, const