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

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]