code: def bin_dense(M, subsampling_factor=3):
      
  
      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
docs: 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