src/Data/DataAlgorithm.h

// $Id$
/* 
 ******************************************************************************
 *                                                                            *
 *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *
 *                                                                            *
 * This software is the property of ACcESS. No part of this code              *
 * may be copied in any form or by any means without the expressed written    *
 * consent of ACcESS.  Copying, use or modification of this software          *
 * by any unauthorised person is illegal unless that person has a software    *
 * license agreement with ACcESS.                                             *
 *                                                                            *
 ******************************************************************************
*/
                                                                           
#if !defined escript_DataAlgorithm_20040714_H
#define escript_DataAlgorithm_20040714_H

#include "escript/Data/DataExpanded.h"
#include "escript/Data/DataTagged.h"
#include "escript/Data/DataConstant.h"
#include "escript/Data/DataArrayView.h"

#include <iostream>
#include <algorithm>
#include <math.h>
#include <limits>
#include <vector>

namespace escript {

/**
   \brief
   Adapt binary algorithms so they may be used in DataArrayView reduction operations.

   Description:
   This functor adapts the given BinaryFunction operation by starting with the
   given inital value applying this operation to successive values, storing the
   rolling result in m_currentValue - which can be accessed or reset by getResult
   and resetResult respectively.
*/
template <class BinaryFunction>
class DataAlgorithmAdapter {
  public:
    DataAlgorithmAdapter(double initialValue):
      m_initialValue(initialValue),
      m_currentValue(initialValue)
    {
    }
    inline void operator()(double value)
    {
      m_currentValue=operation(m_currentValue,value);
      return;
    }
    inline void resetResult()
    {
      m_currentValue=m_initialValue;
    }
    inline double getResult() const
    {
      return m_currentValue;
    }
  private:
    //
    // the initial operation value
    double m_initialValue;
    //
    // the current operation value
    double m_currentValue;
    //
    // The operation to perform
    BinaryFunction operation;
};

/**
   \brief
   Return the maximum value of the two given values.
*/
struct FMax : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return std::max(x,y);
  }
};

/**
   \brief
   Return the minimum value of the two given values.
*/
struct FMin : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return std::min(x,y);
  }
};

/**
   \brief
   Return the absolute maximum value of the two given values.
*/
struct AbsMax : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return std::max(fabs(x),fabs(y));
  }
};

/**
   \brief
   Return the absolute minimum value of the two given values.
*/
struct AbsMin : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return std::min(fabs(x),fabs(y));
  }
};

/**
   \brief
   Return the length between the two given values.
*/
struct Length : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return std::sqrt(std::pow(x,2)+std::pow(y,2));
  }
};

/**
   \brief
   Return the trace of the two given values.
*/
struct Trace : public std::binary_function<double,double,double>
{
  inline double operator()(double x, double y) const
  {
    return x+y;
  }
};

/**
   \brief
   Perform the given operation upon all values in all data-points in the
   given Data object and return the final result.

   Calls DataArrayView::reductionOp
*/
template <class UnaryFunction>
inline
double
algorithm(DataExpanded& data,
          UnaryFunction operation)
{
  int i,j;
  int numDPPSample=data.getNumDPPSample();
  int numSamples=data.getNumSamples();
  int resultVectorLength=numDPPSample*numSamples;
  std::vector<double> resultVector(resultVectorLength);
  DataArrayView dataView=data.getPointDataView();
  // calculate the reduction operation value for each data point
  // storing the result for each data-point in successive entries
  // in resultVector
  //
  // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
  // which maintains state between calls which would be corrupted by parallel execution
  for (i=0;i<numSamples;i++) {
    for (j=0;j<numDPPSample;j++) {
      resultVector[j*numSamples+i]=dataView.reductionOp(data.getPointOffset(i,j),operation);
    }
  }
  // now calculate the reduction operation value across the results
  // for each data-point
  //
  // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
  // which maintains state between calls which would be corrupted by parallel execution
  operation.resetResult();
  for (int l=0;l<resultVectorLength;l++) {
    operation(resultVector[l]);
  }
  return operation.getResult();
}

