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

    void writeCompressedPeaks(SpectrumPtr s, ostream& os)
    {
        // Build arrays to hold peaks prior to compression
        int numPeaks = (int) s->defaultArrayLength;
        double *pD = new double[numPeaks];
        float *pF = new float[numPeaks];

        const BinaryDataArray& mzArray = *s->getMZArray();
        const BinaryDataArray& intensityArray = *s->getIntensityArray();
        for(int j = 0; j < numPeaks; j++)
        {
            pD[j] = mzArray.data[j];
            pF[j] = (float) intensityArray.data[j];
        }

        // compress mz
        uLong sizeM = (uLong) (numPeaks * sizeDoubleMSn);
        uLong comprLenM = compressBound(sizeM);
        Byte *comprM = (Byte*)calloc((uInt)comprLenM, 1);
        int retM = compress(comprM, &comprLenM, (const Bytef*)pD, sizeM);
        
        // compress intensity
        uLong sizeI = (uLong) (numPeaks * sizeFloatMSn);
        uLong comprLenI = compressBound(sizeI);
        Byte *comprI = (Byte*)calloc((uInt)comprLenI, 1);
        int retI = compress(comprI, &comprLenI, (const Bytef*)pF, sizeI);

        // Write the compressed peaks if all is well
        if ((Z_OK == retM) && (Z_OK == retI))
        {
            // write length of compressed array of m/z
            os.write(reinterpret_cast<char *>(&comprLenM), sizeIntMSn);

            // write length of compressed array of intensities
            os.write(reinterpret_cast<char *>(&comprLenI), sizeIntMSn);

            // write compressed array of m/z
            os.write(reinterpret_cast<char *>(comprM), comprLenM);

            // write compressed array of intensities
            os.write(reinterpret_cast<char *>(comprI), comprLenI);
        }

        // Clean up memory
        free(comprM);
        free(comprI);
        delete [] pD;
        delete [] pF;
        
        // In case of error, throw exception AFTER cleaning up memory
        if (Z_OK != retM || Z_OK != retI)
        {
            throw runtime_error("[Serializer_MSn::writeCompressedPeaks] Error compressing peaks.");
        }
    }

    void writeSpectrumText(SpectrumPtr s, ostream& os)
    {
        os << std::setprecision(7); // 124.4567
        
        // Write the scan numbers 
        os << "S\t";
        int scanNum = getScanNumber(s);
        os << scanNum <<  "\t" << scanNum << "\t";

        // Write the precursor mz
        Precursor& precur = s->precursors[0];
        SelectedIon& si = precur.selectedIons[0];
        double mz = si.cvParam(MS_selected_ion_m_z).valueAs<double>();
        os << mz << "\n";
        
        // Write the scan time, if available
        if( s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() )
          os << "I\tRTime\t" << s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() << "\n";

        // Collect charge and mass info
        vector<int> charges;
        vector<double> masses;
        int numChargeStates = 0;
        // for each selected ion
        BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons){
          numChargeStates += getChargeStates(curIon, charges, masses);
        }

        // Write EZ lines if accurate masses are available
        CVParam massParam = si.cvParam(MS_accurate_mass);
        if( !massParam.empty() ){
          for(int i=0; i < numChargeStates; i++){
            os << "I\tEZ\t" << charges[i] << "\t" << masses[i] << "\t0\t0" << endl; // pad last two fields with 0
          }
        }

        // For each charge, write the charge and mass
        for(int i = 0; i < numChargeStates; i++)
        {
            os << "Z\t" << charges[i] << "\t" << masses[i] << "\n"; 
        }

        
        // Write each mz, intensity pair
        const BinaryDataArray& mzArray = *s->getMZArray();
        const BinaryDataArray& intensityArray = *s->getIntensityArray();
        for (size_t p=0; p < s->defaultArrayLength; ++p)
        {
            os << mzArray.data[p] << " " << intensityArray.data[p] << "\n";
        }
    }

