escript/src/DataLazy.cpp


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


#include "DataLazy.h"
#ifdef USE_NETCDF
#include <netcdfcpp.h>
#endif
#ifdef PASO_MPI
#include <mpi.h>
#endif
#include "FunctionSpace.h"
#include "DataTypes.h"
#include "Data.h"

using namespace std;
using namespace boost;

namespace escript
{

const std::string&
opToString(ES_optype op);

namespace
{

// return the FunctionSpace of the result of "left op right"
FunctionSpace
resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
{
    // perhaps this should call interpolate and throw or something?
    // maybe we need an interpolate node -
    // that way, if interpolate is required in any other op we can just throw a 
    // programming error exception.


    if (left->getFunctionSpace()!=right->getFunctionSpace())
    {
        throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
    }
    return left->getFunctionSpace();
}

// return the shape of the result of "left op right"
DataTypes::ShapeType
resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
{
    if (left->getShape()!=right->getShape())
    {
        throw DataException("Shapes not the same - shapes must match for lazy data.");
    }
    return left->getShape();
}

size_t
resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
{
   switch(op)
   {
//   case IDENTITY: return left->getLength();
   case ADD:    // the length is preserved in these ops
   case SUB:
   case MUL:
   case DIV: return left->getLength();
   default: 
    throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");

   }
}

int
calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
{
   switch(op)
   {
   case IDENTITY: return 0;
   case ADD:    // the length is preserved in these ops
   case SUB:
   case MUL:
   case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());
   default: 
    throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
   }
}

string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};
int ES_opcount=5;

// void performOp(ValueType& v, int startOffset, ES_optype op, int m_samplesize)
// {
//    switch(op)
//    {
//    case ADD:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
//      startOffset+m_samplesize,plus<double>());
//        break;    
//    case SUB:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
//      startOffset+m_samplesize,minus<double>());
//        break;
//    case MUL:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
//      startOffset+m_samplesize,multiplies<double>());
//        break;
//    case DIV:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
//      startOffset+m_samplesize,divides<double>());
//        break;
//    default: 
//  throw DataException("Programmer Error - attempt to performOp() for operator "+opToString(op)+".");
//    }
// 
// }

}   // end anonymous namespace


const std::string&
opToString(ES_optype op)
{
  if (op<0 || op>=ES_opcount) 
  {
    op=UNKNOWNOP;
  }
  return ES_opstrings[op];
}


DataLazy::DataLazy(DataAbstract_ptr p)
    : parent(p->getFunctionSpace(),p->getShape()),
    m_op(IDENTITY)
{
   if (p->isLazy())
   {
    // TODO: fix this.   We could make the new node a copy of p?
    // I don't want identity of Lazy.
    // Question: Why would that be so bad?
    // Answer: We assume that the child of ID is something we can call getVector on
    throw DataException("Programmer error - attempt to create identity from a DataLazy.");
   }
   else
   {
    m_id=dynamic_pointer_cast<DataReady>(p);
   }
   m_length=p->getLength();
   m_buffsRequired=0;
   m_samplesize=getNumDPPSample()*getNoValues();
cout << "(1)Lazy created with " << m_samplesize << endl;
}

DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
    : parent(resultFS(left,right,op), resultShape(left,right,op)),
    m_left(left),
    m_right(right),
    m_op(op)
{
   m_length=resultLength(m_left,m_right,m_op);
   m_samplesize=getNumDPPSample()*getNoValues();
   m_buffsRequired=calcBuffs(m_left, m_right, m_op);
cout << "(2)Lazy created with " << m_samplesize << endl;
}

DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
    : parent(resultFS(left,right,op), resultShape(left,right,op)),
    m_op(op)
{
   if (left->isLazy())
   {
    m_left=dynamic_pointer_cast<DataLazy>(left);
   }
   else
   {
    m_left=DataLazy_ptr(new DataLazy(left));
   }
   if (right->isLazy())
   {
    m_right=dynamic_pointer_cast<DataLazy>(right);
   }
   else
   {
    m_right=DataLazy_ptr(new DataLazy(right));
   }

   m_length=resultLength(m_left,m_right,m_op);
   m_samplesize=getNumDPPSample()*getNoValues();
   m_buffsRequired=calcBuffs(m_left, m_right,m_op);
cout << "(3)Lazy created with " << m_samplesize << endl;
}


DataLazy::~DataLazy()
{
}


int
DataLazy::getBuffsRequired() const
{
    return m_buffsRequired;
}

