C++ MantidVecPtr类代码示例

/** Reads from the third line of the input file to the end assuming it contains
*  2D data
*  @param firstLine :: the second line in the file
*  @return a workspace containing the loaded data
*  @throw NotFoundError if there is compulsulary data is missing from the file
*  @throw invalid_argument if there is an inconsistency in the header
* information
*/
const MatrixWorkspace_sptr LoadRKH::read2D(const std::string &firstLine) {
  g_log.information()
      << "file appears to contain 2D information, reading in 2D data mode\n";

  MatrixWorkspace_sptr outWrksp;
  MantidVec axis0Data;
  Progress prog(read2DHeader(firstLine, outWrksp, axis0Data));
  const size_t nAxis1Values = outWrksp->getNumberHistograms();

  for (size_t i = 0; i < nAxis1Values; ++i) {
    // set the X-values to the common bin values we read above
    MantidVecPtr toPass;
    toPass.access() = axis0Data;
    outWrksp->setX(i, toPass);

    // now read in the Y values
    MantidVec &YOut = outWrksp->dataY(i);
    for (double &value : YOut) {
      m_fileIn >> value;
    }
    prog.report("Loading Y data");
  } // loop on to the next spectrum

  // the error values form one big block after the Y-values
  for (size_t i = 0; i < nAxis1Values; ++i) {
    MantidVec &EOut = outWrksp->dataE(i);
    for (double &value : EOut) {
      m_fileIn >> value;
    }
    prog.report("Loading error estimates");
  } // loop on to the next spectrum

  return outWrksp;
}

/** Creates the output workspace, setting the axes according to the input
 * binning parameters
 *  @param[in]  inputWorkspace The input workspace
 *  @param[in]  binParams The bin parameters from the user
 *  @param[out] newAxis        The 'vertical' axis defined by the given
 * parameters
 *  @return A pointer to the newly-created workspace
 */
API::MatrixWorkspace_sptr
SofQW::setUpOutputWorkspace(API::MatrixWorkspace_const_sptr inputWorkspace,
                            const std::vector<double> &binParams,
                            std::vector<double> &newAxis) {
  // Create vector to hold the new X axis values
  MantidVecPtr xAxis;
  xAxis.access() = inputWorkspace->readX(0);
  const int xLength = static_cast<int>(xAxis->size());
  // Create a vector to temporarily hold the vertical ('y') axis and populate
  // that
  const int yLength = static_cast<int>(
      VectorHelper::createAxisFromRebinParams(binParams, newAxis));

  // Create the output workspace
  MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
      inputWorkspace, yLength - 1, xLength, xLength - 1);
  // Create a numeric axis to replace the default vertical one
  Axis *const verticalAxis = new BinEdgeAxis(newAxis);
  outputWorkspace->replaceAxis(1, verticalAxis);

  // Now set the axis values
  for (int i = 0; i < yLength - 1; ++i) {
    outputWorkspace->setX(i, xAxis);
  }

  // Set the axis units
  verticalAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  verticalAxis->title() = "|Q|";

  // Set the X axis title (for conversion to MD)
  outputWorkspace->getAxis(0)->title() = "Energy transfer";

  return outputWorkspace;
}

/**
 * @returns The bin bounadries for the new workspace
 */
MantidVecPtr CreateSimulationWorkspace::createBinBoundaries() const {
  const std::vector<double> rbparams = getProperty("BinParams");
  MantidVecPtr binsPtr;
  MantidVec &newBins = binsPtr.access();
  const int numBoundaries =
      Mantid::Kernel::VectorHelper::createAxisFromRebinParams(rbparams,
                                                              newBins);
  if (numBoundaries <= 2) {
    throw std::invalid_argument(
        "Error in BinParams - Gave invalid number of bin boundaries: " +
        boost::lexical_cast<std::string>(numBoundaries));
  }
  return binsPtr;
}