template <class UnaryFunction>
inline
double
algorithm(DataTagged& data,
          UnaryFunction operation)
{
  const DataTagged::DataMapType& lookup=data.getTagLookup();
  DataTagged::DataMapType::const_iterator i;
  DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
  DataArrayView& dataView=data.getPointDataView();
  std::vector<double> resultVector;
  int resultVectorLength;
  // perform the operation on each tagged value
  for (i=lookup.begin();i!=lookupEnd;i++) {
    resultVector.push_back(dataView.reductionOp(i->second,operation));
  }
  // perform the operation on the default value
  resultVector.push_back(data.getDefaultValue().reductionOp(operation));
  // now calculate the reduction operation value across the results
  // for each tagged value
  //
  // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
  // which maintains state between calls which would be corrupted by parallel execution
  resultVectorLength=resultVector.size();
  operation.resetResult();
  for (int l=0;l<resultVectorLength;l++) {
    operation(resultVector[l]);
  }
  return operation.getResult();
}

template <class UnaryFunction>
inline
double
algorithm(DataConstant& data,
          UnaryFunction operation)
{
  return data.getPointDataView().reductionOp(operation);
}

/**
   \brief
   Perform the given data-point reduction operation on all data-points
   in data, storing results in corresponding data-points of result.

   Objects data and result must be of the same type, and have the same number
   of data points, but where data has data points of rank n, result must have
   data points of rank 0.

   Calls DataArrayView::reductionOp
*/
template <class UnaryFunction>
inline
void
dp_algorithm(DataExpanded& data,
             DataExpanded& result,
             UnaryFunction operation)
{
  int i,j;
  int numSamples=data.getNumSamples();
  int numDPPSample=data.getNumDPPSample();
  DataArrayView dataView=data.getPointDataView();
  DataArrayView resultView=result.getPointDataView();
  // perform the operation on each data-point and assign
  // this to the corresponding element in result
  //
  // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
  // which maintains state between calls which would be corrupted by parallel execution
  for (i=0;i<numSamples;i++) {
    for (j=0;j<numDPPSample;j++) {
      resultView.getData(result.getPointOffset(i,j)) =
        dataView.reductionOp(data.getPointOffset(i,j),operation);
    }
  }
}

template <class UnaryFunction>
inline
void
dp_algorithm(DataTagged& data,
             DataTagged& result,
             UnaryFunction operation)
{
  const DataTagged::DataMapType& lookup=data.getTagLookup();
  DataTagged::DataMapType::const_iterator i;
  DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
  DataArrayView dataView=data.getPointDataView();
  DataArrayView resultView=result.getPointDataView();
  // perform the operation on each tagged data value
  // and assign this to the corresponding element in result
  //
  // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
  // which maintains state between calls which would be corrupted by parallel execution
  for (i=lookup.begin();i!=lookupEnd;i++) {
    resultView.getData(i->second) =
      dataView.reductionOp(i->second,operation);
  }
  // perform the operation on the default data value
  // and assign this to the default element in result
  resultView.getData(0) =
    data.getDefaultValue().reductionOp(operation);
}

template <class UnaryFunction>
inline
void
dp_algorithm(DataConstant& data,
             DataConstant& result,
             UnaryFunction operation)
{
  // perform the operation on the data value
  // and assign this to the element in result
  result.getPointDataView().getData(0) =
    data.getPointDataView().reductionOp(operation);
}

} // end of namespace

