C++ adcontrols::MassSpectrum类代码示例

bool
MappedSpectrum::transform( adcontrols::MassSpectrum& ms ) const
{
    auto& prop = ms.getMSProperty();

    int32_t nDelay = int32_t( ( delay_ / sampInterval_ ) + 0.5 );
    auto si = SamplingInfo( sampInterval_, delay_, nDelay, nSamples_, num_average_, 0 /* mode */ );
    si.fSampInterval( sampInterval_ );
    prop.setSamplingInfo( si );
    prop.setNumAverage( num_average_ );
    prop.setTrigNumber( trig_number_, trig_number_origin_ );
    prop.setTimeSinceEpoch( timeSinceEpoch_.second );

    ms.resize( data_.size() );
    ms.setCentroid( adcontrols::CentroidNative );
    //auto scanlaw = prop.scanLaw();
    
    for ( size_t idx = 0; idx < data_.size(); ++idx ) {
        double tof = this->time( idx );
        ms.setTime( idx, tof );
        //if ( scanlaw )
        //    ms.setMass( idx, ms.compute_mass( data_[ idx ].first ) ); 
        ms.setIntensity( idx, data_[ idx ].second );
    }
    return true;
}

bool
DataprocHandler::doCentroid( adcontrols::MSPeakInfo& pkInfo
                             , adcontrols::MassSpectrum& res
                             , const adcontrols::MassSpectrum& profile
                             , const adcontrols::CentroidMethod& m )
{
    adcontrols::CentroidProcess peak_detector;
    bool result = false;
    
    res.clone( profile, false );

    if ( peak_detector( m, profile ) ) {
        result = peak_detector.getCentroidSpectrum( res );
        pkInfo = peak_detector.getPeakInfo();
    }

    if ( profile.numSegments() > 0 ) {
        for ( size_t fcn = 0; fcn < profile.numSegments(); ++fcn ) {
            adcontrols::MassSpectrum centroid;
            result |= peak_detector( profile.getSegment( fcn ) );
            pkInfo.addSegment( peak_detector.getPeakInfo() );
            peak_detector.getCentroidSpectrum( centroid );
            res.addSegment( centroid );
        }
    }
    return result;
}

bool
QuanSampleProcessor::doCentroid( adcontrols::MSPeakInfo& pkInfo
                                 , adcontrols::MassSpectrum& res
                                 , const adcontrols::MassSpectrum& profile
                                 , const adcontrols::CentroidMethod& m )
{
    adcontrols::CentroidProcess peak_detector;
    bool result = false;
    
    res.clone( profile, false );
    
    if ( peak_detector( m, profile ) ) {
        result = peak_detector.getCentroidSpectrum( res );
        pkInfo = peak_detector.getPeakInfo();
    }
    
    if ( profile.numSegments() > 0 ) {
        for ( size_t fcn = 0; fcn < profile.numSegments(); ++fcn ) {
            auto centroid = std::make_shared< adcontrols::MassSpectrum >();
            result |= peak_detector( profile.getSegment( fcn ) );
            pkInfo.addSegment( peak_detector.getPeakInfo() );
            peak_detector.getCentroidSpectrum( *centroid );
            res << std::move( centroid ); 
        }
    }
    return result;
}

        void operator () ( const adcontrols::MassSpectrum& data
                           , adportable::counting::counting_result& result
                           , std::vector< double >& processed ) {

            assert ( method.algo_ == adcontrols::threshold_method::Differential );
            assert ( findPositive == ( method.slope == adcontrols::threshold_method::CrossUp ) );

            double level = method.threshold_level;

            adportable::counting::peak_finder< findPositive > finder;

            finder( data.getIntensityArray(), data.getIntensityArray() + data.size(), result.indices2(), level );
        }

