C++ DataPoint::SetObservable方法代码示例

//Use accept/reject method to create data
int Foam::GenerateData( int DataAmount )
{
	//if( newDataSet->GetDataNumber() < DataAmount ) newDataSet->ReserveDataSpace( DataAmount );
	for( int dataIndex = 0; dataIndex < DataAmount; ++dataIndex )
	{
		//Generate the discrete observables, and select the correct Foam generator to use
		int combinationIndex = 0;
		int incrementValue = 1;
		DataPoint * temporaryDataPoint = new DataPoint(allNames);
		for( int discreteIndex = int(discreteNames.size() - 1); discreteIndex >= 0; --discreteIndex )
		{
			//Create the discrete observable
			Observable * temporaryObservable = generationBoundary->GetConstraint( *discreteNames_ref[unsigned(discreteIndex)] )->CreateObservable(rootRandom);
			double currentValue = temporaryObservable->GetValue();
			temporaryDataPoint->SetObservable( *discreteNames_ref2[unsigned(discreteIndex)], temporaryObservable );
			delete temporaryObservable;

			//Calculate the index
			for (unsigned int valueIndex = 0; valueIndex < discreteValues[unsigned(discreteIndex)].size(); ++valueIndex )
			{
				if ( fabs( discreteValues[unsigned(discreteIndex)][valueIndex] - currentValue ) < DOUBLE_TOLERANCE )
				{
					combinationIndex += ( incrementValue * int(valueIndex) );
					incrementValue *= int(discreteValues[unsigned(discreteIndex)].size());
					break;
				}
			}
		}

		//Use the index calculated to select a Foam generator and generate an event with it
		Double_t* generatedEvent = new Double_t[ continuousNames.size() ];
		foamGenerators[unsigned(combinationIndex)]->MakeEvent();
		foamGenerators[unsigned(combinationIndex)]->GetMCvect(generatedEvent);

		//Store the continuous observables
		for (unsigned int continuousIndex = 0; continuousIndex < continuousNames.size(); ++continuousIndex )
		{
			string unit = generationBoundary->GetConstraint( *continuousNames_ref[continuousIndex] )->GetUnit();
			double newValue = minima[continuousIndex] + ( ranges[continuousIndex] * generatedEvent[continuousIndex] );
			temporaryDataPoint->SetObservable( *continuousNames_ref2[continuousIndex], newValue, unit );
			//cout << continuousNames[continuousIndex] << "\t" << newValue << "\t" << 0.0 << "\t" << unit << endl;
		}

		delete[] generatedEvent;
		//	Store the event
		temporaryDataPoint->SetDiscreteIndex(dataIndex);
		newDataSet->AddDataPoint(temporaryDataPoint);
	}

	//cout << "Destroying Generator(s)" << endl;
	//this->RemoveGenerator();
	return DataAmount;
}

//Return a list of data points
//Each should take the data average value of each continuous observable
//Each should represent one combination of possible discrete values
vector<DataPoint*> PhaseSpaceBoundary::GetDiscreteCombinations() const
{
	if( storedCombinationID == uniqueID && !StoredCombinations.empty() )
	{
		return StoredCombinations;
	}

	while( !StoredCombinations.empty() )
	{
		if( StoredCombinations.back() != NULL ) delete StoredCombinations.back();
		StoredCombinations.pop_back();
	}

	//Calculate all possible combinations of discrete observables
	vector<string> thisAllNames = this->GetAllNames();
	vector<vector<double> > discreteValues;
	vector<string> discreteNames, continuousNames;
	vector<vector<double> > discreteCombinations = StatisticsFunctions::DiscreteCombinations( &thisAllNames, this, discreteNames, continuousNames, discreteValues );

	(void) continuousNames;

	//Create the data points to return
	vector<DataPoint*> newDataPoints;

	DataPoint* tempPoint = this->GetMidPoint();
	if( tempPoint == NULL ) return newDataPoints;

	for( unsigned int combinationIndex = 0; combinationIndex < discreteCombinations.size(); ++combinationIndex )
	{
		DataPoint* templateDataPoint = new DataPoint( *tempPoint );

		//Output the discrete values for this combination
		for( unsigned int discreteIndex = 0; discreteIndex < discreteNames.size(); ++discreteIndex )
		{
			//Set the data point
			Observable* oldValue = templateDataPoint->GetObservable( discreteNames[discreteIndex] );

			Observable* newValue = new Observable( oldValue->GetName(), discreteCombinations[combinationIndex][discreteIndex], oldValue->GetUnit() );

			templateDataPoint->SetObservable( discreteNames[discreteIndex], newValue );
			delete newValue;
		}

		newDataPoints.push_back( templateDataPoint );
	}

	delete tempPoint;

	StoredCombinations = newDataPoints;

	storedCombinationID = uniqueID; 

	return newDataPoints;
}

