docs: def alignment_to_contacts(
      sam_merged,
      assembly,
      output_dir,
      output_file_network=DEFAULT_NETWORK_FILE_NAME,
      output_file_chunk_data=DEFAULT_CHUNK_DATA_FILE_NAME,
      parameters=DEFAULT_PARAMETERS,
  ):
      """Generates a network file (in edgelist form) from an
      alignment in sam or bam format. Contigs are virtually split into
      'chunks' of nearly fixed size (by default between 500 and 1000 bp)
      to reduce size bias. The chunks are the network nodes and the edges
      are the contact counts.
  
      The network is in a strict barebone form so that it can be reused and
      imported quickly into other applications etc. Verbose information about
      every single node in the network is written on a 'chunk data' file,
      by default called 'idx_contig_hit_size_cov.txt'
  
      Parameters
      ----------
      sam_merged : file, str or pathlib.Path
          The alignment file in SAM/BAM format to be processed.
      assembly : file, str or pathlib.Path
          The initial assembly acting as the alignment file's reference genome.
      output_dir : str or pathlib.Path
          The output directory to write the network and chunk data into.
      output_dir_file_network : str or pathlib.Path, optional
          The specific file name for the output network file. Default is
          network.txt
      output_file_chunk_data : str or pathlib.Path, optional
          The specific file name for the output chunk data file. Default is
          idx_contig_hit_size_cov.txt
      parameters : dict, optional
          A dictionary of parameters for converting the alignment file into a
          network. These are:
          -size_chunk_threshold: the size (in bp) under which chunks are
          discarded. Default is 500.
          -mapq_threshold: the mapping quality under which alignments are
          discarded. Default is 10.
          -chunk_size: the default chunk size (in bp) when applicable, save
          smaller contigs or tail-ends. Default is 1000.
          -read_size: the size of reads used for mapping. Default is 65.
          -self_contacts: whether to count alignments between a chunk and
          itself. Default is False.
          -normalized: whether to normalize contacts by their coverage.
          Default is False.
  
      Returns
      -------
      chunk_complete_data : dict
          A dictionary where the keys are chunks in (contig, position) form and
          the values are their id, name, total contact count, size and coverage.
      all_contacts : dict
          A counter dictionary where the keys are chunk pairs and the values are
          their contact count.
      """
  
      all_contacts = collections.Counter()
      all_chunks = collections.Counter()
  
      #   Initialize parameters
      chunk_size = int(parameters["chunk_size"])
      mapq_threshold = int(parameters["mapq_threshold"])
      size_chunk_threshold = int(parameters["size_chunk_threshold"])
      read_size = int(parameters["read_size"])
      self_contacts = parameters["self_contacts"]
      normalized = parameters["normalized"]
  
      logger.info("Establishing chunk list...")
      chunk_complete_data = dict()
  
      #   Get all information about all chunks from all contigs
      #   (this gets updated at the end)
      global_id = 1
      for record in SeqIO.parse(assembly, "fasta"):
          length = len(record.seq)
  
          n_chunks = length // chunk_size
          n_chunks += (length % chunk_size) >= size_chunk_threshold
  
          for i in range(n_chunks):
  
              if (i + 1) * chunk_size <= length:
                  size = chunk_size
              else:
                  size = length % chunk_size
  
              chunk_name = "{}_{}".format(record.id, i)
              chunk_complete_data[chunk_name] = {
                  "id": global_id,
                  "hit": 0,
                  "size": size,
                  "coverage": 0,
              }
              global_id += 1
  
      logger.info("Opening alignment files...")
  
      current_read = None
  
      # Read the BAM file to detect contacts.
      with pysam.AlignmentFile(sam_merged, "rb") as alignment_merged_handle:
  
          names = alignment_merged_handle.references
          lengths = alignment_merged_handle.lengths
          names_and_lengths = {
              name: length for name, length in zip(names, lengths)
          }
  
          logger.info("Reading contacts...")
  
          # Since the BAM file is supposed to be sorted and interleaved,
          # pairs should be always grouped with one below the other (the exact
          # order doesn't matter since the network is symmetric, so we simply
          # treat the first one as 'forward' and the second one as 'reverse')
  
          # We keep iterating until two consecutive reads have the same name,
          # discarding ones that don't.
  
          while "Reading forward and reverse alignments alternatively":
              try:
                  my_read = next(alignment_merged_handle)
                  if current_read is None:
                      # First read
                      current_read = my_read
                      continue
  
                  elif current_read.query_name != my_read.query_name:
  
                      # print("{}_{}".format(current_read, my_read))
                      current_read = my_read
                      continue
  
                  read_forward, read_reverse = current_read, my_read
  
              except StopIteration:
                  break
  
              # Get a bunch of info about the alignments to pass the tests below
              read_name_forward = read_forward.query_name
              read_name_reverse = read_reverse.query_name
  
              flag_forward, flag_reverse = read_forward.flag, read_reverse.flag
  
              try:
                  assert read_name_forward == read_name_reverse
              except AssertionError:
                  logger.error(
                      "Reads don't have the same name: " "%s and %s",
                      read_name_forward,
                      read_name_reverse,
                  )
                  raise
  
              # To check if a flag contains 4
              # (digit on the third position from the right in base 2),
              # 4 = unmapped in SAM spec
              def is_unmapped(flag):
                  return np.base_repr(flag, padding=3)[-3] == "1"
  
              if is_unmapped(flag_forward) or is_unmapped(flag_reverse):
                  # print("Detected unmapped read on one end, skipping")
                  continue
  
              contig_name_forward = read_forward.reference_name
              contig_name_reverse = read_reverse.reference_name
  
              len_contig_for = names_and_lengths[contig_name_forward]
              len_contig_rev = names_and_lengths[contig_name_reverse]
  
              position_forward = read_forward.reference_start
              position_reverse = read_reverse.reference_start
  
              mapq_forward = read_forward.mapping_quality
              mapq_reverse = read_reverse.mapping_quality
  
              # Some more tests: checking for size, map quality, map status etc.
              mapq_test = min(mapq_forward, mapq_reverse) > mapq_threshold
  
              min_length = min(len_contig_for, len_contig_rev)
              length_test = min_length > size_chunk_threshold
  
              # Trickest test:
              #
              #
              #                contig
              #    pos1                          pos2
              #     ^                             ^
              # |-------|-------|-------|-------|---|
              # <-------><------><------><------><-->            <->
              #   chunk   chunk                  tail   size_chunk_threshold
              #
              # Test is passed if tail >= size_chunk_threshold (pos2)
              # or if the position is a non-tail chunk (pos1)
  
              if position_forward < chunk_size * (len_contig_for // chunk_size):
                  current_chunk_forward_size = chunk_size
              else:
                  current_chunk_forward_size = len_contig_for % chunk_size
  
              if position_reverse < chunk_size * (len_contig_rev // chunk_size):
                  current_chunk_reverse_size = chunk_size
              else:
                  current_chunk_reverse_size = len_contig_rev % chunk_size
  
              min_chunk_size = min(
                  current_chunk_forward_size, current_chunk_reverse_size
              )
  
              chunk_test = min_chunk_size >= size_chunk_threshold
  
              if mapq_test and length_test and chunk_test:
  
                  chunk_forward = position_forward // chunk_size
                  chunk_reverse = position_reverse // chunk_size
  
                  chunk_name_forward = "{}_{}".format(
                      contig_name_forward, chunk_forward
                  )
                  chunk_name_reverse = "{}_{}".format(
                      contig_name_reverse, chunk_reverse
                  )
  
                  if self_contacts or chunk_name_forward != chunk_name_reverse:
  
                      contact = tuple(
                          sorted((chunk_name_forward, chunk_name_reverse))
                      )
  
                      all_contacts[contact] += 1
  
                      chunk_key_forward = (
                          chunk_name_forward,
                          current_chunk_forward_size,
                      )
                      all_chunks[chunk_key_forward] += 1
  
                      chunk_key_reverse = (
                          chunk_name_reverse,
                          current_chunk_reverse_size,
                      )
                      all_chunks[chunk_key_reverse] += 1
  
      logger.info("Writing chunk data...")
  
      # Now we can update the chunk dictionary
      # with the info we gathered from the BAM file
  
      output_chunk_data_path = os.path.join(output_dir, output_file_chunk_data)
  
      with open(output_chunk_data_path, "w") as chunk_data_file_handle:
  
          for name in sorted(chunk_complete_data.keys()):
  
              chunk_data = chunk_complete_data[name]
              size = chunk_data["size"]
              chunk = (name, chunk_data["size"])
              hit = all_chunks[chunk]
              coverage = hit * read_size * 1.0 / size
              try:
                  chunk_complete_data[name]["hit"] = hit
                  chunk_complete_data[name]["coverage"] = coverage
              except KeyError:
                  logger.error(
                      "A mismatch was detected between the reference "
                      "genome and the genome used for the alignment "
                      "file, some sequence names were not found"
                  )
                  raise
  
              idx = chunk_complete_data[name]["id"]
              line = "{}\t{}\t{}\t{}\t{}\n".format(
                  idx, name, hit, size, coverage
              )
              chunk_data_file_handle.write(line)
  
      # Lastly, generate the network proper
  
      logger.info("Writing network...")
  
      output_network_path = os.path.join(output_dir, output_file_network)
  
      with open(output_network_path, "w") as network_file_handle:
  
          for chunks in sorted(all_contacts.keys()):
  
              chunk_name1, chunk_name2 = chunks
              contact_count = all_contacts[chunks]
  
              if normalized:
                  coverage1 = chunk_complete_data[chunk_name1]["coverage"]
                  coverage2 = chunk_complete_data[chunk_name2]["coverage"]
                  mean_coverage = np.sqrt(coverage1 * coverage2)
                  effective_count = contact_count * 1.0 / mean_coverage
              else:
                  effective_count = contact_count
  
              try:
                  idx1 = chunk_complete_data[chunk_name1]["id"]
                  idx2 = chunk_complete_data[chunk_name2]["id"]
                  line = "{}\t{}\t{}\n".format(idx1, idx2, effective_count)
                  network_file_handle.write(line)
              except KeyError as e:
                  logger.warning("Mismatch detected: %s", e)
  
      return chunk_complete_data, all_contacts