1 |
// $Id$ |
2 |
/* |
3 |
************************************************************ |
4 |
* Copyright 2006 by ACcESS MNRF * |
5 |
* * |
6 |
* http://www.access.edu.au * |
7 |
* Primary Business: Queensland, Australia * |
8 |
* Licensed under the Open Software License version 3.0 * |
9 |
* http://www.opensource.org/licenses/osl-3.0.php * |
10 |
* * |
11 |
************************************************************ |
12 |
*/ |
13 |
|
14 |
#include "DataExpanded.h" |
15 |
#include "DataException.h" |
16 |
#include "DataConstant.h" |
17 |
#include "DataTagged.h" |
18 |
|
19 |
#include <boost/python/extract.hpp> |
20 |
|
21 |
using namespace std; |
22 |
using namespace boost::python; |
23 |
using namespace boost; |
24 |
|
25 |
namespace escript { |
26 |
|
27 |
DataExpanded::DataExpanded(const boost::python::numeric::array& value, |
28 |
const FunctionSpace& what) |
29 |
: DataAbstract(what) |
30 |
{ |
31 |
DataArrayView::ShapeType tempShape; |
32 |
// |
33 |
// extract the shape of the python numarray |
34 |
for (int i=0; i<value.getrank(); i++) { |
35 |
tempShape.push_back(extract<int>(value.getshape()[i])); |
36 |
} |
37 |
// |
38 |
// initialise the data array for this object |
39 |
initialise(tempShape,what.getNumSamples(),what.getNumDPPSample()); |
40 |
// |
41 |
// copy the given value to every data point |
42 |
copy(value); |
43 |
} |
44 |
|
45 |
DataExpanded::DataExpanded(const DataExpanded& other) |
46 |
: DataAbstract(other.getFunctionSpace()), |
47 |
m_data(other.m_data) |
48 |
{ |
49 |
// |
50 |
// create the view for the data |
51 |
DataArrayView temp(m_data.getData(),other.getPointDataView().getShape()); |
52 |
setPointDataView(temp); |
53 |
} |
54 |
|
55 |
DataExpanded::DataExpanded(const DataConstant& other) |
56 |
: DataAbstract(other.getFunctionSpace()) |
57 |
{ |
58 |
// |
59 |
// initialise the data array for this object |
60 |
initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample()); |
61 |
// |
62 |
// DataConstant only has one value, copy this to every data point |
63 |
copy(other.getPointDataView()); |
64 |
} |
65 |
|
66 |
DataExpanded::DataExpanded(const DataTagged& other) |
67 |
: DataAbstract(other.getFunctionSpace()) |
68 |
{ |
69 |
// |
70 |
// initialise the data array for this object |
71 |
initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample()); |
72 |
// |
73 |
// for each data point in this object, extract and copy the corresponding data |
74 |
// value from the given DataTagged object |
75 |
int i,j; |
76 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
77 |
DataArrayView::ValueType::size_type numCols=m_data.getNumCols(); |
78 |
#pragma omp parallel for private(i,j) schedule(static) |
79 |
for (i=0;i<numRows;i++) { |
80 |
for (j=0;j<numCols;j++) { |
81 |
try { |
82 |
getPointDataView().copy(getPointOffset(i,j), |
83 |
other.getPointDataView(), |
84 |
other.getPointOffset(i,j)); |
85 |
} |
86 |
catch (std::exception& e) { |
87 |
cout << e.what() << endl; |
88 |
} |
89 |
} |
90 |
} |
91 |
} |
92 |
|
93 |
DataExpanded::DataExpanded(const DataExpanded& other, |
94 |
const DataArrayView::RegionType& region) |
95 |
: DataAbstract(other.getFunctionSpace()) |
96 |
{ |
97 |
// |
98 |
// get the shape of the slice |
99 |
DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region)); |
100 |
// |
101 |
// initialise this Data object to the shape of the slice |
102 |
initialise(shape,other.getNumSamples(),other.getNumDPPSample()); |
103 |
// |
104 |
// copy the data |
105 |
DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region); |
106 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
107 |
DataArrayView::ValueType::size_type numCols=m_data.getNumCols(); |
108 |
int i,j; |
109 |
#pragma omp parallel for private(i,j) schedule(static) |
110 |
for (i=0;i<numRows;i++) { |
111 |
for (j=0;j<numCols;j++) { |
112 |
try { |
113 |
getPointDataView().copySlice(getPointOffset(i,j), |
114 |
other.getPointDataView(), |
115 |
other.getPointOffset(i,j), |
116 |
region_loop_range); |
117 |
} |
118 |
catch (std::exception& e) { |
119 |
cout << e.what() << endl; |
120 |
} |
121 |
} |
122 |
} |
123 |
} |
124 |
|
125 |
DataExpanded::DataExpanded(const DataArrayView& value, |
126 |
const FunctionSpace& what) |
127 |
: DataAbstract(what) |
128 |
{ |
129 |
// |
130 |
// get the shape of the given data value |
131 |
DataArrayView::ShapeType tempShape=value.getShape(); |
132 |
// |
133 |
// initialise this Data object to the shape of the given data value |
134 |
initialise(tempShape,what.getNumSamples(),what.getNumDPPSample()); |
135 |
// |
136 |
// copy the given value to every data point |
137 |
copy(value); |
138 |
} |
139 |
|
140 |
DataExpanded::DataExpanded(const FunctionSpace& what, |
141 |
const DataArrayView::ShapeType &shape, |
142 |
const DataArrayView::ValueType &data) |
143 |
: DataAbstract(what) |
144 |
{ |
145 |
// |
146 |
// create the view of the data |
147 |
initialise(shape,what.getNumSamples(),what.getNumDPPSample()); |
148 |
// |
149 |
// copy the data in the correct format |
150 |
m_data.getData()=data; |
151 |
} |
152 |
|
153 |
DataExpanded::~DataExpanded() |
154 |
{ |
155 |
} |
156 |
|
157 |
void |
158 |
DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape) |
159 |
{ |
160 |
if (getPointDataView().getRank()!=0) { |
161 |
stringstream temp; |
162 |
temp << "Error - Can only reshape Data with data points of rank 0. " |
163 |
<< "This Data has data points with rank: " |
164 |
<< getPointDataView().getRank(); |
165 |
throw DataException(temp.str()); |
166 |
} |
167 |
// |
168 |
// create the new DataBlocks2D data array, and a corresponding DataArrayView |
169 |
DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape)); |
170 |
DataArrayView newView(newData.getData(),shape); |
171 |
// |
172 |
// Copy the original data to every value for the new shape |
173 |
int i,j; |
174 |
int nRows=m_data.getNumRows(); |
175 |
int nCols=m_data.getNumCols(); |
176 |
#pragma omp parallel for private(i,j) schedule(static) |
177 |
for (i=0;i<nRows;i++) { |
178 |
for (j=0;j<nCols;j++) { |
179 |
// NOTE: An exception may be thown from this call if |
180 |
// DOASSERT is on which of course will play |
181 |
// havoc with the omp threads. Run single threaded |
182 |
// if using DOASSERT. |
183 |
newView.copy(newData.index(i,j),m_data(i,j)); |
184 |
} |
185 |
} |
186 |
// swap the new data array for the original |
187 |
m_data.Swap(newData); |
188 |
// set the corresponding DataArrayView |
189 |
DataArrayView temp(m_data.getData(),shape); |
190 |
setPointDataView(temp); |
191 |
} |
192 |
|
193 |
DataAbstract* |
194 |
DataExpanded::getSlice(const DataArrayView::RegionType& region) const |
195 |
{ |
196 |
return new DataExpanded(*this,region); |
197 |
} |
198 |
|
199 |
void |
200 |
DataExpanded::setSlice(const DataAbstract* value, |
201 |
const DataArrayView::RegionType& region) |
202 |
{ |
203 |
const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value); |
204 |
if (tempDataExp==0) { |
205 |
throw DataException("Programming error - casting to DataExpanded."); |
206 |
} |
207 |
// |
208 |
// get shape of slice |
209 |
DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region)); |
210 |
DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region); |
211 |
// |
212 |
// check shape |
213 |
if (getPointDataView().getRank()!=region.size()) { |
214 |
throw DataException("Error - Invalid slice region."); |
215 |
} |
216 |
if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) { |
217 |
throw DataException (value->getPointDataView().createShapeErrorMessage( |
218 |
"Error - Couldn't copy slice due to shape mismatch.",shape)); |
219 |
} |
220 |
// |
221 |
// copy the data from the slice into this object |
222 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
223 |
DataArrayView::ValueType::size_type numCols=m_data.getNumCols(); |
224 |
int i, j; |
225 |
#pragma omp parallel for private(i,j) schedule(static) |
226 |
for (i=0;i<numRows;i++) { |
227 |
for (j=0;j<numCols;j++) { |
228 |
getPointDataView().copySliceFrom(getPointOffset(i,j), |
229 |
tempDataExp->getPointDataView(), |
230 |
tempDataExp->getPointOffset(i,j), |
231 |
region_loop_range); |
232 |
} |
233 |
} |
234 |
} |
235 |
|
236 |
void |
237 |
DataExpanded::copy(const DataArrayView& value) |
238 |
{ |
239 |
// |
240 |
// copy a single value to every data point in this object |
241 |
int nRows=m_data.getNumRows(); |
242 |
int nCols=m_data.getNumCols(); |
243 |
int i,j; |
244 |
#pragma omp parallel for private(i,j) schedule(static) |
245 |
for (i=0;i<nRows;i++) { |
246 |
for (j=0;j<nCols;j++) { |
247 |
// NOTE: An exception may be thown from this call if |
248 |
// DOASSERT is on which of course will play |
249 |
// havoc with the omp threads. Run single threaded |
250 |
// if using DOASSERT. |
251 |
getPointDataView().copy(m_data.index(i,j),value); |
252 |
} |
253 |
} |
254 |
} |
255 |
|
256 |
void |
257 |
DataExpanded::copy(const boost::python::numeric::array& value) |
258 |
{ |
259 |
// |
260 |
// first convert the numarray into a DataArray object |
261 |
DataArray temp(value); |
262 |
// |
263 |
// check the input shape matches this shape |
264 |
if (!getPointDataView().checkShape(temp.getView().getShape())) { |
265 |
throw DataException(getPointDataView().createShapeErrorMessage( |
266 |
"Error - (DataExpanded) Cannot copy due to shape mismatch.", |
267 |
temp.getView().getShape())); |
268 |
} |
269 |
// |
270 |
// now copy over the data |
271 |
copy(temp.getView()); |
272 |
} |
273 |
|
274 |
void |
275 |
DataExpanded::initialise(const DataArrayView::ShapeType& shape, |
276 |
int noSamples, |
277 |
int noDataPointsPerSample) |
278 |
{ |
279 |
// |
280 |
// resize data array to the required size |
281 |
m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape)); |
282 |
// |
283 |
// create the data view of the data array |
284 |
DataArrayView temp(m_data.getData(),shape); |
285 |
setPointDataView(temp); |
286 |
} |
287 |
|
288 |
string |
289 |
DataExpanded::toString() const |
290 |
{ |
291 |
stringstream temp; |
292 |
// |
293 |
// create a temporary view as the offset will be changed |
294 |
DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset()); |
295 |
for (int i=0;i<m_data.getNumRows();i++) { |
296 |
for (int j=0;j<m_data.getNumCols();j++) { |
297 |
tempView.setOffset(m_data.index(i,j)); |
298 |
stringstream suffix; |
299 |
suffix << "(" << i << "," << j << ")"; |
300 |
temp << tempView.toString(suffix.str()); |
301 |
if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) { |
302 |
temp << endl; |
303 |
} |
304 |
} |
305 |
} |
306 |
return temp.str(); |
307 |
} |
308 |
|
309 |
DataArrayView::ValueType::size_type |
310 |
DataExpanded::getPointOffset(int sampleNo, |
311 |
int dataPointNo) const |
312 |
{ |
313 |
return m_data.index(sampleNo,dataPointNo); |
314 |
} |
315 |
|
316 |
DataArrayView |
317 |
DataExpanded::getDataPoint(int sampleNo, |
318 |
int dataPointNo) |
319 |
{ |
320 |
DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo)); |
321 |
return temp; |
322 |
} |
323 |
|
324 |
DataArrayView::ValueType::size_type |
325 |
DataExpanded::getLength() const |
326 |
{ |
327 |
return m_data.size(); |
328 |
} |
329 |
|
330 |
void |
331 |
DataExpanded::setRefValue(int ref, |
332 |
const DataArray& value) |
333 |
{ |
334 |
// |
335 |
// Get the number of samples and data-points per sample. |
336 |
int numSamples = getNumSamples(); |
337 |
int numDPPSample = getNumDPPSample(); |
338 |
|
339 |
// |
340 |
// Determine the sample number which corresponds to this reference number. |
341 |
int sampleNo = -1; |
342 |
int tempRef = -1; |
343 |
for (int n=0; n<numSamples; n++) { |
344 |
tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n); |
345 |
if (tempRef == ref) { |
346 |
sampleNo = n; |
347 |
break; |
348 |
} |
349 |
} |
350 |
if (sampleNo == -1) { |
351 |
throw DataException("DataExpanded::setRefValue error: invalid ref number supplied."); |
352 |
} |
353 |
|
354 |
for (int n=0; n<numDPPSample; n++) { |
355 |
// |
356 |
// Get *each* data-point in the sample in turn. |
357 |
DataArrayView pointView = getDataPoint(sampleNo, n); |
358 |
// |
359 |
// Assign the values in the DataArray to this data-point. |
360 |
pointView.copy(value.getView()); |
361 |
} |
362 |
} |
363 |
|
364 |
void |
365 |
DataExpanded::getRefValue(int ref, |
366 |
DataArray& value) |
367 |
{ |
368 |
// |
369 |
// Get the number of samples and data-points per sample. |
370 |
int numSamples = getNumSamples(); |
371 |
int numDPPSample = getNumDPPSample(); |
372 |
|
373 |
// |
374 |
// Determine the sample number which corresponds to this reference number |
375 |
int sampleNo = -1; |
376 |
int tempRef = -1; |
377 |
for (int n=0; n<numSamples; n++) { |
378 |
tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n); |
379 |
if (tempRef == ref) { |
380 |
sampleNo = n; |
381 |
break; |
382 |
} |
383 |
} |
384 |
if (sampleNo == -1) { |
385 |
throw DataException("DataExpanded::getRefValue error: invalid ref number supplied."); |
386 |
} |
387 |
|
388 |
// |
389 |
// Get the *first* data-point associated with this sample number. |
390 |
DataArrayView pointView = getDataPoint(sampleNo, 0); |
391 |
|
392 |
// |
393 |
// Load the values from this data-point into the DataArray |
394 |
value.getView().copy(pointView); |
395 |
} |
396 |
|
397 |
int |
398 |
DataExpanded::archiveData(ofstream& archiveFile, |
399 |
const DataArrayView::ValueType::size_type noValues) const |
400 |
{ |
401 |
return(m_data.archiveData(archiveFile, noValues)); |
402 |
} |
403 |
|
404 |
int |
405 |
DataExpanded::extractData(ifstream& archiveFile, |
406 |
const DataArrayView::ValueType::size_type noValues) |
407 |
{ |
408 |
return(m_data.extractData(archiveFile, noValues)); |
409 |
} |
410 |
|
411 |
void |
412 |
DataExpanded::copyAll(const boost::python::numeric::array& value) { |
413 |
// |
414 |
// Get the number of samples and data-points per sample. |
415 |
int numSamples = getNumSamples(); |
416 |
int numDataPointsPerSample = getNumDPPSample(); |
417 |
int dataPointRank = getPointDataView().getRank(); |
418 |
ShapeType dataPointShape = getPointDataView().getShape(); |
419 |
// |
420 |
// check rank: |
421 |
if (value.getrank()!=dataPointRank+1) |
422 |
throw DataException("Rank of numarray does not match Data object rank"); |
423 |
if (value.getshape()[0]!=numSamples*numDataPointsPerSample) |
424 |
throw DataException("leading dimension of numarray is too small"); |
425 |
// |
426 |
int dataPoint = 0; |
427 |
for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
428 |
for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
429 |
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
430 |
if (dataPointRank==0) { |
431 |
dataPointView()=extract<double>(value[dataPoint]); |
432 |
} else if (dataPointRank==1) { |
433 |
for (int i=0; i<dataPointShape[0]; i++) { |
434 |
dataPointView(i)=extract<double>(value[dataPoint][i]); |
435 |
} |
436 |
} else if (dataPointRank==2) { |
437 |
for (int i=0; i<dataPointShape[0]; i++) { |
438 |
for (int j=0; j<dataPointShape[1]; j++) { |
439 |
dataPointView(i,j)=extract<double>(value[dataPoint][i][j]); |
440 |
} |
441 |
} |
442 |
} else if (dataPointRank==3) { |
443 |
for (int i=0; i<dataPointShape[0]; i++) { |
444 |
for (int j=0; j<dataPointShape[1]; j++) { |
445 |
for (int k=0; k<dataPointShape[2]; k++) { |
446 |
dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]); |
447 |
} |
448 |
} |
449 |
} |
450 |
} else if (dataPointRank==4) { |
451 |
for (int i=0; i<dataPointShape[0]; i++) { |
452 |
for (int j=0; j<dataPointShape[1]; j++) { |
453 |
for (int k=0; k<dataPointShape[2]; k++) { |
454 |
for (int l=0; l<dataPointShape[3]; l++) { |
455 |
dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]); |
456 |
} |
457 |
} |
458 |
} |
459 |
} |
460 |
} |
461 |
dataPoint++; |
462 |
} |
463 |
} |
464 |
} |
465 |
void |
466 |
DataExpanded::eigenvalues(DataAbstract* ev) |
467 |
{ |
468 |
int sampleNo,dataPointNo; |
469 |
int numSamples = getNumSamples(); |
470 |
int numDataPointsPerSample = getNumDPPSample(); |
471 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
472 |
if (temp_ev==0) { |
473 |
throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error)."); |
474 |
} |
475 |
DataArrayView& thisView=getPointDataView(); |
476 |
DataArrayView& evView=ev->getPointDataView(); |
477 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
478 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
479 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
480 |
DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo), |
481 |
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
482 |
} |
483 |
} |
484 |
} |
485 |
void |
486 |
DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol) |
487 |
{ |
488 |
int numSamples = getNumSamples(); |
489 |
int numDataPointsPerSample = getNumDPPSample(); |
490 |
int sampleNo,dataPointNo; |
491 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
492 |
if (temp_ev==0) { |
493 |
throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error)."); |
494 |
} |
495 |
DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V); |
496 |
if (temp_V==0) { |
497 |
throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error)."); |
498 |
} |
499 |
DataArrayView& thisView=getPointDataView(); |
500 |
DataArrayView& evView=ev->getPointDataView(); |
501 |
DataArrayView& VView=V->getPointDataView(); |
502 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
503 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
504 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
505 |
DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo), |
506 |
evView,ev->getPointOffset(sampleNo,dataPointNo), |
507 |
VView,V->getPointOffset(sampleNo,dataPointNo), |
508 |
tol); |
509 |
} |
510 |
} |
511 |
} |
512 |
|
513 |
} // end of namespace |