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])
  #            else:
  #                diagonal_inter.append(D[i])
  #        inter_contacts.append(np.average(np.array(diagonal_inter)))
          intra_contacts.append(np.average(np.array(diagonal_intra)))
  
      intra_contacts.extend(inter_contacts)
  
      return [positions, np.array(intra_contacts)]
docs: Compute a distance law trend using the contact averages of equal distances.
      Specific positions can be supplied if needed.