code: 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)]