ExpressionEvaluator.cp

#include "ExpressionEvaluator.h"
#include "CubicData.h"
#include "TrilinearInterpolation.h"
#include "DefinedField.h"
#include "QuantumProgress.h"

#include <cassert>
#include <ctime>
#include <iostream>
using namespace std;

vecN3 ExpressionEvaluator::globalDefaultResolution = {32,32,32};

ExpressionEvaluator::ExpressionEvaluator(GeneralCubicData *data)
	: buffered(false), notOwned(false), theDefinedField(NULL)
{
	if(data)
	{
		buffered = notOwned = true;
		bufferedData = data;
		bufferedData->IncRefCount();
	}
	xRange = yRange = zRange = fullRange;
}

ExpressionEvaluator::ExpressionEvaluator(DefinedField *defined)
	: theDefinedField(defined)
{
	buffered = notOwned = true;
	bufferedData = defined->GetData();
	bufferedData->IncRefCount();
	defined->AddEvaluator(this);
	defined->GetRanges(xRange,yRange,zRange);
}

ExpressionEvaluator::~ExpressionEvaluator()
{
	if(buffered)
		bufferedData->DecRefCount();
	if(theDefinedField)
		theDefinedField->RemoveEvaluator(this);
}

bool ExpressionEvaluator::AlwaysCreateField()
{
	return false;
}

bool ExpressionEvaluator::IsComplex()
{
	if(buffered)
		return bufferedData->C != NULL;
	else
		return false;			
}

int ExpressionEvaluator::Dimensions()
{
	if(buffered)
		return bufferedData->C ? bufferedData->C->n : bufferedData->R->n;
	else
		return 1;
}

bool ExpressionEvaluator::GetDefaultResolution(vecN3& res)
{
	if(buffered)
	{
		if(bufferedData->C)
		{
			res.x = bufferedData->C->cx;
			res.y = bufferedData->C->cy;
			res.z = bufferedData->C->cz;
		}
		else
		{
			res.x = bufferedData->R->cx;
			res.y = bufferedData->R->cy;
			res.z = bufferedData->R->cz;
		}
		return true;
	}
	else
		return false;
	//else
		//res = globalDefaultResolution; // should we ever return this?
	//return true;
}

void ExpressionEvaluator::SetGlobalDefaultResolution(vecN3 res)
{
	globalDefaultResolution = res;
}

void ExpressionEvaluator::CheckTypeMismatch(bool complex, int n)
{
	if(complex != IsComplex() || n != Dimensions())
		throw TypeMismatch();
}

vecR3 ExpressionEvaluator::MapPos(vecR3 pos)
{
	pos.x = -1 + 2 * (pos.x - xRange.first)/(xRange.second - xRange.first);
	pos.y = -1 + 2 * (pos.y - yRange.first)/(yRange.second - yRange.first);
	pos.z = -1 + 2 * (pos.z - zRange.first)/(zRange.second - zRange.first);
	return pos;
}
	
vecN3 intPos;	//###
	
void ExpressionEvaluator::EvaluateReal(vecR3 pos, int n, number values[])
{
	if((buffered || AlwaysCreateField()) && !IsComplex())
	{
		PreBuffer();
		if(n != bufferedData->R->n)
			throw TypeMismatch();
		pos = MapPos(pos);
		for(int i=0;i<n;i++)
		{
			values[i] = TrilinearInterpolation(bufferedData->R, pos, i);
		}
	}
	else
		throw TypeMismatch();
}


void ExpressionEvaluator::EvaluateComplex(vecR3 pos, int n, complex<number> values[])
{
	if((buffered || AlwaysCreateField()) && IsComplex())
	{
		PreBuffer();
		if(n != bufferedData->C->n)
			throw TypeMismatch();
		pos = MapPos(pos);
		for(int i=0;i<n;i++)
		{
			values[i] = TrilinearInterpolation(bufferedData->C, pos, i);
		}
	}
	else
	{
		number tmp[maxDimensions];
		
		EvaluateReal(pos,n,tmp);
		for(int i=0;i<n;i++)
			values[i] = tmp[i];
	}
}