bool
TimeDigitalHistogram::translate( adcontrols::MassSpectrum& sp
                                 , const TimeDigitalHistogram& hgrm )
{
    sp.setCentroid( adcontrols::CentroidNative );

    using namespace adcontrols::metric;

    // ext_trig_delay should be managed before came here.  (ex. histogram::move())

    double ext_trig_delay = hgrm.this_protocol_.delay_pulses().at( adcontrols::TofProtocol::EXT_ADC_TRIG ).first;
    
    adcontrols::MSProperty prop;
    adcontrols::SamplingInfo info( hgrm.xIncrement()
                                   , hgrm.initialXOffset() + ext_trig_delay
                                   , int32_t( ( hgrm.initialXOffset() + ext_trig_delay ) / hgrm.xIncrement() )  // delay
                                   , uint32_t( hgrm.actualPoints() ) // this is for acq. time range calculation
                                   , uint32_t( hgrm.trigger_count() )
                                   , hgrm.this_protocol_.mode() /* mode */);
    
    prop.setAcceleratorVoltage( 0 ); // empty

    prop.setSamplingInfo( info );
        
    prop.setTimeSinceInjection( hgrm.initialXTimeSeconds() );
    prop.setTimeSinceEpoch( hgrm.timeSinceEpoch().first );
    prop.setNumAverage( uint32_t( hgrm.trigger_count() ) );
    prop.setTrigNumber( uint32_t( hgrm.serialnumber().first ) );

    prop.setDataInterpreterClsid( "adcontrols::TimeDigitalHistogram" );

    TimeDigitalHistogram::device_data data( hgrm.this_protocol_ );
    std::string ar;
    adportable::binary::serialize<>()( data, ar );
    prop.setDeviceData( ar.data(), ar.size() );
        
    sp.setMSProperty( prop );

    size_t size = hgrm.size();

    sp.resize( size );
    size_t idx = 0;
    for ( auto it = hgrm.begin(); it != hgrm.end(); ++it, ++idx ) {
        sp.setTime( idx, it->first );
        sp.setIntensity( idx, it->second );
    }
    
    return true;
}

bool
rawdata::getSpectrum( int fcn, size_t pos, adcontrols::MassSpectrum& ms, uint32_t objid ) const
{
    auto it = std::find_if( conf_.begin(), conf_.end(), [=]( const adutils::AcquiredConf::data& c ){
            if ( objid == 0 )
                return c.trace_method == signalobserver::eTRACE_SPECTRA && c.trace_id == L"MS.PROFILE";
            else
                return c.objid == objid;
        });
    if ( it == conf_.end() )
        return false;

    adcontrols::translate_state state;
    uint64_t npos = npos0_ + pos;

    if ( fcn < 0 && fcnIdx_.size() > 1 ) {

        // find 'index' from <pos, fcn> array that indicates first 'pos' after protocol has been switched
        auto index = std::lower_bound( fcnIdx_.begin(), fcnIdx_.end(), npos
                                       , [] ( const std::pair< size_t, int >& a, size_t npos ) { return a.first < npos; } );
        if ( index == fcnIdx_.end() )
            return false;
        while ( index != fcnIdx_.begin() && index->first > npos )
            --index;
        int rep = int( npos - index->first );  // id within a replicates (rep'licates is the offset from (fcn=0,rep=0) spectrum)
        while ( index != fcnIdx_.begin() && index->second != 0 ) // find fcn=0
            --index;

        // read all protocols
        while ( ( state = fetchSpectrum( it->objid, it->dataInterpreterClsid, index->first + rep, ms, it->trace_id ) )
                == adcontrols::translate_indeterminate )
                if ( ++index == fcnIdx_.end() )  // move forward to next protocol (rep'licates is the offset for actual spectrum)
                    break;
    } else {

        state = fetchSpectrum( it->objid, it->dataInterpreterClsid, npos, ms, it->trace_id );

    }

    if ( ms.getMSProperty().dataInterpreterClsid() == 0 ) {
        // workaround for batchproc::import
        adcontrols::MSProperty prop = ms.getMSProperty();
        prop.setDataInterpreterClsid( adportable::utf::to_utf8( it->dataInterpreterClsid ).c_str() );
    }
    if ( fcn < 0 )
        return state == adcontrols::translate_complete;

    return state == adcontrols::translate_complete || state == adcontrols::translate_indeterminate;
}

