From a symmetrical matrix M of size n, return a vector d whose each
    component d[i] is the total sum of a square of 2*window+1 size centered on
    the i-th main diagonal element. Edge elements may be extrapolated based on
    the square size reduction (i.e. for window = 4, the first component will be
    equal to the first diagonal pixel multiplied by 81, the second one will be
    equal to the first 2x2 square on the diagonal multiplied by 81/4, etc.),
    except in the case of circular genomes. If they aren't, d will be of size
    n - 2*(window-1) instead of n.

def domainogram(M, window=None, circ=False, extrapolate=True):
    """From a symmetrical matrix M of size n, return a vector d whose each
    component d[i] is the total sum of a square of 2*window+1 size centered on
    the i-th main diagonal element. Edge elements may be extrapolated based on
    the square size reduction (i.e. for window = 4, the first component will be
    equal to the first diagonal pixel multiplied by 81, the second one will be
    equal to the first 2x2 square on the diagonal multiplied by 81/4, etc.),
    except in the case of circular genomes. If they aren't, d will be of size
    n - 2*(window-1) instead of n.
    """

    # Sanity checks
    if not type(M) is np.ndarray:
        M = np.array(M)

    if M.shape[0] != M.shape[1]:
        raise ValueError("Matrix is not square.")

    try:
        n = min(M.shape)
    except AttributeError:
        n = M.size

    # Default window argument
    if window is None:
        window = max(n // 100, 5)

    if window >= n:
        raise ValueError("Please choose a smaller window size.")

    if circ:
        d = [(np.sum(M[-i + window:, -i + window:]) +
             np.sum(M[:i - window + 1, :i - window + 1])
             for i in range(window))]
    elif extrapolate:
        d = [(np.sum(M[0:2 * i + 1, 0:2 * i + 1]) *
              ((2 * window + 1) ** 2.0) / ((2 * i + 1)**2.0))
             for i in range(window)]
    else:
        d = []

    d += [np.sum(M[i - window:i + window + 1, i - window:i + window + 1])
          for i in range(window, n - window)]

    if circ:
        d += [M[i:, i:].sum() + M[:n - i, n - i].sum()
              for i in range(n - window, n)]
    elif extrapolate:
        d += [M[i - window:, i - window:].sum() * ((2 * window + 1)**2.0) /
              ((2 * (n - i) + 1)**2.0) for i in range(n - window, n)]

    return np.array(d)