vecN3 ExpressionEvaluator::GetRequestedResolution(vecN3 res,bool forceRes)
{
	if(res.x == 0)
	{
		if(!GetDefaultResolution(res))
		{
			//res.x = res.y = res.z = 32;
			res = globalDefaultResolution;
				// ### 
		}
	}
	else if(!forceRes)
		GetDefaultResolution(res);
		
	return res;
}


GeneralCubicData* ExpressionEvaluator::DoCreateField(vecN3 res, bool forceRes)
{
	res = GetRequestedResolution(res,forceRes);

	StartDeterminateProgress();
	cout << "evaluating field at " 
			<< res.x << "x" << res.y << "x" << res.z << ".\n";
	cout << "[                  ]" << endl;
	clock_t startTime = clock();
	
	const int progressMax = 20;
	int progressDone = 0;
		
	bool c = IsComplex();
	int n = Dimensions();

	number xmin = xRange.first;
	number xsiz = xRange.size();
	number ymin = yRange.first;
	number ysiz = yRange.size();
	number zmin = zRange.first;
	number zsiz = zRange.size();

		
	GeneralCubicData *d = new GeneralCubicData();
	if(c)
		d->C = new CubicData< complex<number> > (res.x, res.y, res.z, n);
	else
		d->R = new CubicData< number > (res.x, res.y, res.z, n);
	
	number theMin = 0, theMax = 0;
	bool	first = true;
		
	try
	{
		for(int i=0;i<res.x;i++)
		{
			for(int j=0;j<res.y;j++)
				for(int k=0;k<res.z;k++)
				{
					intPos.x = i, intPos.y = j, intPos.z = k;
					vecR3 pos;
					pos.x = xmin + xsiz*number(i)/(res.x-1);
					pos.y = ymin + ysiz*number(j)/(res.y-1);
					pos.z = zmin + zsiz*number(k)/(res.z-1);
					if(c)
						EvaluateComplex(pos, n, &d->C->data[i][j][n*k]);
					else
						EvaluateReal(pos, n, &d->R->data[i][j][n*k]);
												
					if(!c && n == 1)
					{
						number val = d->R->data[i][j][k];
						if(first)
						{
							theMin = theMax = val;
							first = false;
						}
						else
						{
							if(val < theMin)
								theMin = val;
							if(val > theMax)
								theMax = val;
						}
					}
				}
			int progressShould = progressMax*(i+1)/res.x;
			while(progressDone < progressShould)
				cout << '*', progressDone++;
			cout << flush;
			Progress(double(i+1)/res.x);
			CheckAbort();
		}
	}
	catch(AbortException)
	{
		d->DecRefCount();
		throw;
	}
	cout << endl;
	d->UpdateInfo(theMin, theMax);

	cout << "evaluation done in "
		 << double(clock() - startTime)/CLOCKS_PER_SEC
		 << " seconds.\n";
	EndDeterminateProgress();

	return d;
}

void ExpressionEvaluator::PreBuffer()
{
	if(buffered)
		return;
	bufferedData = DoCreateField();
	buffered = true;
}

GeneralCubicData* ExpressionEvaluator::CreateField(vecN3 res, bool forceRes)
{	
	if(buffered)
	{
		res = GetRequestedResolution(res,forceRes);
		
		bool resMatch = true;
		if(bufferedData->C)
			resMatch = (bufferedData->C->cx == res.x 
						&& bufferedData->C->cy == res.y
						&& bufferedData->C->cz == res.z);
		else
			resMatch = (bufferedData->R->cx == res.x 
						&& bufferedData->R->cy == res.y
						&& bufferedData->R->cz == res.z);
		if(resMatch)
		{
			bufferedData->IncRefCount();
			return bufferedData;
		}
	}
	
	return DoCreateField(res, forceRes);
}

void ExpressionEvaluator::OrphanDefinedEvaluator()
{
	theDefinedField = NULL;
}
void ExpressionEvaluator::UpdateDefinedEvaluator(GeneralCubicData *data)
{
	assert(bufferedData);
	bufferedData->DecRefCount();
	bufferedData = data;
	bufferedData->IncRefCount();
}

void ExpressionEvaluator::GetRanges(Range& x, Range& y, Range& z)
{
	x = xRange;
	y = yRange;
	z = zRange;
}