bool
DataprocHandler::doMSCalibration( adcontrols::MSCalibrateResult& res
                                  , adcontrols::MassSpectrum& centroid
                                  , const adcontrols::MSCalibrateMethod& m
                                  , const adcontrols::MSAssignedMasses& assigned )
{
    using adcontrols::MSProperty;

    const double tolerance = m.massToleranceDa();
    const double threshold = centroid.getMaxIntensity() * m.minimumRAPercent() / 100;
    res.tolerance( tolerance );  // set tolerance in result
    res.threshold( threshold );  // set threshold in result

    std::map< size_t, size_t > mode_map;
    for ( adcontrols::MSAssignedMasses::vector_type::const_iterator it = assigned.begin(); it != assigned.end(); ++it ) 
        mode_map[ it->mode() ]++;
    // std::map<size_t, size_t>::iterator itMax = std::max_element( mode_map.begin(), mode_map.end() );
    // int mode = static_cast<int>(itMax->first);

    mass_calibrator calibrator( assigned, centroid.getMSProperty() );
    adcontrols::MSCalibration calib;
    if ( ! calibrator.polfit( calib, m.polynomialDegree() + 1 ) )
        return false;

    res.references( m.references() );
    res.calibration( calib );
    centroid.setCalibration( calib, true ); // m/z assign based on manually determined peaks

    // continue auto-assign
    assign_masses assign( tolerance, threshold );
    adcontrols::MSAssignedMasses assignedMasses;
	adcontrols::segment_wrapper< adcontrols::MassSpectrum > segments( centroid );
	for ( size_t n = 0; n < segments.size(); ++n ) {
		assign( assignedMasses, segments[n], m.references(), 0, static_cast<int>(n) );
	}

    mass_calibrator calibrator2( assignedMasses, centroid.getMSProperty() );
    if ( calibrator2.polfit( calib, m.polynomialDegree() + 1 ) ) {
        for ( auto it: assignedMasses )
            it.mass( calibrator2.compute_mass( it.time(), it.mode(), calib ) );

        centroid.setCalibration( calib, true );
        res.calibration( calib );
        res.assignedMasses( assignedMasses );

        return true;
    }
    return false;
}

bool
DataprocHandler::doMSCalibration( adcontrols::MSCalibrateResult& res
                                 , adcontrols::MassSpectrum& centroid
                                 , const adcontrols::MSCalibrateMethod& m )
{
    using adcontrols::MSProperty;

    res.calibration( centroid.calibration() );
    res.references( m.references() );
    double tolerance = m.massToleranceDa();
    double threshold = adcontrols::segments_helper::max_intensity( centroid ) * m.minimumRAPercent() / 100;
    res.tolerance( tolerance );
    res.threshold( threshold );

    assign_masses assigner( tolerance, threshold );

    adcontrols::MSAssignedMasses assignedMasses;
    
	adcontrols::segment_wrapper< adcontrols::MassSpectrum > segments( centroid );
    int n = 0;
	for ( auto seg: segments )
		assigner( assignedMasses, seg, res.references(), seg.mode(), n++ );

	res.assignedMasses( assignedMasses ); // set peak assign result

    // annotate each peak on spectrum
    doAnnotateAssignedPeaks( centroid, assignedMasses );

    mass_calibrator calibrator( assignedMasses, centroid.getMSProperty() );

    adcontrols::MSCalibration calib;
    if ( calibrator.polfit( calib, m.polynomialDegree() + 1 ) ) {
        for ( auto it: assignedMasses ) {
            double mass = calibrator.compute_mass( it.time(), it.mode(), calib );
            it.mass( mass );
        }
        res.calibration( calib );
        // res.assignedMasses( assignedMasses );
#if defined _DEBUG && 0
        calibresult_validation( res, centroid, threshold );
#endif
        return true;
    }
    return false;
}

//virtual
bool
datafile::getSpectrum( int /* fcn*/, size_t idx, adcontrols::MassSpectrum& ms, uint32_t /* objid */) const
{
    if ( unsigned( idx ) < spcfile_->number_of_subfiles() ) {

        const galactic::spchdr& hdr = *spcfile_->spchdr();
        const galactic::subhdr& sub = *spcfile_->subhdr( idx );

        std::pair< double, double > range = std::make_pair( hdr.ffirst(), hdr.flast() );
        const size_t npts = hdr.fnpts();
        ms.resize( npts );
        ms.setAcquisitionMassRange( range.first, range.second );
        
        for ( size_t i = 0; i < npts; ++i ) {
            ms.setMass( int(i), i * double(( range.second - range.first )) / ( npts - 1 ) + range.first );
            ms.setIntensity( int(i), sub[i] );
        }
        return true;
    }
	return false;
}