/** Creates the output workspace, its size, units, etc.
*  @param binParams the bin boundary specification using the same same syntax as param the Rebin algorithm
*  @return A pointer to the newly-created workspace
*/
API::MatrixWorkspace_sptr Q1D2::setUpOutputWorkspace(const std::vector<double> & binParams) const
{
  // Calculate the output binning
  MantidVecPtr XOut;
  size_t sizeOut = static_cast<size_t>(
    VectorHelper::createAxisFromRebinParams(binParams, XOut.access()));

  // Now create the output workspace
  MatrixWorkspace_sptr outputWS = WorkspaceFactory::Instance().create(m_dataWS,1,sizeOut,sizeOut-1);
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("MomentumTransfer");
  outputWS->setYUnitLabel("1/cm");

  // Set the X vector for the output workspace
  outputWS->setX(0, XOut);
  outputWS->isDistribution(true);

  outputWS->getSpectrum(0)->clearDetectorIDs();
  outputWS->getSpectrum(0)->setSpectrumNo(1);

  return outputWS;
}

/** Resizes the workspace to contain the number of spectra/events lists given.
 *  Any existing eventlists will be cleared first.
 *  Spectrum numbers will be set to count from 1
 *  @param numSpectra The number of spectra to resize the workspace to
 */
void EventWorkspace::resizeTo(const std::size_t numSpectra) {
  // Remove all old EventLists and resize the vector
  this->clearData();
  data.resize(numSpectra);
  m_noVectors = numSpectra;
  for (size_t i = 0; i < numSpectra; ++i) {
    data[i] = new EventList(mru, static_cast<specid_t>(i + 1));
  }

  // Put on a default set of X vectors, with one bin of 0 & extremely close to
  // zero
  MantidVecPtr xVals;
  MantidVec &x = xVals.access();
  x.resize(2, 0.0);
  // Move the rhs very,very slightly just incase something doesn't like them
  // being the same
  x[1] = std::numeric_limits<double>::min();
  this->setAllX(xVals);

  // Clearing the MRU list is a good idea too.
  this->clearMRU();
}

/** Convert the workspace units according to a simple output = a * (input^b) relationship
 *  @param outputWS :: the output workspace
 *  @param factor :: the conversion factor a to apply
 *  @param power :: the Power b to apply to the conversion
 */
