C++ SpectrumPtr::getMZIntensityPairs方法代码示例

void verifySpectrumSize(const SpectrumList& sl, size_t index, size_t size)
{
    SpectrumPtr spectrum = sl.spectrum(index, true);
    vector<MZIntensityPair> data;
    spectrum->getMZIntensityPairs(data);
    unit_assert(data.size() == size);
}

PWIZ_API_DECL
msdata::SpectrumPtr SpectrumList_MZWindow::spectrum(size_t index, bool getBinaryData) const
{
    SpectrumPtr spectrum = inner_->spectrum(index, getBinaryData);
    if (!getBinaryData) return spectrum;

    vector<MZIntensityPair> data;
    spectrum->getMZIntensityPairs(data);

    vector<MZIntensityPair>::const_iterator begin = lower_bound(
        data.begin(), data.end(), MZIntensityPair(mzLow_,0), hasLowerMZ);

    vector<MZIntensityPair>::const_iterator end = upper_bound(
        data.begin(), data.end(), MZIntensityPair(mzHigh_,0), hasLowerMZ);

    vector<MZIntensityPair> newData;
    copy(begin, end, back_inserter(newData));

    BinaryDataArrayPtr intensityArray = spectrum->getIntensityArray();
    CVID intensityUnits = intensityArray.get() ? intensityArray->cvParam(MS_intensity_unit).cvid : CVID_Unknown;
    
    SpectrumPtr newSpectrum(new Spectrum(*spectrum));
    newSpectrum->binaryDataArrayPtrs.clear();
    newSpectrum->setMZIntensityPairs(newData, intensityUnits);

    return newSpectrum;
}

void cat(const char* filename)
{
    // open the data file

    ExtendedReaderList readers; // for vendor Reader support
    MSDataFile msd(filename, &readers);

    // verify that we have a SpectrumList

    if (!msd.run.spectrumListPtr.get())
        throw runtime_error("[mscat] No spectra found.");

    SpectrumList& spectrumList = *msd.run.spectrumListPtr;

    // write header

    const size_t columnWidth = 14;

    cout << "#"
         << setw(columnWidth) << "scanNumber"
         << setw(columnWidth) << "msLevel"
         << setw(columnWidth) << "m/z"
         << setw(columnWidth) << "intensity" << endl;

    // iterate through the spectra in the SpectrumList
    for (size_t i=0, size=spectrumList.size(); i!=size; i++)
    {
        // retrieve the spectrum, with binary data
        const bool getBinaryData = true;
        SpectrumPtr spectrum = spectrumList.spectrum(i, getBinaryData);

        // fill in MZIntensityPair vector for convenient access to binary data
        vector<MZIntensityPair> pairs;
        spectrum->getMZIntensityPairs(pairs);

        // iterate through the m/z-intensity pairs
        for (vector<MZIntensityPair>::const_iterator it=pairs.begin(), end=pairs.end(); it!=end; ++it)
        {
            cout << " "
                 << setw(columnWidth) << spectrum->nativeID
                 << setw(columnWidth) << "ms" + spectrum->cvParam(MS_ms_level).value
                 << setw(columnWidth) << fixed << setprecision(4) << it->mz
                 << setw(columnWidth) << fixed << setprecision(2) << it->intensity << endl;
        }
    }
}

void RAMPAdapter::Impl::getScanPeaks(size_t index, std::vector<double>& result) const
{
    // use previous spectrum if possible (it must have binary data)
    if (!lastSpectrum.get() || lastSpectrum->index != index) {
        lastSpectrum = msd_.run.spectrumListPtr->spectrum(index, true); // full read
    } else if (!lastSpectrum->hasBinaryData()) {
        // copy lastSpectrum header, avoids reread of header if format supports it
        lastSpectrum = msd_.run.spectrumListPtr->spectrum(lastSpectrum, true);
    }

    SpectrumPtr spectrum = lastSpectrum;
    result.clear();
    result.resize(spectrum->defaultArrayLength * 2);
    if (spectrum->defaultArrayLength == 0) return;

    spectrum->getMZIntensityPairs(reinterpret_cast<MZIntensityPair*>(&result[0]), 
                                  spectrum->defaultArrayLength);
}

void printSpectrumList(const SpectrumList& sl, ostream& os)
{
    os << "size: " << sl.size() << endl;

    for (size_t i=0, end=sl.size(); i<end; i++)
    {
        SpectrumPtr spectrum = sl.spectrum(i, false);
        vector<MZIntensityPair> data;
        spectrum->getMZIntensityPairs(data);

        os << spectrum->index << " " 
           << spectrum->id << ": ";

        copy(data.begin(), data.end(), ostream_iterator<MZIntensityPair>(os, " "));

        os << endl;
    }
}

PWIZ_API_DECL void ThresholdFilter::operator () (const SpectrumPtr s) const
{
    if (!msLevelsToThreshold.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return;

    // do nothing to empty spectra
    if (s->defaultArrayLength == 0)
        return;

    if (byType == ThresholdingBy_Count ||
        byType == ThresholdingBy_CountAfterTies)
    {
        // if count threshold is greater than number of data points, return as is
        if (s->defaultArrayLength <= threshold)
            return;
        else if (threshold == 0)
        {
            s->getMZArray()->data.clear();
            s->getIntensityArray()->data.clear();
            s->defaultArrayLength = 0;
            return;
        }
    }

    vector<MZIntensityPair> mzIntensityPairs;
    s->getMZIntensityPairs(mzIntensityPairs);

    if (orientation == Orientation_MostIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationMore_Predicate);
    else if (orientation == Orientation_LeastIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationLess_Predicate);
    else
        throw runtime_error("[threshold()] invalid orientation type");

    double tic = accumulate(mzIntensityPairs.begin(), mzIntensityPairs.end(), 0.0, MZIntensityPairIntensitySum());

    if (tic == 0)
    {
        s->getMZArray()->data.clear();
        s->getIntensityArray()->data.clear();
        s->defaultArrayLength = 0;
        return;
    }

    double bpi = orientation == Orientation_MostIntense ? mzIntensityPairs.front().intensity
                                                        : mzIntensityPairs.back().intensity;

    // after the threshold is applied, thresholdItr should be set to the first data point to erase
    vector<MZIntensityPair>::iterator thresholdItr;

    switch (byType)
    {
        case ThresholdingBy_Count:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + (size_t) threshold;

            // iterate backward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (thresholdItr == mzIntensityPairs.begin())
                    break;
                else if (i != (--thresholdItr)->intensity)
                {
                    ++thresholdItr;
                    break;
                }
            }
            break;

        case ThresholdingBy_CountAfterTies:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + ((size_t) threshold)-1;

            // iterate forward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (++thresholdItr == mzIntensityPairs.end() ||
                    i != thresholdItr->intensity)
                    break;
            }
            break;

        case ThresholdingBy_AbsoluteIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationMore_Predicate);
            else
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationLess_Predicate);
            break;

        case ThresholdingBy_FractionOfBasePeakIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    if (x_zero_p > mzRightLimit && makeBounds)
                        break;
                    intData[index_right] += y_val;
                    index_right++;
                    x_zero_p = binnedSpectrum.mzData[index_right];//binnedSpectrum.getPairFromIndex(index_right).first;//
                    y_val =  mzMath::y(x_zero_p, x_zero, y_max, sigma_squared);
                    //counter++;
                }

                //do not padd with zero if bin already exists
                //actually it is never supposed to happen since mz is increasing
                //so automatically create an entry and set it to 0
                //in some case (i suppose...) a previous peak could be larger
                //that is a reason to not reset the intensity to 0 too
                //if (! binnedSpectrum.getIntensityFromIndex(index_right))
                //binnedSpectrum.getOrInitIntensityFromIndex(index_right);
                /*} else {
                    binnedSpectrum.getOrInitIntensityFromIndex(index + 1);
                }*/
            }

            //pair<vector<double>, vector<double> > data = binnedSpectrum.getData();
            //printf("%d, %d\n", intData.size(), binnedSpectrum.mzData.size());
            ptr->setMZIntensityArrays(binnedSpectrum.mzData, intData,  CVID_Unknown);

            //here binning
            /*---comparing data intensities*/
            //printf("binnedSpectrumSize:%d, rawSize:%d\n", data.first.size(), rawSpectrum->getMZArray()->data.size());

            //bin
            //int nbBins = (int) (1700.0 / 0.1);
            map<int, vector<pair<double, double> > > rawIntensitiesSumByBinIndex, mzdbIntensitiesSumByIndex;
            vector<MZIntensityPair> rawMzIntensities;
            rawSpectrum->getMZIntensityPairs(rawMzIntensities);
            for (size_t j = 0; j < rawMzIntensities.size(); ++j) {
                const MZIntensityPair& pair = rawMzIntensities[j];
                int idx = (int) (pair.mz / 1);
                rawIntensitiesSumByBinIndex[idx].push_back(make_pair(pair.mz, pair.intensity));
            }

            for (size_t j= 0; j < intData.size(); ++j) {
                const double& mz = binnedSpectrum.mzData[j];
                const double& intensity = intData[j];
                int idx = (int) (mz / 1);
                mzdbIntensitiesSumByIndex[idx].push_back(make_pair(mz, intensity));
            }

            //iterate on raw map
            double counter = 0;
            double sum = 0;
            for (auto it = rawIntensitiesSumByBinIndex.begin(); it != rawIntensitiesSumByBinIndex.end(); ++it) {
                const int& idx = it->first;
                const vector<pair<double, double> >& vec = it->second;
                double integratedIntensity = 0;
                if (! vec.empty()) {
                    integrate(vec, integratedIntensity);
                }
                if (integratedIntensity) {
                    if (mzdbIntensitiesSumByIndex.find(idx) != mzdbIntensitiesSumByIndex.end()) {
                        const vector<pair<double, double> >& mzdbPoints = mzdbIntensitiesSumByIndex[idx];
                        //printf("v1 %f, %f\n", intensity, mzdbIntensity);
                        double mzDBIntegratedIntensity = 0;
                        integrate(mzdbPoints, mzDBIntegratedIntensity);
                        double error = (integratedIntensity - mzDBIntegratedIntensity) / integratedIntensity;
                        sum += error;
                        counter ++;

void test(bool indexed)
{
    if (os_) *os_ << "test(): indexed=\"" << boolalpha << indexed << "\"\n";

    MSData tiny;
    examples::initializeTiny(tiny);

    Serializer_mzML::Config config;
    config.indexed = indexed;
    Serializer_mzML serializer(config);  

    ostringstream oss;
    serializer.write(oss, tiny);

    if (os_) *os_ << "oss:\n" << oss.str() << endl;

    shared_ptr<istream> is(new istringstream(oss.str()));

    // dummy would normally be read in from file
  
    MSData dummy;

    ParamGroupPtr pg1(new ParamGroup);
    pg1->id = "CommonMS1SpectrumParams";
    pg1->cvParams.push_back(MS_positive_scan);
    pg1->cvParams.push_back(MS_full_scan);
    dummy.paramGroupPtrs.push_back(pg1);

    ParamGroupPtr pg2(new ParamGroup);
    pg2->id = "CommonMS2SpectrumParams";
    pg2->cvParams.push_back(MS_positive_scan);
    pg2->cvParams.push_back(MS_full_scan);
    dummy.paramGroupPtrs.push_back(pg2);

    // so we don't have any dangling references
    dummy.instrumentConfigurationPtrs.push_back(InstrumentConfigurationPtr(new InstrumentConfiguration("LCQDeca")));
    dummy.dataProcessingPtrs.push_back(DataProcessingPtr(new DataProcessing("XcaliburProcessing")));

    SpectrumListPtr sl = SpectrumList_mzML::create(is, dummy, indexed);

    // check easy functions

    unit_assert(sl.get());
    unit_assert(sl->size() == 4);
    unit_assert(sl->find ("S19") == 0);
    unit_assert(sl->findNative("19") == 0);
    unit_assert(sl->find("S20") == 1);
    unit_assert(sl->findNative("20") == 1);
    unit_assert(sl->find("S21") == 2);
    unit_assert(sl->findNative("21") == 2);
    unit_assert(sl->find("S22") == 3);
    unit_assert(sl->findNative("22") == 3);

    unit_assert(sl->findSpotID("A1").empty());
    IndexList spotIndexList = sl->findSpotID("A1,42x42,4242x4242");
    unit_assert(spotIndexList.size() == 1);
    unit_assert(spotIndexList[0] == 3);


    // check scan 19

    SpectrumPtr s = sl->spectrum(0); // read without binary data
    unit_assert(s.get());
    unit_assert(s->id == "S19");
    unit_assert(s->nativeID == "19");
    unit_assert(s->spotID.empty());
    unit_assert(s->cvParam(MS_ms_level).valueAs<int>() == 1);
    unit_assert(s->binaryDataArrayPtrs.empty());

    unit_assert(sl->spectrumIdentity(0).index == 0);
    unit_assert(sl->spectrumIdentity(0).id == "S19");
    unit_assert(sl->spectrumIdentity(0).nativeID == "19");
    unit_assert(sl->spectrumIdentity(0).spotID.empty());
 
    SpectrumPtr s_cache = sl->spectrum(0); // cache read
    unit_assert(s_cache.get() == s.get());

    s = sl->spectrum(0, true); // read with binary data
    unit_assert(s_cache.get() != s.get());

    vector<MZIntensityPair> pairs;
    s->getMZIntensityPairs(pairs);
    unit_assert(pairs.size() == 15);
    for (int i=0; i<15; i++)
        unit_assert(pairs[i].mz==i && pairs[i].intensity==15-i);

    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.size() == 1);
    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.back()->id == "CommonMS1SpectrumParams");
    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.back()->cvParams.size() == 2);

    // check scan 20

    s = sl->spectrum(1, true);
    unit_assert(s.get());
    unit_assert(s->id == "S20");
    unit_assert(s->nativeID == "20");
    unit_assert(s->spotID.empty());
    unit_assert(s->cvParam(MS_ms_level).valueAs<int>() == 2);

    unit_assert(sl->spectrumIdentity(1).index == 1);
//.........这里部分代码省略.........