static void
calibresult_validation( const adcontrols::MSCalibrateResult& res
                        , const adcontrols::MassSpectrum& centroid
                        , double threshold )
{
    const adcontrols::MSReferences& ref = res.references();
    const adcontrols::MSAssignedMasses& assigned = res.assignedMasses();
    
    std::ofstream of( "massassign.txt" );
    of << "#\tm/z(observed)\ttof(us)\tintensity\t\tformula,\tm/z(exact)\tm/z(calibrated)\terror(mDa)" << std::endl;
    
    adcontrols::MSReferences::vector_type::const_iterator refIt = ref.begin();
    for ( adcontrols::MSAssignedMasses::vector_type::const_iterator it = assigned.begin(); it != assigned.end(); ++it, ++refIt ) {
        const adcontrols::MSAssignedMass& a = *it;
        
        std::string formula = adportable::string::convert( a.formula() );
        of << std::setprecision(8)
           << std::setw(4) << a.idMassSpectrum() << "\t" // id
           << std::setw(15) << std::fixed << centroid.getMass( a.idMassSpectrum() ) << "\t"           // m/z(observed)
           << std::scientific << centroid.getTime( a.idMassSpectrum() ) << "\t"      // tof
           << std::fixed << std::setprecision( 0 ) << centroid.getIntensity( a.idMassSpectrum() ) << "\t" // intensity
           << formula << "\t"
           << std::setprecision(8) << std::fixed   << it->exactMass() << "\t"                             // mass(exact)
           << std::fixed   << a.mass() << "\t"                                    // m/z(calibrated)
           << std::setprecision(1) << ( a.mass() - it->exactMass() ) * 1000 << "\t"  // error(mDa)
           << ( it->enable() ? "used" : "not used" ) 
           << std::endl;
    }
    const std::vector<double>& coeffs = res.calibration().coeffs();
    
    of << "#--------------------------- Calibration coefficients: " << std::endl;
    for ( size_t i = 0; i < coeffs.size(); ++i )
        of << std::scientific << std::setprecision(14) << coeffs[i] << std::endl;
    
    of << "#--------------------------- centroid peak list (#,mass,intensity)--------------------------" << std::endl;
    
    adcontrols::MSReferences::vector_type::const_iterator it = res.references().begin();
    for ( size_t i = 0; i < centroid.size(); ++i ) {
        if ( centroid.getIntensity( i ) > threshold ) {
            
            double mq = adcontrols::MSCalibration::compute( res.calibration().coeffs(), centroid.getTime( i ) );
            double mass = mq * mq;
            
            double error = 0;
            if ( it != res.references().end() && std::abs( it->exactMass() - mass ) < 0.2 ) {
                error = ( it->exactMass() - mass ) * 1000; // mDa
                ++it;
            }
            of << i << "\t"
               << std::setprecision(8) << std::fixed << centroid.getMass( i ) << "\t"
               << std::setprecision(8) << mass << "\t"
                    << std::setprecision(1) << centroid.getIntensityArray()[i] << std::endl;
        }
    }
}

//virtual
bool
datafile::getSpectrum( int /* fcn*/, size_t idx, adcontrols::MassSpectrum& ms, uint32_t /* objid */) const
{
    if ( lrpfile_ && unsigned( idx ) < lrpfile_->number_of_spectra() ) {

        if ( auto msdata = (*lrpfile_)[ idx ] ) {

            std::vector< double > time, intens;
            if ( lrpfile_->getMS( *msdata, time, intens ) ) {

                ms.resize( time.size() );
                ms.setMassArray( time.data() );
                ms.setIntensityArray( intens.data() );

                ms.setAcquisitionMassRange( time.front(), time.back() );

                return true;
            }
        }
    }
    return false;
}

bool
MSChromatogramExtractor::impl::doMSLock( adcontrols::lockmass::mslock& mslock
                                       , const adcontrols::MassSpectrum& centroid
                                       , const adcontrols::MSLockMethod& m )
{
    // TODO: consider how to handle segmented spectrum -- current impl is always process first 
    adcontrols::MSFinder find( m.tolerance( m.toleranceMethod() ), m.algorithm(), m.toleranceMethod() );

    for ( auto& msref : msrefs_ ) {
        size_t idx = find( centroid, msref.second );
        if ( idx != adcontrols::MSFinder::npos ) 
            mslock << adcontrols::lockmass::reference( msref.first, msref.second, centroid.getMass( idx ), centroid.getTime( idx ) );
    }

    if ( mslock.fit() ) {
        // mslock( centroid, true );
        return true;
    }
    return false;
}

bool
fsio::save_mscalibfile( adfs::filesystem& fs, const adcontrols::MassSpectrum& ms )
{
    return save( fs, ms, ms.dataClass(), L"/MSCalibration" );
}