void ConvertUnits::convertQuickly(API::MatrixWorkspace_sptr outputWS, const double& factor, const double& power)
{
  Progress prog(this,0.2,1.0,m_numberOfSpectra);
  int64_t numberOfSpectra_i = static_cast<int64_t>(m_numberOfSpectra); // cast to make openmp happy

  // See if the workspace has common bins - if so the X vector can be common
  // First a quick check using the validator
  CommonBinsValidator sameBins;
  bool commonBoundaries = false;
  if ( sameBins.isValid(outputWS) == "" )
  {
    commonBoundaries =  WorkspaceHelpers::commonBoundaries(outputWS);
    // Only do the full check if the quick one passes
    if (commonBoundaries)
    {
      // Calculate the new (common) X values
      MantidVec::iterator iter;
      for (iter = outputWS->dataX(0).begin(); iter != outputWS->dataX(0).end(); ++iter)
      {
        *iter = factor * std::pow(*iter,power);
      }

      MantidVecPtr xVals;
      xVals.access() = outputWS->dataX(0);

      PARALLEL_FOR1(outputWS)
      for (int64_t j = 1; j < numberOfSpectra_i; ++j)
      {
        PARALLEL_START_INTERUPT_REGION
        outputWS->setX(j,xVals);
        prog.report("Convert to " + m_outputUnit->unitID());
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
      if (!m_inputEvents) // if in event mode the work is done
        return;
    }
  }

    /**
     * Execute the algorithm
     */
    void ExtractMask::exec()
    {
      MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
      Geometry::Instrument_const_sptr instr = inputWS->getInstrument();

      // convert input to a mask workspace
      DataObjects::MaskWorkspace_const_sptr inputMaskWS
          = boost::dynamic_pointer_cast<const DataObjects::MaskWorkspace>(inputWS);
      bool inputWSIsSpecial = bool(inputMaskWS);
      if (inputWSIsSpecial)
      {
        g_log.notice() << "Input workspace is a MaskWorkspace.\n";
      }

      DataObjects::MaskWorkspace_sptr maskWS;
      // List masked of detector IDs
      std::vector<detid_t> detectorList;

      if (instr)
      {
        const int nHist = static_cast<int>(inputWS->getNumberHistograms());

        // Create a new workspace for the results, copy from the input to ensure that we copy over the instrument and current masking
        maskWS = DataObjects::MaskWorkspace_sptr(new DataObjects::MaskWorkspace(inputWS));
        maskWS->setTitle(inputWS->getTitle());

        Progress prog(this,0.0,1.0,nHist);

        MantidVecPtr xValues;
        xValues.access() = MantidVec(1, 0.0);

        PARALLEL_FOR2(inputWS, maskWS)
        for( int i = 0; i < nHist; ++i )
        {
          PARALLEL_START_INTERUPT_REGION

          bool inputIsMasked(false);
          IDetector_const_sptr inputDet;
          try
          {
            inputDet = inputWS->getDetector(i);
            if (inputWSIsSpecial) {
              inputIsMasked = inputMaskWS->isMaskedIndex(i);
            }
            // special workspaces can mysteriously have the mask bit set
            // but only check if we haven't already decided to mask the spectrum
            if( !inputIsMasked && inputDet->isMasked() )
            {
              inputIsMasked = true;
            }

            if (inputIsMasked)
            {
              detid_t id = inputDet->getID();
              PARALLEL_CRITICAL(name)
              {
                detectorList.push_back(id);
              }
            }
          }
          catch(Kernel::Exception::NotFoundError &)
          {
            inputIsMasked = false;
          }

          maskWS->setMaskedIndex(i, inputIsMasked);

          prog.report();

          PARALLEL_END_INTERUPT_REGION
        }
        PARALLEL_CHECK_INTERUPT_REGION

        // Clear all the "masked" bits on the output masked workspace
        Geometry::ParameterMap & pmap = maskWS->instrumentParameters();
        pmap.clearParametersByName("masked");
      }
      else // no instrument
      {
        // TODO should fill this in
        throw std::runtime_error("No instrument");

//.........这里部分代码省略.........
    }
  }

  size_t iPart = 0;
  if (m_spec_max != Mantid::EMPTY_INT()) {
    uint32_t ifirst = pixel_id[0];
    range_check out_range(m_spec_min, m_spec_max, id_to_wi);
    auto newEnd = std::remove_if(pixel_id.begin(), pixel_id.end(), out_range);
    pixel_id.erase(newEnd, pixel_id.end());
    // check if beginning or end of array was erased
    if (ifirst != pixel_id[0])
      iPart = m_numPixels - pixel_id.size();
    m_numPixels = pixel_id.size();
    if (m_numPixels == 0) {
      file->close();
      m_fileMutex.unlock();
      g_log.warning() << "No pixels from " << bankName << std::endl;
      return;
    };
  }
  // Load the TOF vector
  std::vector<float> tof;
  file->readData(m_axisField, tof);
  size_t m_numBins = tof.size() - 1;
  if (tof.size() <= 1) {
    file->close();
    m_fileMutex.unlock();
    g_log.warning() << "Invalid " << m_axisField << " data in " << bankName
                    << std::endl;
    return;
  }

  // Make a shared pointer
  MantidVecPtr Xptr;
  MantidVec &X = Xptr.access();
  X.resize(tof.size(), 0);
  X.assign(tof.begin(), tof.end());

  // Load the data. Coerce ints into double.
  std::string errorsField = "";
  std::vector<double> data;
  file->openData(m_dataField);
  file->getDataCoerce(data);
  if (file->hasAttr("errors"))
    file->getAttr("errors", errorsField);
  file->closeData();

  // Load the errors
  bool hasErrors = !errorsField.empty();
  std::vector<double> errors;
  if (hasErrors) {
    try {
      file->openData(errorsField);
      file->getDataCoerce(errors);
      file->closeData();
    } catch (...) {
      g_log.information() << "Error loading the errors field, '" << errorsField
                          << "' for bank " << bankName
                          << ". Will use sqrt(counts). " << std::endl;
      hasErrors = false;
    }
  }

  /*if (data.size() != m_numBins * m_numPixels)
  { file->close(); m_fileMutex.unlock(); g_log.warning() << "Invalid size of '"
  << m_dataField << "' data in " << bankName << std::endl; return; }

/** Execute the algorithm.
 */
void ResampleX::exec() {
  // generically having access to the input workspace is a good idea
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
  bool inPlace = (inputWS == outputWS); // Rebinning in-place
  m_isDistribution = inputWS->isDistribution();
  m_isHistogram = inputWS->isHistogramData();
  int numSpectra = static_cast<int>(inputWS->getNumberHistograms());

  // the easy parameters
  m_useLogBinning = getProperty("LogBinning");
  m_numBins = getProperty("NumberBins");
  m_preserveEvents = getProperty("PreserveEvents");

  // determine the xmin/xmax for the workspace
  vector<double> xmins = getProperty("XMin");
  vector<double> xmaxs = getProperty("XMax");
  string error = determineXMinMax(inputWS, xmins, xmaxs);
  if (!error.empty())
    throw std::runtime_error(error);

  bool common_limits = true;
  {
    double xmin_common = xmins[0];
    double xmax_common = xmaxs[0];
    for (size_t i = 1; i < xmins.size(); ++i) {
      if (xmins[i] != xmin_common) {
        common_limits = false;
        break;
      }
      if (xmaxs[i] != xmax_common) {
        common_limits = false;
        break;
      }
    }
  }
  if (common_limits) {
    g_log.debug() << "Common limits between all spectra\n";
  } else {
    g_log.debug() << "Does not have common limits between all spectra\n";
  }

  // start doing actual work
  EventWorkspace_const_sptr inputEventWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
  if (inputEventWS != NULL) {
    if (m_preserveEvents) {
      EventWorkspace_sptr outputEventWS =
          boost::dynamic_pointer_cast<EventWorkspace>(outputWS);
      if (inPlace) {
        g_log.debug() << "Rebinning event workspace in place\n";
      } else {
        g_log.debug() << "Rebinning event workspace out of place\n";

        // copy the event workspace to a new EventWorkspace
        outputEventWS = boost::dynamic_pointer_cast<EventWorkspace>(
            API::WorkspaceFactory::Instance().create(
                "EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
        // copy geometry over.
        API::WorkspaceFactory::Instance().initializeFromParent(
            inputEventWS, outputEventWS, false);
        // copy over the data as well.
        outputEventWS->copyDataFrom((*inputEventWS));
      }

      if (common_limits) {
        // get the delta from the first since they are all the same
        MantidVecPtr xValues;
        double delta =
            this->determineBinning(xValues.access(), xmins[0], xmaxs[0]);
        g_log.debug() << "delta = " << delta << "\n";
        outputEventWS->setAllX(xValues);
      } else {
        // initialize progress reporting.
        Progress prog(this, 0.0, 1.0, numSpectra);

        // do the rebinning
        PARALLEL_FOR2(inputEventWS, outputWS)
        for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
          PARALLEL_START_INTERUPT_REGION
          MantidVec xValues;
          double delta = this->determineBinning(xValues, xmins[wkspIndex],
                                                xmaxs[wkspIndex]);
          g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta
                        << " xmin=" << xmins[wkspIndex]
                        << " xmax=" << xmaxs[wkspIndex] << "\n";
          outputEventWS->getSpectrum(wkspIndex)->setX(xValues);
          prog.report(name()); // Report progress
          PARALLEL_END_INTERUPT_REGION
        }
        PARALLEL_CHECK_INTERUPT_REGION
      }

      this->setProperty(
          "OutputWorkspace",
          boost::dynamic_pointer_cast<MatrixWorkspace>(outputEventWS));
    }    // end if (m_preserveEvents)
    else // event workspace -> matrix workspace
//.........这里部分代码省略.........

/// Exec function
void CreateWorkspace::exec()
{
    // Contortions to get at the vector in the property without copying it
    const Property * const dataXprop = getProperty("DataX");
    const Property * const dataYprop = getProperty("DataY");
    const Property * const dataEprop = getProperty("DataE");
    const std::vector<double>& dataX = *dynamic_cast<const ArrayProperty<double>*>(dataXprop);
    const std::vector<double>& dataY = *dynamic_cast<const ArrayProperty<double>*>(dataYprop);
    const std::vector<double>& dataE = *dynamic_cast<const ArrayProperty<double>*>(dataEprop);

    const int nSpec = getProperty("NSpec");
    const std::string xUnit = getProperty("UnitX");
    const std::string vUnit = getProperty("VerticalAxisUnit");
    const std::vector<std::string> vAxis = getProperty("VerticalAxisValues");
    std::string parentWorkspace = getPropertyValue("ParentWorkspace");

    if ( ( vUnit != "SpectraNumber" ) && ( static_cast<int>(vAxis.size()) != nSpec ) )
    {
        throw std::invalid_argument("Number of y-axis labels must match number of histograms.");
    }

    // Verify length of vectors makes sense with NSpec
    if ( ( dataY.size() % nSpec ) != 0 )
    {
        throw std::invalid_argument("Length of DataY must be divisible by NSpec");
    }
    const std::size_t ySize = dataY.size() / nSpec;

    // Check whether the X values provided are to be re-used for (are common to) every spectrum
    const bool commonX( dataX.size() == ySize || dataX.size() == ySize+1 );

    std::size_t xSize;
    MantidVecPtr XValues;
    if ( commonX )
    {
        xSize = dataX.size();
        XValues.access() = dataX;
    }
    else
    {
        if ( dataX.size() % nSpec != 0 )
        {
            throw std::invalid_argument("Length of DataX must be divisible by NSpec");
        }

        xSize = static_cast<int>(dataX.size()) / nSpec;
        if ( xSize < ySize || xSize > ySize + 1 )
        {
            throw std::runtime_error("DataX width must be as DataY or +1");
        }

    }

    const bool dataE_provided = !dataE.empty();
    if ( dataE_provided && dataY.size() != dataE.size() )
    {
        throw std::runtime_error("DataE (if provided) must be the same size as DataY");
    }

    MatrixWorkspace_sptr parentWS;
    if (!parentWorkspace.empty())
    {
        try
        {
            parentWS = boost::dynamic_pointer_cast<MatrixWorkspace>( AnalysisDataService::Instance().retrieve(parentWorkspace) );
        }
        catch(...)
        {
            g_log.warning("Parent workspace not found");
            // ignore parent workspace
        }
    }

    // Create the OutputWorkspace
    MatrixWorkspace_sptr outputWS;
    if (parentWS)
    {
        // if parent is defined use it to initialise the workspace
        outputWS = WorkspaceFactory::Instance().create(parentWS, nSpec, xSize, ySize);
    }
    else
    {
        // otherwise create a blank workspace
        outputWS = WorkspaceFactory::Instance().create("Workspace2D", nSpec, xSize, ySize);
    }

    Progress progress(this,0,1,nSpec);

    PARALLEL_FOR1(outputWS)
    for ( int i = 0; i < nSpec; i++ )
    {
        PARALLEL_START_INTERUPT_REGION

        const std::vector<double>::difference_type xStart = i*xSize;
        const std::vector<double>::difference_type xEnd = xStart + xSize;
        const std::vector<double>::difference_type yStart = i*ySize;
        const std::vector<double>::difference_type yEnd = yStart + ySize;

        // Just set the pointer if common X bins. Otherwise, copy in the right chunk (as we do for Y).
//.........这里部分代码省略.........

	/**
	* Create histogram workspace
	* @returns Created workspace
	*/
	MatrixWorkspace_sptr LoadFITS::initAndPopulateHistogramWorkspace()
	{		
		MantidVecPtr x;
		x.access().resize(m_allHeaderInfo.size() + 1);

		// Init time bins
		double binCount = 0;
		for(size_t i=0;i<m_allHeaderInfo.size() + 1; ++i)
		{
			x.access()[i] = binCount;
			if(i != m_allHeaderInfo.size()) binCount += m_allHeaderInfo[i].timeBin;
		}

		size_t spectraCount = 0;
		if(m_allHeaderInfo[0].numberOfAxis > 0) spectraCount += m_allHeaderInfo[0].axisPixelLengths[0];

		// Presumably 2 axis, but futureproofing.
		for(int i=1;i<m_allHeaderInfo[0].numberOfAxis;++i)
		{
			spectraCount *= m_allHeaderInfo[0].axisPixelLengths[i];
		}

		MatrixWorkspace_sptr retVal(new DataObjects::Workspace2D);
		retVal->initialize(spectraCount, m_allHeaderInfo.size()+1, m_allHeaderInfo.size());
				
		IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
 
		try 
		{
			std::string directoryName = Kernel::ConfigService::Instance().getInstrumentDirectory();
			directoryName = directoryName + "/IMAT_Definition.xml";
			
			loadInst->setPropertyValue("Filename", directoryName);
			loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", retVal);
			loadInst->execute();
		} 
		catch (std::exception & ex) 
		{
			g_log.information("Cannot load the instrument definition. " + string(ex.what()) );
		}

		int bitsPerPixel = m_allHeaderInfo[0].bitsPerPixel; // assumes all files have the same, which they should. 
		vector<vector<double> > yVals(spectraCount, std::vector<double>(m_binChunkSize));
		vector<vector<double> > eVals(spectraCount, std::vector<double>(m_binChunkSize));
		
		// allocate memory to contain the data section of the file:
		void * bufferAny = NULL;        
		bufferAny = malloc ((bitsPerPixel/8)*spectraCount);

		if (bufferAny == NULL) 
		{
			throw std::runtime_error("FITS loader couldn't allocate enough memory to run. Try a smaller chunk size.");	
		}

		size_t steps = static_cast<size_t>(ceil(m_allHeaderInfo.size()/m_binChunkSize));
		Progress prog(this,0.0,1.0,steps);
		
		// Load a chunk of files at a time into workspace
		try
		{
			for(size_t i=0; i<m_allHeaderInfo.size(); i+=m_binChunkSize)
			{
				loadChunkOfBinsFromFile(retVal, yVals, eVals, bufferAny, x, spectraCount, bitsPerPixel, i);
				prog.report(name());
			}
		}
		catch(...)
		{
			// Exceptions should be handled internally, but catch here to free any memory. Belt and braces.
			free(bufferAny);
			g_log.error("FITS Loader unable to correctly parse files.");
			throw std::runtime_error("FITS loader unable to correctly parse files.");
		}

		// Memory no longer needed
		free (bufferAny);  
		
		retVal->mutableRun().addProperty("Filename", m_allHeaderInfo[0].filePath);

		// Set the Unit of the X Axis
		try
		{
			retVal->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
		}
		catch ( Exception::NotFoundError & )
		{
			retVal->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
			Unit_sptr unit = retVal->getAxis(0)->unit();
			boost::shared_ptr<Units::Label> label = boost::dynamic_pointer_cast<Units::Label>(unit);
			label->setLabel("TOF", "TOF");
		}

		retVal->setYUnit("Counts");
		retVal->setTitle("Test Workspace");

		return retVal;
//.........这里部分代码省略.........