void
DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const
{
  if (m_op==IDENTITY)   // copy the contents into the vector
  {
cout << "Begin ID" << endl;
cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;
    const ValueType& vec=m_id->getVector();
    size_t srcOffset=m_id->getPointOffset(sampleNo, 0);
cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;
    for (size_t i=0;i<m_samplesize;++i,++srcOffset,++offset)
    {
cout << "Trying offset=" << offset << " srcOffset=" << srcOffset << endl;
    v[offset]=vec[srcOffset];   
    }
cout << "End ID" << endl;
    return;
  }
  size_t rightoffset=offset+m_samplesize;
  m_left->resolveSample(v,sampleNo,offset);
  m_right->resolveSample(v,sampleNo,rightoffset);
//  for (int i=0;i<getNumDPPSample();++i)
  {
    switch(m_op)
    {
    case ADD:       // since these are pointwise ops, pretend each sample is one point
    tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),plus<double>());
    break;
    case SUB:       
    tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),minus<double>());
    break;
    case MUL:       
    tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),multiplies<double>());
    break;
    case DIV:       
    tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),divides<double>());
    break;
    default:
    throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
    }
  }
}

DataReady_ptr
DataLazy::resolve()
{
  // This is broken!     We need to have a buffer per thread!
  // so the allocation of v needs to move inside the loop somehow

cout << "Sample size=" << m_samplesize << endl;
cout << "Buffers=" << m_buffsRequired << endl;

  ValueType v(m_samplesize*max(1,m_buffsRequired));
cout << "Buffer created with size=" << v.size() << endl;
  ValueType dummy(getNoValues());
  DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
  ValueType& resvec=result->getVector();
  DataReady_ptr resptr=DataReady_ptr(result);
  int sample;
  #pragma omp parallel for private(sample) schedule(static)
  for (sample=0;sample<getNumSamples();++sample)
  {
    resolveSample(v,sample,0);
    for (int i=0;i<m_samplesize;++i)    // copy values into the output vector
    {
    resvec[i]=v[i];
    }
  }
  return resptr;
}

std::string
DataLazy::toString() const
{
  return "Lazy evaluation object. No details available.";
}

// Note that in this case, deepCopy does not make copies of the leaves.
// Hopefully copy on write (or whatever we end up using) will take care of this.
DataAbstract* 
DataLazy::deepCopy()
{
  if (m_op==IDENTITY)
  {
    return new DataLazy(m_left);    // we don't need to copy the child here
  }
  return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op); 
}


DataTypes::ValueType::size_type
DataLazy::getLength() const
{
  return m_length;
}


DataAbstract*
DataLazy::getSlice(const DataTypes::RegionType& region) const
{
  throw DataException("getSlice - not implemented for Lazy objects.");
}

DataTypes::ValueType::size_type 
DataLazy::getPointOffset(int sampleNo,
                 int dataPointNo) const
{
  throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
}