void
MSPropertyForm::render( std::ostream& o, const adcontrols::MassSpectrum& ms )
{
    using adportable::utf;
    using namespace adcontrols::metric;

    static const char * const polarities [] = { "Indeterminant", "Positive", "Negative", "Mixed" };

    adcontrols::segment_wrapper< const adcontrols::MassSpectrum > segments( ms );

    std::vector< std::pair< std::string, std::string > > temp;
    temp.push_back( std::make_pair( "Spectrum", ( ms.isCentroid() ? "Centroid" : "Profile" ) ) );
    temp.push_back( std::make_pair( "Polarity", polarities[ ms.polarity() ] ) );

    o << "<hr>";
    o << "<table border=\"1\" cellpadding=\"1\">";
    o << "<caption>Descriptions</caption>";
    o << "<tr>";
    o << "<th>key</th>"
      << "<th>description</th>"
      << "</tr>";
    for ( auto it: temp ) {
        o << "<tr>"
          << "<td>" << it.first << "</td>"
          << "<td>" << it.second << "</td>"
          << "</tr>";
    }
    for ( size_t i = 0; i < ms.getDescriptions().size(); ++i ) {
        auto desc = ms.getDescriptions()[ i ];
        o << "<tr>"
          << "<td>" << utf::to_utf8( desc.key() ) << "</td>"
          << "<td>" << utf::to_utf8( desc.text() ) << "</td>"
          << "</tr>";
    }
    o << "</table>";

    std::pair< double, double > massrange = ms.getAcquisitionMassRange();
    o << "Acquisition mass range: " << boost::format( "%.3lf -- %.3lf" ) % massrange.first % massrange.second;

    o << "<table border=\"1\" cellpadding=\"4\">";
    o << "<caption>Acquisition and calibration parameter(s)</caption>";
    o << "<tr>";
    o << "<th>#seg</th>"
      << "<th>samp. interval(ps)</th>"
      << "<th>sampling start (&mu;s)</th>"
      << "<th>sampling end (&mu;s)</th>"
      << "<th>num. samples</th>"
      << "<th>num. average</th>"
      << "<th>mode</th>"
      << "<th>classid</th>"
      << "</tr>";
    int n = 0;
    for ( auto& m: segments ) {
        const adcontrols::MSCalibration& calib =  m.calibration();
        size_t nrowspan = calib.coeffs().empty() ? 1 : 2;
        const adcontrols::MSProperty& prop = m.getMSProperty();
        const adcontrols::MSProperty::SamplingInfo& info = prop.getSamplingInfo();
        double start_delay = scale_to<double, micro>( info.sampInterval * info.nSamplingDelay, pico );
        double time_end = scale_to<double, micro>( info.sampInterval * ( info.nSamplingDelay + info.nSamples ), pico );
        o << "<tr>"
          << boost::format( "<td rowspan=\"%1%\">" ) % nrowspan << n++ << "</td>"
          << "<td>" << info.sampInterval << "</td>"
          << "<td>" << boost::format( "%.4lf&mu;s" ) % start_delay << "</td>"
          << "<td>" << boost::format( "%.4lf&mu;s" ) % time_end << "</td>"
          << "<td>" << info.nSamples << "</td>"
          << "<td>" << info.nAverage << "</td>"
          << "<td>" << info.mode << "</td>"
		  << "<td>" << prop.dataInterpreterClsid() << "</td>"
          << "</tr>";

        if ( ! calib.coeffs().empty() ) {
            //-----------------------------
            o << "<tr>";
            o << "<td colspan=7><b>Calibration ID:</b><i>" << utf::to_utf8( calib.calibId() ) << "</i>"
              << "     " << calib.date();
            o << "<hr>";
            o << "&radic;<span style=\"text-decoration: overline\">&nbsp;<i>m/z</i></span> = "
              << boost::format( "%.14lf" ) % calib.coeffs()[0] << " + ";
            for ( size_t i = 1; i < calib.coeffs().size(); ++i )
                o << boost::format( "\t%.14lf &times; t<sup>%d</sup>" ) % calib.coeffs()[i] % i;

            if ( ! calib.t0_coeffs().empty() ) {
                o << "<br>"
                  << "T<sub>0</sub> = " << calib.t0_coeffs()[0];
                for ( size_t i = 1; i < calib.t0_coeffs().size(); ++i )
                    o << boost::format( "\t%.14lf &times; &radic;<span style=\"text-decoration: overline\">&nbsp;<i>m/z</i></span><sup>%d</sup>" )
                        % calib.coeffs()[i] % i;
                if ( calib.algorithm() == adcontrols::MSCalibration::MULTITURN_NORMALIZED )
                    o << "\t(MULTITURN_NORMAILZED algorithm)";
            }
            o << "</tr>";
            //-----------------------------            
		}
    }
    o << "</table>";
    o << "<hr>";

    try {
        // device (averager) dependent data (require data interpreter)
        std::vector < std::pair< std::string, std::string > > textv;
//.........这里部分代码省略.........

bool
fsio::load_mscalibfile( adfs::filesystem& fs, adcontrols::MassSpectrum& ms )
{
    return load( fs, ms, ms.dataClass(), L"/MSCalibration" );
}