PWIZ_API_DECL SpectrumPtr SpectrumList_MetadataFixer::spectrum(size_t index, bool getBinaryData) const
{
    // always get binary data
    SpectrumPtr s = inner_->spectrum(index, true);

    BinaryDataArrayPtr mzArray = s->getMZArray();
    BinaryDataArrayPtr intensityArray = s->getIntensityArray();
    if (!mzArray.get() || !intensityArray.get())
        return s;

    vector<double>& mzs = mzArray->data;
    vector<double>& intensities = intensityArray->data;

    double tic = 0;
    if (!mzs.empty())
    {
        double bpmz, bpi = -1;
        for (size_t i=0, end=mzs.size(); i < end; ++i)
        {
            tic += intensities[i];
            if (bpi < intensities[i])
            {
                bpi = intensities[i];
                bpmz = mzs[i];
            }
        }

        replaceCvParam(*s, MS_base_peak_intensity, bpi, MS_number_of_counts);
        replaceCvParam(*s, MS_base_peak_m_z, bpmz, MS_m_z);
        replaceCvParam(*s, MS_lowest_observed_m_z, mzs.front(), MS_m_z);
        replaceCvParam(*s, MS_highest_observed_m_z, mzs.back(), MS_m_z);
    }

    replaceCvParam(*s, MS_TIC, tic, MS_number_of_counts);

    return s;
}

    void writeSpectrumBinary(SpectrumPtr s, int version, bool compress, ostream& os)
    {
      // todo
      // MSnScanInfo header(s);
      // header.writeSpectrumHeader(version, compress, os);

        int scanNum = getScanNumber(s);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn); // Yes, there are two

        Precursor& precur = s->precursors[0];
        SelectedIon& si = precur.selectedIons[0];
        double mz = si.cvParam(MS_selected_ion_m_z).valueAs<double>();
        os.write(reinterpret_cast<char *>(&mz), sizeDoubleMSn);

        float rt = (float) s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds();
        os.write(reinterpret_cast<char *>(&rt), sizeFloatMSn);
        
        if (version >= 2)
        {
            float basePeakIntensity = s->cvParam(MS_base_peak_intensity).valueAs<float>();
            os.write(reinterpret_cast<char *>(&basePeakIntensity), sizeFloatMSn);

            double basePeakMZ = s->cvParam(MS_base_peak_m_z).valueAs<double>();
            os.write(reinterpret_cast<char *>(&basePeakMZ), sizeDoubleMSn);
            
            // We don't have this information, but we need to write something,
            // so pad with 0's. (version 2 specific data)
            double conversionFactorA = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorA), sizeDoubleMSn); 
            double conversionFactorB = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorB), sizeDoubleMSn); 

            double tic = s->cvParam(MS_total_ion_current).valueAs<double>();
            os.write(reinterpret_cast<char *>(&tic), sizeDoubleMSn);

            // TODO
            float ionInjectionTime = (float)0;
            os.write(reinterpret_cast<char *>(&ionInjectionTime), sizeFloatMSn);
        }
        
        vector<int> charges;
        vector<double> masses;
        int numChargeStates = 0;
        BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons)
        {
            numChargeStates += getChargeStates(curIon, charges, masses);
        }
        os.write(reinterpret_cast<char *>(&numChargeStates), sizeIntMSn);
        
        bool hasAccurateMass = false;
        if (version == 3)
        {
          int numEzStates = 0;
          CVParam massParam = si.cvParam(MS_accurate_mass);
          if (!massParam.empty())
          {
            numEzStates = numChargeStates;
            hasAccurateMass = true;
          }
          os.write(reinterpret_cast<char *>(&numEzStates), sizeIntMSn);
        }

        int numPeaks = (int) s->defaultArrayLength;
        os.write(reinterpret_cast<char *>(&numPeaks), sizeIntMSn);

        // end spectrum header info

        // Write out each charge state and corresponding mass
        for(int i = 0; i < numChargeStates; i++)
        {
            os.write(reinterpret_cast<char *>(&(charges[i])), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&(masses[i])), sizeDoubleMSn);
        }
    
        // if there are accurate masses, write out EZ entries
        if( hasAccurateMass ){
          float blank = 0;  // we don't have rTime or area, pad with zeros
          for(int i=0; i < numChargeStates; i++){
            os.write(reinterpret_cast<char *>(&charges[i]), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&masses[i]), sizeDoubleMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
          }
        }

        // Do we need to write compressed m/z, intensity arrays?
        if (compress)
        {
            writeCompressedPeaks(s, os);
        }
        else
        {
            // No need to compress, just write out the arrays
            const BinaryDataArray& mzArray = *s->getMZArray();
            const BinaryDataArray& intensityArray = *s->getIntensityArray();
            for(int i = 0; i < numPeaks; i++)
            {
                double mzPeak = mzArray.data[i];
                os.write(reinterpret_cast<char *>(&mzPeak), sizeDoubleMSn);
//.........这里部分代码省略.........

int test()
{
    int failedTests = 0;

    // create the spectrum list
    SpectrumListSimple* sl = new SpectrumListSimple;
    SpectrumListPtr originalList(sl);

    for (size_t i=0; i < testScanSummerSize; ++i)
    {
        TestScanSummerCalculator& t = testScanSummerCalculators[i];
        SpectrumPtr s(new Spectrum);
        s->set(MS_MSn_spectrum);
        s->set(MS_ms_level,2);
        s->index = i;
        s->scanList.scans.push_back(Scan());
        Scan& scanRef = s->scanList.scans[0];
        scanRef.set(MS_scan_start_time,t.rTime,UO_second);
        s->precursors.push_back(Precursor(t.inputPrecursorMZ));
        
        vector<double> inputMZArray = parseDoubleArray(t.inputMZArray);
        vector<double> inputIntensityArray = parseDoubleArray(t.inputIntensityArray);
        s->setMZIntensityArrays(inputMZArray, inputIntensityArray, MS_number_of_detector_counts);
        s->defaultArrayLength = inputMZArray.size();

        scanRef.scanWindows.push_back(ScanWindow());
        scanRef.scanWindows[0].set(MS_scan_window_lower_limit,inputMZArray[0]);
        scanRef.scanWindows[0].set(MS_scan_window_upper_limit,inputMZArray[inputMZArray.size()-1]);

        sl->spectra.push_back(s);
    }


    vector<double> goldMZArray = parseDoubleArray(goldStandard[0].inputMZArray);
    vector<double> goldIntensityArray = parseDoubleArray(goldStandard[0].inputIntensityArray);

    // run spectral summation
    try
    {

        SpectrumListPtr calculator(new SpectrumList_ScanSummer(originalList,0.05,10));

        for (size_t i=0; i < calculator->size(); ++i) 
        {
            SpectrumPtr s = calculator->spectrum(i,true);
            vector<double>& mzs = s->getMZArray()->data;
            vector<double>& intensities = s->getIntensityArray()->data;
            Precursor& precursor = s->precursors[0];
            SelectedIon& selectedIon = precursor.selectedIons[0];
            double precursorMZ = selectedIon.cvParam(MS_selected_ion_m_z).valueAs<double>();
            double rTime = s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds();

            vector<double> goldMZArray = parseDoubleArray(goldStandard[i].inputMZArray);
            vector<double> goldIntensityArray = parseDoubleArray(goldStandard[i].inputIntensityArray);

            unit_assert(mzs.size() == goldMZArray.size());
            unit_assert(intensities.size() == goldIntensityArray.size());
            unit_assert(precursorMZ == goldStandard[i].inputPrecursorMZ);
            unit_assert(rTime == goldStandard[i].rTime);

            for (size_t j=0; j < mzs.size(); ++j)
            {
                unit_assert_equal(mzs[j], goldMZArray[j], 1e-5);
                unit_assert_equal(intensities[j], goldIntensityArray[j], 1e-5);
            }

        }

        
    
            
    }
    catch (exception& e)
    {
        cerr << "Test failed:\n" << e.what() << endl;
        ++failedTests;
    }
    return failedTests;
}

//.........这里部分代码省略.........
    double startMz, stopMz;
    msExperiment->getAcquisitionMassRange(startMz, stopMz);
    scan.scanWindows.push_back(ScanWindow(startMz, stopMz, MS_m_z));

    // decide whether to use Points or Peaks to populate data arrays
    bool doCentroid = msLevelsToCentroid.contains(msLevel);

    bool continuousData = spectrum->getDataIsContinuous();
    if (continuousData && !doCentroid)
        result->set(MS_profile_spectrum);
    else
    {
        result->set(MS_centroid_spectrum);
        doCentroid = continuousData;
    }

    /*if (experimentType == MRM)
    {
        MRMTransitions^ transitions = msExperiment->MRMTransitions;
        double q1mz = transitions[ie.transition]->Q1Mass;//ie.transition->first;
        double q3mz = transitions[ie.transition]->Q3Mass;
        double intensity = points[ie.transition]->Y;
        result->defaultArrayLength = 1;//ie.transition->second.size();

        Precursor precursor;
        SelectedIon selectedIon;

        selectedIon.set(MS_selected_ion_m_z, q1mz, MS_m_z);

        precursor.activation.set(MS_CID); // assume CID

        precursor.selectedIons.push_back(selectedIon);
        result->precursors.push_back(precursor);

        if (getBinaryData)
        {
            mzArray.resize(result->defaultArrayLength, q3mz);
            intensityArray.resize(result->defaultArrayLength, intensity);
        }
    }
    else*/
    {
        if (spectrum->getHasPrecursorInfo())
        {
            double selectedMz, intensity;
            int charge;
            spectrum->getPrecursorInfo(selectedMz, intensity, charge);

            Precursor precursor;
            SelectedIon selectedIon;

            selectedIon.set(MS_selected_ion_m_z, selectedMz, MS_m_z);
            if (charge > 0)
                selectedIon.set(MS_charge_state, charge);

            precursor.activation.set(MS_beam_type_collision_induced_dissociation); // assume beam-type CID since all ABI instruments that write WIFFs are either QqTOF or QqLIT

            precursor.selectedIons.push_back(selectedIon);

            // fix isolation width
            if (spectrum->getHasIsolationInfo())
            {
                double centerMz, lowerLimit, upperLimit;
                spectrum->getIsolationInfo(centerMz, lowerLimit, upperLimit);

                precursor.set(MS_isolation_window_target_m_z, centerMz, MS_m_z);
                precursor.set(MS_isolation_window_upper_offset, upperLimit);
                precursor.set(MS_isolation_window_lower_offset, lowerLimit);
            }


            result->precursors.push_back(precursor);
        }

        //virtual bool getHasIsolationInfo() const;
        //virtual void getIsolationInfo(double& centerMz, double& lowerLimit, double& upperLimit) const;


        //result->set(MS_lowest_observed_m_z, spectrum->getMinX(), MS_m_z);
        //result->set(MS_highest_observed_m_z, spectrum->getMaxX(), MS_m_z);
        result->set(MS_base_peak_intensity, spectrum->getBasePeakY(), MS_number_of_detector_counts);
        result->set(MS_base_peak_m_z, spectrum->getBasePeakX(), MS_m_z);
        result->set(MS_total_ion_current, spectrum->getSumY(), MS_number_of_detector_counts);

        if (getBinaryData)
        {
            result->setMZIntensityArrays(std::vector<double>(), std::vector<double>(), MS_number_of_detector_counts);
            BinaryDataArrayPtr mzArray = result->getMZArray();
            BinaryDataArrayPtr intensityArray = result->getIntensityArray();

            spectrum->getData(doCentroid, mzArray->data, intensityArray->data);
            if (doCentroid)
                result->set(MS_profile_spectrum); // let SpectrumList_PeakPicker know this was a profile spectrum
        }

        result->defaultArrayLength = spectrum->getDataSize(doCentroid);
    }

    return result;
}

PWIZ_API_DECL SpectrumPtr SpectrumList_PeakPicker::spectrum(size_t index, bool getBinaryData) const
{
    SpectrumPtr s;
    
    switch (mode_)
    {
        case 1:
            s = dynamic_cast<detail::SpectrumList_Thermo*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 2:
            s = dynamic_cast<detail::SpectrumList_Bruker*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 3:
            s = dynamic_cast<detail::SpectrumList_ABI*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 4:
            s = dynamic_cast<detail::SpectrumList_Agilent*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 5:
            s = dynamic_cast<detail::SpectrumList_ABI_T2D*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 0:
        default:
            s = inner_->spectrum(index, true);
            break;
    }

    if (!msLevelsToPeakPick_.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return s;

    vector<CVParam>& cvParams = s->cvParams;
    vector<CVParam>::iterator itr = std::find(cvParams.begin(), cvParams.end(), MS_profile_spectrum);

    // return non-profile spectra as-is
    // (could have been acquired as centroid, or vendor may have done the centroiding)
    if (itr == cvParams.end())
        return s;

    // make sure the spectrum has binary data
    if (!s->getMZArray().get() || !s->getIntensityArray().get())
        s = inner_->spectrum(index, true);

    // replace profile term with centroid term
    vector<CVParam>& cvParams2 = s->cvParams;
    *(std::find(cvParams2.begin(), cvParams2.end(), MS_profile_spectrum)) = MS_centroid_spectrum;

    try
    {
        vector<double>& mzs = s->getMZArray()->data;
        vector<double>& intensities = s->getIntensityArray()->data;
        vector<double> xPeakValues, yPeakValues;
        algorithm_->detect(mzs, intensities, xPeakValues, yPeakValues);
        mzs.swap(xPeakValues);
        intensities.swap(yPeakValues);
        s->defaultArrayLength = mzs.size();
    }
    catch(std::exception& e)
    {
        throw std::runtime_error(std::string("[SpectrumList_PeakPicker::spectrum()] Error picking peaks: ") + e.what());
    }

    s->dataProcessingPtr = dp_;
    return s;
}

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(),
//.........这里部分代码省略.........

      void mzDataCache::addKeyValue(const string& s, SpectrumPtr& spec) {
			vector<double>& mz = spec->getMZArray()->data;
			vector<double>& intensity = spec->getIntensityArray()->data;
			db.set((s + "mz").c_str(), (s + "mz").size() * sizeof (char), (char*) &mz[0], mz.size() * sizeof (double));
			db.set((s + "int").c_str(), (s + "int").size() * sizeof (char), (char*) &intensity[0], intensity.size() * sizeof (double));
		}