docs: def matchToFeatures(fiContainer, specContainer, specfiles=None, fMassKey='mz',
                      sMassKey='obsMz', isotopeErrorList=(0),
                      precursorTolerance=5, toleranceUnit='ppm',
                      rtExpansionUp=0.10, rtExpansionDown=0.05, matchCharge=True,
                      scoreKey='pep', largerBetter=False):
      """Annotate :class:`Fi <maspy.core.Fi>` (Feature items) by matching
      :class:`Si <maspy.core.Si>` (Spectrum items) or :class:`Sii
      <maspy.core.Sii>` (Spectrum identification items).
  
      :param fiContainer: :class:`maspy.core.FeatureContainer`, contains ``Fi``.
      :param specContainer: :class:`maspy.core.MsrunContainer` or
          :class:`maspy.core.SiiContainer`, contains ``Si`` or ``Sii``.
      :param specfiles: filenames of ms-run files, if specified consider only
          items from those files
      :type specfiles: str, list or None
      :param fMassKey: mass attribute key in :attr:`Fi.__dict__`
      :param sMassKey: mass attribute key in :attr:`Si.__dict__` or
          :attr:`Sii.__dict__` (eg 'obsMz', 'excMz')
      :param isotopeErrorList: allowed isotope errors relative to the spectrum
          mass, for example "0" or "1". If no feature has been matched with
          isotope error 0, the spectrum mass is increased by the mass difference
          of carbon isotopes 12 and 13 and matched again. The different isotope
          error values are tested in the specified order therefore "0" should
          normally be the first value of the list.
      :type isotopeErrorList: list or tuple of int
      :param precursorTolerance: the largest allowed mass deviation of ``Si`` or
          ``Sii`` relative to ``Fi``
      :param toleranceUnit: defines how the ``precursorTolerance`` is applied to
          the mass value of ``Fi``. ``"ppm": mass * (1 +/- tolerance*1E-6)`` or
          ``"da": mass +/- value``
      :param rtExpansionUp: relative upper expansion of ``Fi`` retention time
          area. ``limitHigh = Fi.rtHigh + (Fi.rtHigh - Fi.rtLow) * rtExpansionUp``
      :param rtExpansionDown: relative lower expansion of ``Fi`` retention time
          area. ``limitLow = Fi.rtLow - (Fi.rtHigh - Fi.rtLow) * rtExpansionDown``
      :param matchCharge: bool, True if ``Fi`` and ``Si`` or ``Sii`` must have the
          same ``charge`` state to be matched.
      :param scoreKey: ``Sii`` attribute name used for scoring the identification
          reliability
      :param largerBetter: bool, True if higher score value means a better
          identification reliability
  
      .. note:
          Concerning the feature retention area expansion. If ``Si`` or ``Sii`` is
          matched to multiple ``Fi`` the rt expansion is removed and the matching
          is repeated.
  
      .. note:
          If the ``specContainer`` is a ``SiiContainer`` then matched ``Fi`` are
          annotated with :attr:`Sii.peptide`, if multiple ``Sii`` are matched to
          ``Fi`` the one with the best score is used.
  
      #TODO: this function is nested pretty badly and should maybe be rewritten
      #TODO: replace tolerance unit "ppm" by tolerance mode "relative" and change
          repsective calculations
      """
      isotopeErrorList = aux.toList(isotopeErrorList)
  
      if specContainer.__class__.__name__ == 'MsrunContainer':
          listKeySpecIds = 'siIds'
      else:
          listKeySpecIds = 'siiIds'
      specContainerSpecfiles = [_ for _ in viewkeys(specContainer.info)]
  
      if specfiles is not None:
          specfiles = aux.toList(specfiles)
      else:
          specfiles = [_ for _ in viewkeys(fiContainer.info)]
      specfiles = list(set(specfiles).intersection(set(specContainerSpecfiles)))
  
      for specfile in specfiles:
          multiMatchCounter = int()
          isotopeErrorMatchCounter = int()
          specArrays = specContainer.getArrays([sMassKey, 'rt', 'charge',
                                                'msLevel'], specfiles=specfile
                                                )
          featureArrays = fiContainer.getArrays(['rtHigh', 'rtLow', 'charge',
                                                 fMassKey], specfiles=specfile,
                                                 sort=fMassKey
                                                )
          featureArrays['rtHighExpanded'] = (featureArrays['rtHigh'] +
                                             (featureArrays['rtHigh'] -
                                              featureArrays['rtLow']) *
                                             rtExpansionUp
                                             )
          featureArrays['rtLowExpanded'] = (featureArrays['rtLow'] -
                                            (featureArrays['rtHigh'] -
                                             featureArrays['rtLow']) *
                                            rtExpansionDown
                                            )
  
          specFeatureDict = dict() ## key = scanNr, value = set(featureKeys)
          featureSpecDict = dict() ## key = featureKey, value = set(scanNrs)
  
          for specPos, specId in enumerate(specArrays['id']):
              specZ = specArrays['charge'][specPos]
              if specZ is None:
                  continue
              specMass = specArrays[sMassKey][specPos]
              specRt = specArrays['rt'][specPos]
  
              matchComplete = False
              isotopeErrorPos = 0
  
              while not matchComplete:
                  isotopeError = isotopeErrorList[isotopeErrorPos]
  
                  # calculate mass limits for each isotope error
                  if toleranceUnit.lower() == 'ppm':
                      specMassHigh = ((specMass + isotopeError * 1.003355 / specZ)
                                      * (1 + precursorTolerance*1E-6)
                                      )
                      specMassLow = ((specMass + isotopeError * 1.003355 / specZ)
                                     * (1 - precursorTolerance*1E-6)
                                     )
                  elif toleranceUnit.lower() == 'da':
                      specMassHigh = ((specMass + isotopeError * 1.003355 / specZ)
                                      + precursorTolerance
                                      )
                      specMassLow  = ((specMass + isotopeError * 1.003355 / specZ)
                                      - precursorTolerance
                                      )
  
                  posL = bisect.bisect_left(featureArrays[fMassKey],
                                            specMassLow
                                            )
                  posR = bisect.bisect_right(featureArrays[fMassKey],
                                             specMassHigh
                                             )
  
                  if matchCharge:
                      chargeMask = (featureArrays['charge'][posL:posR] == specZ)
  
                  fRtHighKey = 'rtHighExpanded'
                  fRtLowKey = 'rtLowExpanded'
                  for fRtHighKey, fRtLowKey in [('rtHighExpanded',
                                                 'rtLowExpanded'),
                                                ('rtHigh', 'rtLow')
                                                ]:
                      rtMask = ((featureArrays[fRtLowKey][posL:posR] <= specRt) &
                                (featureArrays[fRtHighKey][posL:posR] >= specRt)
                                )
                      if matchCharge:
                          matchedFeatureIds = featureArrays['id'][posL:posR][rtMask & chargeMask]
                      else:
                          matchedFeatureIds = featureArrays['id'][posL:posR][rtMask]
  
                      if len(matchedFeatureIds) <= 1:
                          break
  
                  # if exactly one feature has been matched,
                  if len(matchedFeatureIds) > 0:
                      if len(matchedFeatureIds) == 1:
                          matchComplete = True
                          if isotopeErrorList[isotopeErrorPos] != 0:
                              isotopeErrorMatchCounter += 1
                      else:
                          #Stop if Spectrum can be matched to multiple features
                          multiMatchCounter += 1
                          break
  
                  isotopeErrorPos += 1
                  if isotopeErrorPos >= len(isotopeErrorList):
                      #Stop if all allowed isotope errors have been tested
                      break
  
              if matchComplete:
                  for featureId in matchedFeatureIds:
                      getattr(fiContainer.container[specfile][featureId],
                              listKeySpecIds
                              ).append(specId)
                      fiContainer.container[specfile][featureId].isMatched = True
                      specFeatureDict[specId] = featureId
                      featureSpecDict[featureId] = specId
  
          stats = dict()
          stats['totalFeatures'] = len(featureArrays['id'])
          stats['matchedFeatures'] = len(featureSpecDict)
          stats['relMatchedFeatures'] = round(100*stats['matchedFeatures']/stats['totalFeatures'], 1)
          stats['totalSpectra'] = len(specArrays['id'][(specArrays['msLevel'] != 1)])
          stats['matchedSpectra'] = len(specFeatureDict)
          stats['relMatchedSpectra'] = round(100*stats['matchedSpectra']/stats['totalSpectra'], 1)
  
          print('------', specfile, '------')
          print('Annotated features:\t\t\t', stats['matchedFeatures'], '/', stats['totalFeatures'], '=', stats['relMatchedFeatures'], '%')
          print('Spectra matched to features:\t\t', stats['matchedSpectra'], '/', stats['totalSpectra'], '=', stats['relMatchedSpectra'], '%')
          if multiMatchCounter != 0:
                  print('Discarded because of multiple matches:\t', multiMatchCounter)
          if isotopeErrorMatchCounter != 0:
                  print('Isotope error matched spectra:\t\t', isotopeErrorMatchCounter)
  
          #annotate feature with sii information (peptide, sequence, score)
          if isinstance(specContainer, maspy.core.SiiContainer):
              for featureId in viewkeys(featureSpecDict):
                  matches = list()
                  for specId in fiContainer.container[specfile][featureId].siiIds:
                      _sii = specContainer.getValidItem(specfile, specId)
                      score = getattr(_sii, scoreKey)
                      peptide = _sii.peptide
                      sequence = _sii.sequence
                      matches.append([score, peptide, sequence])
                  matches.sort(reverse=largerBetter)
  
                  fiContainer.container[specfile][featureId].isAnnotated = True
                  fiContainer.container[specfile][featureId].score = matches[0][0]
                  fiContainer.container[specfile][featureId].peptide = matches[0][1]
                  fiContainer.container[specfile][featureId].sequence = matches[0][2]