<commit_before>/** @file CPropertyExtractor.cpp
 *
 *  @author <NAME>
 */

#include "../include/CPropertyExtractor.h"
#include "../include/Functions.h"
#include "../include/Util.h"

using namespace std;

namespace TBTK{

CPropertyExtractor::CPropertyExtractor(ChebyshevSolver *cSolver,
					int numCoefficients,
					int energyResolution,
					bool useGPUToCalculateCoefficients,
					bool useGPUToGenerateGreensFunctions,
					bool useLookupTable,
					double lowerBound,
					double upperBound){
	this->cSolver = cSolver;
	this->numCoefficients = numCoefficients;
	this->energyResolution = energyResolution;
	this->lowerBound = lowerBound;
	this->upperBound = upperBound;
	this->useGPUToCalculateCoefficients = useGPUToCalculateCoefficients;
	this->useGPUToGenerateGreensFunctions = useGPUToGenerateGreensFunctions;
	this->useLookupTable = useLookupTable;

	if(useLookupTable){
		cSolver->generateLookupTable(numCoefficients, energyResolution, lowerBound, upperBound);
		if(useGPUToGenerateGreensFunctions)
			cSolver->loadLookupTableGPU();
	}
	else if(useGPUToGenerateGreensFunctions){
		cout << "Error in PropertyExtractorChebyshev: useLookupTable cannot be false if useGPUToGenerateGreensFunction is true.\n";
		exit(1);
	}
}

CPropertyExtractor::~CPropertyExtractor(){
	if(useGPUToGenerateGreensFunctions)
		cSolver->destroyLookupTableGPU();
}

complex<double>* CPropertyExtractor::calculateGreensFunction(Index to, Index from, ChebyshevSolver::GreensFunctionType type){
	vector<Index> toIndices;
	toIndices.push_back(to);

	return calculateGreensFunctions(toIndices, from, type);
}

complex<double>* CPropertyExtractor::calculateGreensFunctions(vector<Index> &to, Index from, ChebyshevSolver::GreensFunctionType type){
	complex<double> *coefficients = new complex<double>[numCoefficients*to.size()];

	if(useGPUToCalculateCoefficients){
		cSolver->calculateCoefficientsGPU(to, from, coefficients, numCoefficients);
	}
	else{
		cSolver->calculateCoefficients(to, from, coefficients, numCoefficients);
//		cout << "Error in PropertyExtractorChebyshev::calculateGreensFunctions: CPU generation of coefficients not yet supported.\n";
//		exit(1);
	}

	complex<double> *greensFunction = new complex<double>[energyResolution*to.size()];

	if(useGPUToGenerateGreensFunctions){
		for(unsigned int n = 0; n < to.size(); n++){
			cSolver->generateGreensFunctionGPU(greensFunction,
								&(coefficients[n*numCoefficients]),
								type);
		}
	}
	else{
		if(useLookupTable){
			#pragma omp parallel for
			for(unsigned int n = 0; n < to.size(); n++){
				cSolver->generateGreensFunction(greensFunction,
								&(coefficients[n*numCoefficients]),
								type);
			}
		}
		else{
			#pragma omp parallel for
			for(unsigned int n = 0; n < to.size(); n++){
				cSolver->generateGreensFunction(greensFunction,
								&(coefficients[n*numCoefficients]),
								numCoefficients,
								energyResolution,
								lowerBound,
								upperBound,
								type);
			}
		}
	}

	delete [] coefficients;

	return greensFunction;
}

complex<double> CPropertyExtractor::calculateExpectationValue(Index to, Index from){
	const complex<double> i(0, 1);

	complex<double> expectationValue = 0.;

	complex<double> *greensFunction = calculateGreensFunction(to, from, ChebyshevSolver::GreensFunctionType::NonPrincipal);
	Model::Statistics statistics = cSolver->getModel()->getStatistics();

	const double dE = (upperBound - lowerBound)/energyResolution;
	for(int e = 0; e < energyResolution; e++){
		double weight;
		if(statistics == Model::Statistics::FermiDirac){
			weight = Functions::fermiDiracDistribution((e/(double)energyResolution - 1.)*(upperBound - lowerBound),
									cSolver->getModel()->getChemicalPotential(),
									cSolver->getModel()->getTemperature());
		}
		else{
			weight = Functions::boseEinsteinDistribution((e/(double)energyResolution - 1.)*(upperBound - lowerBound),
									cSolver->getModel()->getChemicalPotential(),
									cSolver->getModel()->getTemperature());
		}

		expectationValue += weight*conj(i*greensFunction[e])*dE;
	}

	delete [] greensFunction;

	return expectationValue;
}

double* CPropertyExtractor::calculateDensity(Index pattern, Index ranges){
	for(unsigned int n = 0; n < pattern.indices.size(); n++){
		if(pattern.indices.at(n) >= 0)
			ranges.indices.at(n) = 1;
	}

	int densityArraySize = 1;
	for(unsigned int n = 0; n < ranges.indices.size(); n++){
		if(pattern.indices.at(n) < IDX_SUM_ALL)
			densityArraySize *= ranges.indices.at(n);<commit_msg>Added factor 1/M_PI to calculation of expectation value in CPropertyExtractor<commit_after>/** @file CPropertyExtractor.cpp
 *
 *  @author <NAME>
 */

#include "../include/CPropertyExtractor.h"
#include "../include/Functions.h"
#include "../include/Util.h"

using namespace std;

namespace TBTK{

CPropertyExtractor::CPropertyExtractor(ChebyshevSolver *cSolver,
					int numCoefficients,
					int energyResolution,
					bool useGPUToCalculateCoefficients,
					bool useGPUToGenerateGreensFunctions,
					bool useLookupTable,
					double lowerBound,
					double upperBound){
	this->cSolver = cSolver;
	this->numCoefficients = numCoefficients;
	this->energyResolution = energyResolution;
	this->lowerBound = lowerBound;
	this->upperBound = upperBound;
	this->useGPUToCalculateCoefficients = useGPUToCalculateCoefficients;
	this->useGPUToGenerateGreensFunctions = useGPUToGenerateGreensFunctions;
	this->useLookupTable = useLookupTable;

	if(useLookupTable){
		cSolver->generateLookupTable(numCoefficients, energyResolution, lowerBound, upperBound);
		if(useGPUToGenerateGreensFunctions)
			cSolver->loadLookupTableGPU();
	}
	else if(useGPUToGenerateGreensFunctions){
		cout << "Error in PropertyExtractorChebyshev: useLookupTable cannot be false if useGPUToGenerateGreensFunction is true.\n";
		exit(1);
	}
}

CPropertyExtractor::~CPropertyExtractor(){
	if(useGPUToGenerateGreensFunctions)
		cSolver->destroyLookupTableGPU();
}

complex<double>* CPropertyExtractor::calculateGreensFunction(Index to, Index from, ChebyshevSolver::GreensFunctionType type){
	vector<Index> toIndices;
	toIndices.push_back(to);

	return calculateGreensFunctions(toIndices, from, type);
}

complex<double>* CPropertyExtractor::calculateGreensFunctions(vector<Index> &to, Index from, ChebyshevSolver::GreensFunctionType type){
	complex<double> *coefficients = new complex<double>[numCoefficients*to.size()];

	if(useGPUToCalculateCoefficients){
		cSolver->calculateCoefficientsGPU(to, from, coefficients, numCoefficients);
	}
	else{
		cSolver->calculateCoefficients(to, from, coefficients, numCoefficients);
//		cout << "Error in PropertyExtractorChebyshev::calculateGreensFunctions: CPU generation of coefficients not yet supported.\n";
//		exit(1);
	}

	complex<double> *greensFunction = new complex<double>[energyResolution*to.size()];

	if(useGPUToGenerateGreensFunctions){
		for(unsigned int n = 0; n < to.size(); n++){
			cSolver->generateGreensFunctionGPU(greensFunction,
								&(coefficients[n*numCoefficients]),
								type);
		}
	}
	else{
		if(useLookupTable){
			#pragma omp parallel for
			for(unsigned int n = 0; n < to.size(); n++){
				cSolver->generateGreensFunction(greensFunction,
								&(coefficients[n*numCoefficients]),
								type);
			}
		}
		else{
			#pragma omp parallel for
			for(unsigned int n = 0; n < to.size(); n++){
				cSolver->generateGreensFunction(greensFunction,
								&(coefficients[n*numCoefficients]),
								numCoefficients,
								energyResolution,
								lowerBound,
								upperBound,
								type);
			}
		}
	}

	delete [] coefficients;

	return greensFunction;
}

complex<double> CPropertyExtractor::calculateExpectationValue(Index to, Index from){
	const complex<double> i(0, 1);

	complex<double> expectationValue = 0.;

	complex<double> *greensFunction = calculateGreensFunction(to, from, ChebyshevSolver::GreensFunctionType::NonPrincipal);
	Model::Statistics statistics = cSolver->getModel()->getStatistics();

	const double dE = (upperBound - lowerBound)/energyResolution;
	for(int e = 0; e < energyResolution; e++){
		double weight;
		if(statistics == Model::Statistics::FermiDirac){
			weight = Functions::fermiDiracDistribution((e/(double)energyResolution - 1.)*(upperBound - lowerBound),
									cSolver->getModel()->getChemicalPotential(),
									cSolver->getModel()->getTemperature());
		}
		else{
			weight = Functions::boseEinsteinDistribution((e/(double)energyResolution - 1.)*(upperBound - lowerBound),
									cSolver->getModel()->getChemicalPotential(),
									cSolver->getModel()->getTemperature());
		}

		expectationValue += weight*conj(i*greensFunction[e])*dE/M_PI;
	}

	delete [] greensFunction;

	return expectationValue;
}

double* CPropertyExtractor::calculateDensity(Index pattern, Index ranges){
	for(unsigned int n = 0; n < pattern.indices.size(); n++){
		if(pattern.indices.at(n) >= 0)
			ranges.indices.at(n) = 1;
	}

	int densityArraySize = 1;
	for(unsigned int n = 0; n < ranges.indices.size(); n++){
		if(pattern.indices.at(n) < IDX_SUM_ALL)
			densityArraySize *= ranges.indices.at(n);
{'fraction_non_alphanumeric': 0.12468221509892782, 'fraction_numerical': 0.0036476179949154414, 'mean_word_length': 5.283333333333333, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 83, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4148183', 'n_tokens_mistral': 3338, 'n_tokens_neox': 2828, 'n_words': 526}