vector<DataPoint*> RapidFitIntegrator::initGSLDataPoints( unsigned int number, vector<double> maxima, vector<double> minima, DataPoint* templateDataPoint, vector<string> doIntegrate )
{
	vector<DataPoint*> doEval_points;
	pthread_mutex_lock( &GSL_DATAPOINT_MANAGEMENT_THREADLOCK );

#ifdef __RAPIDFIT_USE_GSL

	unsigned int nDim = (unsigned) minima.size();

	unsigned int npoint = number;
	vector<double> * integrationPoints = new vector<double>[ nDim ];

	//pthread_mutex_lock( &gsl_mutex );
	cout << "Allocating GSL Integration Tool. nDim " << nDim << endl;
	gsl_qrng * q = NULL;
	try
	{
		//q = gsl_qrng_alloc( gsl_qrng_sobol, (unsigned)nDim );
		q = gsl_qrng_alloc( gsl_qrng_niederreiter_2, (unsigned)nDim );
	}
	catch(...)
	{
		cout << "Can't Allocate Integration Tool for GSL Integral." << endl;
		cout << " Dim: " << nDim << endl;
		exit(-742);
	}

	if( q == NULL )
	{
		cout << "Can't Allocate Integration Tool for GSL Integral." << endl;
		cout << " Dim: " << nDim << endl;
		exit(-741);
	}

	double* v = new double[ nDim ];
	for (unsigned int i = 0; i < npoint; i++)
	{
		gsl_qrng_get( q, v );
		for( unsigned int j = 0; j < nDim; j++)
		{
			integrationPoints[j].push_back( v[j] );
		}
	}
	delete v;
	gsl_qrng_free(q);
	//cout << "Freed GSL Integration Tool" << endl;
	//pthread_mutex_unlock( &gsl_mutex );

	/*vector<double> minima_v, maxima_v;
	  for( unsigned int i=0; i< nDim; ++i )
	  {
	  minima_v.push_back( minima[i] );
	  maxima_v.push_back( maxima[i] );
	  }*/
	//cout << "Constructing Functions" << endl;
	//IntegratorFunction* quickFunction = new IntegratorFunction( functionToWrap, NewDataPoint, doIntegrate, dontIntegrate, NewBoundary, componentIndex, minima_v, maxima_v );

	for (unsigned int i = 0; i < integrationPoints[0].size(); ++i)
	{
		double* point = new double[ nDim ];
		for ( unsigned int j = 0; j < nDim; j++)
		{
			//cout << doIntegrate[j] << " " << maxima[j] << " " << minima[j] << " " << integrationPoints[j][i] << endl;
			point[j] = integrationPoints[j][i]*(maxima[j]-minima[j])+minima[j];
		}

		DataPoint* thisPoint = new DataPoint( *templateDataPoint );

		for( unsigned int k=0; k< doIntegrate.size(); ++k )
		{
			thisPoint->SetObservable( doIntegrate[k], point[k], "noUnitsHere" );
		}

		doEval_points.push_back( thisPoint );
		//result += quickFunction->DoEval( point );
		//delete[] point;
	}

	delete[] integrationPoints;
#endif

	(void) maxima; (void) minima; (void) number;
	pthread_mutex_unlock( &GSL_DATAPOINT_MANAGEMENT_THREADLOCK );
	return doEval_points;
}