// // The startOffset is where to write results in the output vector v
// void
// DataLazy::processSample(ValueType& v, int sampleNo, size_t startOffset)
// {
//     m_left.processSample(v,sampleNo,startOffset);
//     m_right.processSample(v,sampleNo,startOffset+getSampleSize());
//     int i;
//     #pragma omp parallel for private(i) schedule(static)
//     for (i=0;i<getSampleSize();++i)
//     {
//      performOp(v,startOffset+i*m_pointsize,ES_optype,m_samplesize);
//     }
// }


}   // end namespace
1	jfenwick	1865
2			/*******************************************************
3			*
4			* Copyright (c) 2003-2008 by University of Queensland
5			* Earth Systems Science Computational Center (ESSCC)
6			* http://www.uq.edu.au/esscc
7			*
8			* Primary Business: Queensland, Australia
9			* Licensed under the Open Software License version 3.0
10			* http://www.opensource.org/licenses/osl-3.0.php
11			*
12			*******************************************************/
13
14
15			#include "DataLazy.h"
16			#ifdef USE_NETCDF
17			#include <netcdfcpp.h>
18			#endif
19			#ifdef PASO_MPI
20			#include <mpi.h>
21			#endif
22			#include "FunctionSpace.h"
23			#include "DataTypes.h"
24	jfenwick	1868	#include "Data.h"
25	jfenwick	1865
26			using namespace std;
27	jfenwick	1868	using namespace boost;
28	jfenwick	1865
29			namespace escript
30			{
31
32			const std::string&
33			opToString(ES_optype op);
34
35			namespace
36			{
37
38			// return the FunctionSpace of the result of "left op right"
39			FunctionSpace
40			resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
41			{
42			// perhaps this should call interpolate and throw or something?
43			// maybe we need an interpolate node -
44			// that way, if interpolate is required in any other op we can just throw a
45			// programming error exception.
46	jfenwick	1879
47
48			if (left->getFunctionSpace()!=right->getFunctionSpace())
49			{
50			throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
51			}
52	jfenwick	1865	return left->getFunctionSpace();
53			}
54
55			// return the shape of the result of "left op right"
56			DataTypes::ShapeType
57			resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
58			{
59	jfenwick	1879	if (left->getShape()!=right->getShape())
60			{
61			throw DataException("Shapes not the same - shapes must match for lazy data.");
62			}
63			return left->getShape();
64	jfenwick	1865	}
65
66			size_t
67			resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
68			{
69			switch(op)
70			{
71	jfenwick	1879	// case IDENTITY: return left->getLength();
72	jfenwick	1868	case ADD: // the length is preserved in these ops
73			case SUB:
74			case MUL:
75			case DIV: return left->getLength();
76	jfenwick	1865	default:
77			throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
78
79			}
80			}
81
82	jfenwick	1879	int
83			calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
84			{
85			switch(op)
86			{
87			case IDENTITY: return 0;
88			case ADD: // the length is preserved in these ops
89			case SUB:
90			case MUL:
91			case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());
92			default:
93			throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
94			}
95			}
96
97	jfenwick	1868	string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};
98			int ES_opcount=5;
99	jfenwick	1865
100	jfenwick	1879	// void performOp(ValueType& v, int startOffset, ES_optype op, int m_samplesize)
101			// {
102			// switch(op)
103			// {
104			// case ADD: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
105			// startOffset+m_samplesize,plus<double>());
106			// break;
107			// case SUB: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
108			// startOffset+m_samplesize,minus<double>());
109			// break;
110			// case MUL: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
111			// startOffset+m_samplesize,multiplies<double>());
112			// break;
113			// case DIV: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
114			// startOffset+m_samplesize,divides<double>());
115			// break;
116			// default:
117			// throw DataException("Programmer Error - attempt to performOp() for operator "+opToString(op)+".");
118			// }
119			//
120			// }
121
122	jfenwick	1865	} // end anonymous namespace
123
124
125			const std::string&
126			opToString(ES_optype op)
127			{
128			if (op<0 \|\| op>=ES_opcount)
129			{
130			op=UNKNOWNOP;
131			}
132			return ES_opstrings[op];
133			}
134
135
136			DataLazy::DataLazy(DataAbstract_ptr p)
137			: parent(p->getFunctionSpace(),p->getShape()),
138			m_op(IDENTITY)
139			{
140	jfenwick	1879	if (p->isLazy())
141			{
142			// TODO: fix this. We could make the new node a copy of p?
143			// I don't want identity of Lazy.
144			// Question: Why would that be so bad?
145			// Answer: We assume that the child of ID is something we can call getVector on
146			throw DataException("Programmer error - attempt to create identity from a DataLazy.");
147			}
148			else
149			{
150			m_id=dynamic_pointer_cast<DataReady>(p);
151			}
152			m_length=p->getLength();
153			m_buffsRequired=0;
154			m_samplesize=getNumDPPSample()*getNoValues();
155			cout << "(1)Lazy created with " << m_samplesize << endl;
156	jfenwick	1865	}
157
158	jfenwick	1879	DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
159	jfenwick	1865	: parent(resultFS(left,right,op), resultShape(left,right,op)),
160			m_left(left),
161			m_right(right),
162			m_op(op)
163			{
164	jfenwick	1879	m_length=resultLength(m_left,m_right,m_op);
165			m_samplesize=getNumDPPSample()*getNoValues();
166			m_buffsRequired=calcBuffs(m_left, m_right, m_op);
167			cout << "(2)Lazy created with " << m_samplesize << endl;
168	jfenwick	1865	}
169
170	jfenwick	1879	DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
171			: parent(resultFS(left,right,op), resultShape(left,right,op)),
172			m_op(op)
173			{
174			if (left->isLazy())
175			{
176			m_left=dynamic_pointer_cast<DataLazy>(left);
177			}
178			else
179			{
180			m_left=DataLazy_ptr(new DataLazy(left));
181			}
182			if (right->isLazy())
183			{
184			m_right=dynamic_pointer_cast<DataLazy>(right);
185			}
186			else
187			{
188			m_right=DataLazy_ptr(new DataLazy(right));
189			}
190
191			m_length=resultLength(m_left,m_right,m_op);
192			m_samplesize=getNumDPPSample()*getNoValues();
193			m_buffsRequired=calcBuffs(m_left, m_right,m_op);
194			cout << "(3)Lazy created with " << m_samplesize << endl;
195			}
196
197
198	jfenwick	1865	DataLazy::~DataLazy()
199			{
200			}
201
202	jfenwick	1879
203			int
204			DataLazy::getBuffsRequired() const
205	jfenwick	1865	{
206	jfenwick	1879	return m_buffsRequired;
207			}
208
209			void
210			DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
211			{
212			if (m_op==IDENTITY) // copy the contents into the vector
213	jfenwick	1865	{
214	jfenwick	1879	cout << "Begin ID" << endl;
215			cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;
216			const ValueType& vec=m_id->getVector();
217			size_t srcOffset=m_id->getPointOffset(sampleNo, 0);
218			cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;
219			for (size_t i=0;i<m_samplesize;++i,++srcOffset,++offset)
220			{
221			cout << "Trying offset=" << offset << " srcOffset=" << srcOffset << endl;
222			v[offset]=vec[srcOffset];
223			}
224			cout << "End ID" << endl;
225			return;
226	jfenwick	1865	}
227	jfenwick	1879	size_t rightoffset=offset+m_samplesize;
228			m_left->resolveSample(v,sampleNo,offset);
229			m_right->resolveSample(v,sampleNo,rightoffset);
230			// for (int i=0;i<getNumDPPSample();++i)
231	jfenwick	1865	{
232	jfenwick	1879	switch(m_op)
233			{
234			case ADD: // since these are pointwise ops, pretend each sample is one point
235			tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),plus<double>());
236			break;
237			case SUB:
238			tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),minus<double>());
239			break;
240			case MUL:
241			tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),multiplies<double>());
242			break;
243			case DIV:
244			tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),divides<double>());
245			break;
246	jfenwick	1868	default:
247			throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
248	jfenwick	1879	}
249	jfenwick	1868	}
250	jfenwick	1865	}
251
252	jfenwick	1879	DataReady_ptr
253			DataLazy::resolve()
254			{
255			// This is broken! We need to have a buffer per thread!
256			// so the allocation of v needs to move inside the loop somehow
257
258			cout << "Sample size=" << m_samplesize << endl;
259			cout << "Buffers=" << m_buffsRequired << endl;
260
261			ValueType v(m_samplesize*max(1,m_buffsRequired));
262			cout << "Buffer created with size=" << v.size() << endl;
263			ValueType dummy(getNoValues());
264			DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
265			ValueType& resvec=result->getVector();
266			DataReady_ptr resptr=DataReady_ptr(result);
267			int sample;
268			#pragma omp parallel for private(sample) schedule(static)
269			for (sample=0;sample<getNumSamples();++sample)
270			{
271			resolveSample(v,sample,0);
272			for (int i=0;i<m_samplesize;++i) // copy values into the output vector
273			{
274			resvec[i]=v[i];
275			}
276			}
277			return resptr;
278			}
279
280	jfenwick	1865	std::string
281			DataLazy::toString() const
282			{
283			return "Lazy evaluation object. No details available.";
284			}
285
286			// Note that in this case, deepCopy does not make copies of the leaves.
287			// Hopefully copy on write (or whatever we end up using) will take care of this.
288			DataAbstract*
289			DataLazy::deepCopy()
290			{
291			if (m_op==IDENTITY)
292			{
293			return new DataLazy(m_left); // we don't need to copy the child here
294			}
295			return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
296			}
297
298
299			DataTypes::ValueType::size_type
300			DataLazy::getLength() const
301			{
302	jfenwick	1879	return m_length;
303	jfenwick	1865	}
304
305
306			DataAbstract*
307			DataLazy::getSlice(const DataTypes::RegionType& region) const
308			{
309	jfenwick	1879	throw DataException("getSlice - not implemented for Lazy objects.");
310	jfenwick	1865	}
311
312			DataTypes::ValueType::size_type
313			DataLazy::getPointOffset(int sampleNo,
314			int dataPointNo) const
315			{
316			throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
317			}
318
319	jfenwick	1879	// // The startOffset is where to write results in the output vector v
320			// void
321			// DataLazy::processSample(ValueType& v, int sampleNo, size_t startOffset)
322			// {
323			// m_left.processSample(v,sampleNo,startOffset);
324			// m_right.processSample(v,sampleNo,startOffset+getSampleSize());
325			// int i;
326			// #pragma omp parallel for private(i) schedule(static)
327			// for (i=0;i<getSampleSize();++i)
328			// {
329			// performOp(v,startOffset+i*m_pointsize,ES_optype,m_samplesize);
330			// }
331			// }
332
333
334	jfenwick	1865	} // end namespace