Sum over each block of given subsampling factor, returns a matrix whose
    dimensions are this much as small (e.g. a 27x27 matrix binned with a
    subsampling factor equal to 3 will return a 9x9 matrix whose each component
    is the sum of the corresponding 3x3 block in the original matrix).
    Remaining columns and rows are summed likewise and added to the end of the
    new matrix.

    :note: Will not work for numpy verisons below 1.7

def bin_dense(M, subsampling_factor=3):
    """Sum over each block of given subsampling factor, returns a matrix whose
    dimensions are this much as small (e.g. a 27x27 matrix binned with a
    subsampling factor equal to 3 will return a 9x9 matrix whose each component
    is the sum of the corresponding 3x3 block in the original matrix).
    Remaining columns and rows are summed likewise and added to the end of the
    new matrix.

    :note: Will not work for numpy verisons below 1.7
    """

    m = min(M.shape)
    n = (m // subsampling_factor) * subsampling_factor

    if n == 0:
        return np.array([M.sum()])

    N = np.array(M[:n, :n], dtype=np.float64)
    N = N.reshape(n // subsampling_factor, subsampling_factor,
                  n // subsampling_factor, subsampling_factor).sum(axis=(1, 3))
    if m > n:
        remaining_row = M[n:, :n]
        remaining_col = M[:n, n:]
        remaining_square = M[n:m, n:m]
        R = remaining_row.reshape(m % subsampling_factor,
                                  m // subsampling_factor,
                                  subsampling_factor).sum(axis=(0, 2))
        C = remaining_col.T.reshape(m % subsampling_factor,
                                    m // subsampling_factor,
                                    subsampling_factor).sum(axis=(0, 2)).T
        S = remaining_square.sum()
        N = np.append(N, [R], axis=0)
        result = np.append(N, np.array([list(C) + [S]]).T, axis=1)
    else:
        result = N

    return result