//Constructor with correct argument
MakeFoam::MakeFoam( IPDF * InputPDF, PhaseSpaceBoundary * InputBoundary, DataPoint * InputPoint ) : finishedCells(), centerPoints(), centerValues(), cellIntegrals(), integratePDF(InputPDF)
{
	//Make the container to hold the possible cells
	queue<PhaseSpaceBoundary*> possibleCells;
	PhaseSpaceBoundary* firstCell = InputBoundary;
	possibleCells.push(firstCell);

	//Make a list of observables to integrate over
	vector<string> doIntegrate, dontIntegrate;
	vector<string> pdfDontIntegrate = InputPDF->GetDoNotIntegrateList();
	StatisticsFunctions::DoDontIntegrateLists( InputPDF, InputBoundary, &(pdfDontIntegrate), doIntegrate, dontIntegrate );

	//Continue until all possible cells have been examined
	while ( !possibleCells.empty() )
	{
		//Remove the next possible cell
		PhaseSpaceBoundary* currentCell = possibleCells.front();
		possibleCells.pop();

		//Set up the histogram storage
		vector< vector<double> > histogramBinHeights, histogramBinMiddles, histogramBinMaxes;
		double normalisation = 0.0;
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
			vector<double> binHeights, binMiddles, binMaxes;
			IConstraint * temporaryConstraint = currentCell->GetConstraint( doIntegrate[observableIndex] );
			double minimum = temporaryConstraint->GetMinimum();
			double delta = ( temporaryConstraint->GetMaximum() - minimum ) / (double)HISTOGRAM_BINS;

			//Loop over bins
			for (int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
			{
				binHeights.push_back(0.0);
				binMiddles.push_back( minimum + ( delta * ( binIndex + 0.5 ) ) );
				binMaxes.push_back( minimum + ( delta * ( binIndex + 1.0 ) ) );
			}

			histogramBinHeights.push_back(binHeights);
			histogramBinMiddles.push_back(binMiddles);
			histogramBinMaxes.push_back(binMaxes);
		}

		//MC sample the cell, make projections, sort of
		for (int sampleIndex = 0; sampleIndex < MAXIMUM_SAMPLES; ++sampleIndex )
		{
			//Create a data point within the current cell
			DataPoint samplePoint( InputPoint->GetAllNames() );
			for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
			{
				//Generate random values to explore integrable observables
				IConstraint * temporaryConstraint = currentCell->GetConstraint( doIntegrate[observableIndex] );
				samplePoint.SetObservable( doIntegrate[observableIndex], temporaryConstraint->CreateObservable() );
			}
			for (unsigned int observableIndex = 0; observableIndex < dontIntegrate.size(); ++observableIndex )
			{
				//Use given values for unintegrable observables
				Observable * temporaryObservable = new Observable( *(InputPoint->GetObservable( dontIntegrate[observableIndex] )) );
				samplePoint.SetObservable( dontIntegrate[observableIndex], temporaryObservable );
				delete temporaryObservable;
			}

			//Evaluate the function at this point
			double sampleValue = InputPDF->Evaluate( &samplePoint );
			normalisation += sampleValue;

			//Populate the histogram
			for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
			{
				double observableValue = samplePoint.GetObservable( doIntegrate[observableIndex] )->GetValue();

				for ( int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
				{
					if ( observableValue < histogramBinMaxes[observableIndex][unsigned(binIndex)] )
					{
						histogramBinHeights[observableIndex][unsigned(binIndex)] += sampleValue;
						break;
					}
				}
			}
		}

		//Normalise the histograms
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
			for ( int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
			{
				histogramBinHeights[observableIndex][unsigned(binIndex)] /= normalisation;
			}
		}

		//Find the maximum gradient
		vector<double> midPoints;
		string maximumGradientObservable, unit;
		double maximumGradient=0.;
		double lowPoint=0.;
		double splitPoint=0.;
		double highPoint=0.;
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
//.........这里部分代码省略.........