#endif
1	jgs	102	// $Id$
2	jgs	82	/*
3			******************************************************************************
4			* *
5			* COPYRIGHT ACcESS 2004 - All Rights Reserved *
6			* *
7			* This software is the property of ACcESS. No part of this code *
8			* may be copied in any form or by any means without the expressed written *
9			* consent of ACcESS. Copying, use or modification of this software *
10			* by any unauthorised person is illegal unless that person has a software *
11			* license agreement with ACcESS. *
12			* *
13			******************************************************************************
14			*/
15
16	jgs	102	#if !defined escript_DataAlgorithm_20040714_H
17	jgs	82	#define escript_DataAlgorithm_20040714_H
18
19			#include "escript/Data/DataExpanded.h"
20			#include "escript/Data/DataTagged.h"
21			#include "escript/Data/DataConstant.h"
22			#include "escript/Data/DataArrayView.h"
23
24			#include <iostream>
25			#include <algorithm>
26			#include <math.h>
27			#include <limits>
28	jgs	113	#include <vector>
29	jgs	82
30			namespace escript {
31	jgs	102
32	jgs	82	/**
33			\brief
34	jgs	113	Adapt binary algorithms so they may be used in DataArrayView reduction operations.
35	jgs	82
36			Description:
37	jgs	113	This functor adapts the given BinaryFunction operation by starting with the
38			given inital value applying this operation to successive values, storing the
39			rolling result in m_currentValue - which can be accessed or reset by getResult
40			and resetResult respectively.
41	jgs	82	*/
42	jgs	113	template <class BinaryFunction>
43			class DataAlgorithmAdapter {
44			public:
45			DataAlgorithmAdapter(double initialValue):
46			m_initialValue(initialValue),
47			m_currentValue(initialValue)
48			{
49			}
50			inline void operator()(double value)
51			{
52			m_currentValue=operation(m_currentValue,value);
53			return;
54			}
55			inline void resetResult()
56			{
57			m_currentValue=m_initialValue;
58			}
59			inline double getResult() const
60			{
61			return m_currentValue;
62			}
63			private:
64			//
65			// the initial operation value
66			double m_initialValue;
67			//
68			// the current operation value
69			double m_currentValue;
70			//
71			// The operation to perform
72			BinaryFunction operation;
73			};
74
75			/**
76			\brief
77			Return the maximum value of the two given values.
78			*/
79	jgs	82	struct FMax : public std::binary_function<double,double,double>
80			{
81			inline double operator()(double x, double y) const
82			{
83			return std::max(x,y);
84			}
85			};
86
87			/**
88			\brief
89	jgs	113	Return the minimum value of the two given values.
90	jgs	82	*/
91			struct FMin : public std::binary_function<double,double,double>
92			{
93			inline double operator()(double x, double y) const
94			{
95			return std::min(x,y);
96			}
97			};
98
99			/**
100			\brief
101	jgs	113	Return the absolute maximum value of the two given values.
102	jgs	82	*/
103			struct AbsMax : public std::binary_function<double,double,double>
104			{
105			inline double operator()(double x, double y) const
106			{
107			return std::max(fabs(x),fabs(y));
108			}
109			};
110	jgs	102
111	jgs	82	/**
112			\brief
113	jgs	117	Return the absolute minimum value of the two given values.
114			*/
115			struct AbsMin : public std::binary_function<double,double,double>
116			{
117			inline double operator()(double x, double y) const
118			{
119			return std::min(fabs(x),fabs(y));
120			}
121			};
122
123			/**
124			\brief
125	jgs	113	Return the length between the two given values.
126	jgs	106	*/
127			struct Length : public std::binary_function<double,double,double>
128			{
129			inline double operator()(double x, double y) const
130			{
131			return std::sqrt(std::pow(x,2)+std::pow(y,2));
132			}
133			};
134
135			/**
136			\brief
137	jgs	113	Return the trace of the two given values.
138	jgs	106	*/
139			struct Trace : public std::binary_function<double,double,double>
140			{
141			inline double operator()(double x, double y) const
142			{
143			return x+y;
144			}
145			};
146
147			/**
148			\brief
149	jgs	113	Perform the given operation upon all values in all data-points in the
150			given Data object and return the final result.
151	jgs	82
152	jgs	113	Calls DataArrayView::reductionOp
153	jgs	82	*/
154			template <class UnaryFunction>
155	jgs	102	inline
156			double
157			algorithm(DataExpanded& data,
158			UnaryFunction operation)
159	jgs	82	{
160			int i,j;
161	jgs	113	int numDPPSample=data.getNumDPPSample();
162			int numSamples=data.getNumSamples();
163			int resultVectorLength=numDPPSample*numSamples;
164			std::vector<double> resultVector(resultVectorLength);
165			DataArrayView dataView=data.getPointDataView();
166			// calculate the reduction operation value for each data point
167			// storing the result for each data-point in successive entries
168			// in resultVector
169	jgs	122	//
170			// this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
171			// which maintains state between calls which would be corrupted by parallel execution
172			for (i=0;i<numSamples;i++) {
173			for (j=0;j<numDPPSample;j++) {
174			resultVector[j*numSamples+i]=dataView.reductionOp(data.getPointOffset(i,j),operation);
175	jgs	82	}
176			}
177	jgs	113	// now calculate the reduction operation value across the results
178			// for each data-point
179	jgs	122	//
180			// this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
181			// which maintains state between calls which would be corrupted by parallel execution
182	jgs	113	operation.resetResult();
183			for (int l=0;l<resultVectorLength;l++) {
184			operation(resultVector[l]);
185			}
186	jgs	82	return operation.getResult();
187			}
188
189			template <class UnaryFunction>
190	jgs	102	inline
191			double
192			algorithm(DataTagged& data,
193			UnaryFunction operation)
194	jgs	82	{
195			const DataTagged::DataMapType& lookup=data.getTagLookup();
196			DataTagged::DataMapType::const_iterator i;
197			DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
198			DataArrayView& dataView=data.getPointDataView();
199	jgs	113	std::vector<double> resultVector;
200			int resultVectorLength;
201			// perform the operation on each tagged value
202	jgs	102	for (i=lookup.begin();i!=lookupEnd;i++) {
203	jgs	113	resultVector.push_back(dataView.reductionOp(i->second,operation));
204	jgs	82	}
205	jgs	113	// perform the operation on the default value
206			resultVector.push_back(data.getDefaultValue().reductionOp(operation));
207			// now calculate the reduction operation value across the results
208			// for each tagged value
209	jgs	122	//
210			// this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
211			// which maintains state between calls which would be corrupted by parallel execution
212	jgs	113	resultVectorLength=resultVector.size();
213			operation.resetResult();
214			for (int l=0;l<resultVectorLength;l++) {
215			operation(resultVector[l]);
216			}
217	jgs	82	return operation.getResult();
218			}
219
220			template <class UnaryFunction>
221	jgs	102	inline
222			double
223			algorithm(DataConstant& data,
224			UnaryFunction operation)
225	jgs	82	{
226	jgs	108	return data.getPointDataView().reductionOp(operation);
227	jgs	82	}
228
229	jgs	102	/**
230			\brief
231	jgs	113	Perform the given data-point reduction operation on all data-points
232			in data, storing results in corresponding data-points of result.
233	jgs	82
234	jgs	102	Objects data and result must be of the same type, and have the same number
235	jgs	106	of data points, but where data has data points of rank n, result must have
236			data points of rank 0.
237	jgs	102
238	jgs	113	Calls DataArrayView::reductionOp
239	jgs	102	*/
240			template <class UnaryFunction>
241			inline
242			void
243	jgs	106	dp_algorithm(DataExpanded& data,
244			DataExpanded& result,
245	jgs	102	UnaryFunction operation)
246			{
247			int i,j;
248	jgs	121	int numSamples=data.getNumSamples();
249	jgs	113	int numDPPSample=data.getNumDPPSample();
250			DataArrayView dataView=data.getPointDataView();
251			DataArrayView resultView=result.getPointDataView();
252			// perform the operation on each data-point and assign
253			// this to the corresponding element in result
254	jgs	122	//
255			// this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
256			// which maintains state between calls which would be corrupted by parallel execution
257			for (i=0;i<numSamples;i++) {
258			for (j=0;j<numDPPSample;j++) {
259			resultView.getData(result.getPointOffset(i,j)) =
260			dataView.reductionOp(data.getPointOffset(i,j),operation);
261	jgs	102	}
262			}
263			}
264
265			template <class UnaryFunction>
266			inline
267			void
268	jgs	106	dp_algorithm(DataTagged& data,
269			DataTagged& result,
270	jgs	102	UnaryFunction operation)
271			{
272			const DataTagged::DataMapType& lookup=data.getTagLookup();
273			DataTagged::DataMapType::const_iterator i;
274			DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
275	jgs	113	DataArrayView dataView=data.getPointDataView();
276			DataArrayView resultView=result.getPointDataView();
277			// perform the operation on each tagged data value
278			// and assign this to the corresponding element in result
279	jgs	122	//
280			// this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
281			// which maintains state between calls which would be corrupted by parallel execution
282	jgs	102	for (i=lookup.begin();i!=lookupEnd;i++) {
283	jgs	113	resultView.getData(i->second) =
284			dataView.reductionOp(i->second,operation);
285	jgs	102	}
286	jgs	113	// perform the operation on the default data value
287	jgs	106	// and assign this to the default element in result
288	jgs	113	resultView.getData(0) =
289			data.getDefaultValue().reductionOp(operation);
290	jgs	102	}
291
292			template <class UnaryFunction>
293			inline
294			void
295	jgs	106	dp_algorithm(DataConstant& data,
296			DataConstant& result,
297	jgs	102	UnaryFunction operation)
298			{
299	jgs	113	// perform the operation on the data value
300			// and assign this to the element in result
301	jgs	106	result.getPointDataView().getData(0) =
302	jgs	113	data.getPointDataView().reductionOp(operation);
303	jgs	102	}
304
305	jgs	82	} // end of namespace
306	jgs	122
307	jgs	82	#endif
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision