def distance_diagonal_law(matrix, positions=None):
    n = min(matrix.shape)
    if positions is None:
        return np.array([np.average(np.diagonal(matrix, j)) for j in range(n)])
    else:
        contigs = positions_to_contigs(positions)
    def is_intra(i, j):
        return contigs[i] == contigs[j]
    max_intra_distance = max((len(contigs == u) for u in set(contigs)))
    intra_contacts = []
    inter_contacts = [np.average(np.diagonal(matrix, j))
                      for j in range(max_intra_distance, n)]
    for j in range(max_intra_distance):
        D = np.diagonal(matrix, j)
        for i in range(len(D)):
            diagonal_intra = []
            if is_intra(i, j):
                diagonal_intra.append(D[i])
        intra_contacts.append(np.average(np.array(diagonal_intra)))
    intra_contacts.extend(inter_contacts)
    return [positions, np.array(intra_contacts)]