/[escript]/release/5.2/speckley/test/python/run_readWriteOnSpeckley.py
ViewVC logotype

Contents of /release/5.2/speckley/test/python/run_readWriteOnSpeckley.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6692 - (show annotations)
Mon Jun 25 02:31:06 2018 UTC (3 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 19230 byte(s)
Fix

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import division, print_function
18
19 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 import os
27 import numpy as np
28 import esys.escriptcore.utestselect as unittest
29 from esys.escriptcore.testing import *
30 from esys.escript import *
31 from esys.speckley import *
32
33 try:
34 SPECKLEY_WORKDIR=os.environ['SPECKLEY_WORKDIR']
35 except KeyError:
36 SPECKLEY_WORKDIR='/tmp'
37
38 HAVE_UNZIP = hasFeature('unzip')
39
40 #NE=4 # number elements, must be even
41 #for x in [int(sqrt(mpiSize)),2,3,5,7,1]:
42 # NX=x
43 # NY=mpiSize//x
44 # if NX*NY == mpiSize:
45 # break
46 #
47 #for x in [(int(mpiSize**(1/3.)),int(mpiSize**(1/3.))),(2,3),(2,2),(1,2),(1,1)]:
48 # NXb=x[0]
49 # NYb=x[1]
50 # NZb=mpiSize//(x[0]*x[1])
51 # if NXb*NYb*NZb == mpiSize:
52 # break
53
54 mpiSize = getMPISizeWorld()
55 mpiRank = getMPIRankWorld()
56
57 def adjust(NE, ftype):
58 if ftype == ContinuousFunction:
59 return [i+1 for i in NE]
60 return NE
61
62
63 class WriteBinaryGridTestBase(unittest.TestCase): #subclassing required
64 NX = 4*min(10,mpiSize)
65 NZ = 4
66
67 def generateUniqueData(self, ftype):
68 dim = self.domain.getDim()
69 FSx=ftype(self.domain).getX()
70 NE = adjust(self.NE, ftype)
71 # normalise and scale range of values
72 x = [FSx[i]-inf(FSx[i]) for i in range(dim)]
73 x = [(NE[i]-1)*(x[i]/sup(x[i])) for i in range(dim)]
74 xMax = [int(sup(x[i]))+1 for i in range(dim)]
75 nvals = NE[0]*NE[1]
76 data = x[0] + xMax[0]*x[1]
77 if self.datatype == DATATYPE_INT32:
78 data += 0.05
79 if dim > 2:
80 data += xMax[0]*xMax[1]*x[2]
81 nvals*=NE[2]
82
83 grid = np.array(range(nvals), dtype=self.dtype).reshape(tuple(reversed(NE)))
84 return data, grid
85
86 def writeThenRead(self, data, ftype, fcode):
87 filename = os.path.join(SPECKLEY_WORKDIR, "_wgrid%dd%s"%(self.domain.getDim(),fcode))
88 filename = filename + self.dtype.replace('<','L').replace('>','B')
89 self.domain.writeBinaryGrid(data, filename, self.byteorder, self.datatype)
90 MPIBarrierWorld()
91 result = np.fromfile(filename, dtype=self.dtype).reshape(
92 tuple(reversed(adjust(self.NE,ftype))))
93 return result
94
95 def test_writeGrid2D(self):
96 self.NE = [self.NX*mpiSize, self.NZ]
97 for order in range(2,11):
98 self.domain = Rectangle(order, self.NE[0], self.NE[1], d0=mpiSize)
99 data, ref = self.generateUniqueData(ContinuousFunction)
100 result = self.writeThenRead(data, ContinuousFunction, 'CF')
101 self.assertAlmostEqual(Lsup(ref-result), 0, delta=1e-9,
102 msg="Data doesn't match for "+str(ContinuousFunction(self.domain)))
103
104 def test_writeGrid3D(self):
105 self.NE = [self.NX*mpiSize, self.NX, self.NZ]
106 for order in range(2,11):
107 self.domain = Brick(order, self.NE[0], self.NE[1], self.NE[2], d0=mpiSize)
108 for ftype,fcode in [(ContinuousFunction,'CF')]:
109 data, ref = self.generateUniqueData(ftype)
110 result = self.writeThenRead(data, ftype, fcode)
111 self.assertAlmostEqual(Lsup(ref-result), 0, delta=1e-9,
112 msg="Data doesn't match for "+str(ftype(self.domain)))
113
114 class Test_writeBinaryGridSpeckley_LITTLE_FLOAT32(WriteBinaryGridTestBase):
115 def setUp(self):
116 self.byteorder = BYTEORDER_LITTLE_ENDIAN
117 self.datatype = DATATYPE_FLOAT32
118 self.dtype = "<f4"
119
120 class Test_writeBinaryGridSpeckley_LITTLE_FLOAT64(WriteBinaryGridTestBase):
121 def setUp(self):
122 self.byteorder = BYTEORDER_LITTLE_ENDIAN
123 self.datatype = DATATYPE_FLOAT64
124 self.dtype = "<f8"
125
126 class Test_writeBinaryGridSpeckley_LITTLE_INT32(WriteBinaryGridTestBase):
127 def setUp(self):
128 self.byteorder = BYTEORDER_LITTLE_ENDIAN
129 self.datatype = DATATYPE_INT32
130 self.dtype = "<i4"
131
132 class Test_writeBinaryGridSpeckley_BIG_FLOAT32(WriteBinaryGridTestBase):
133 def setUp(self):
134 self.byteorder = BYTEORDER_BIG_ENDIAN
135 self.datatype = DATATYPE_FLOAT32
136 self.dtype = ">f4"
137
138 class Test_writeBinaryGridSpeckley_BIG_FLOAT64(WriteBinaryGridTestBase):
139 def setUp(self):
140 self.byteorder = BYTEORDER_BIG_ENDIAN
141 self.datatype = DATATYPE_FLOAT64
142 self.dtype = ">f8"
143
144 class Test_writeBinaryGridSpeckley_BIG_INT32(WriteBinaryGridTestBase):
145 def setUp(self):
146 self.byteorder = BYTEORDER_BIG_ENDIAN
147 self.datatype = DATATYPE_INT32
148 self.dtype = ">i4"
149
150
151 class ReadBinaryGridTestBase(unittest.TestCase): #subclassing required
152 """
153 The reader tests work in several stages:
154 1) create numpy array and write to temporary file (ref)
155 2) call readBinaryGrid with that filename
156 3) write the resulting Data object using writeBinaryGrid (test)
157 4) read the result using numpy and compare (ref) and (test)
158 As such, it is important to note that a working writeBinaryGrid() method
159 is assumed!
160 """
161 # set defaults which may be overridden in subclasses
162 NX = 4
163 NZ = 3
164 fspaces = [(ContinuousFunction,'CF')]
165 byteorder = BYTEORDER_NATIVE
166 datatype = DATATYPE_FLOAT64
167 dtype = "f8"
168 shape = ()
169 fill = -42.57
170 first = [0,0,0]
171 multiplier = [1,1,1]
172 reverse = [0,0,0]
173
174 def generateUniqueData(self, ftype):
175 dim = self.domain.getDim()
176 NE = adjust(self.Ndata, ftype)
177 nvals=NE[0]*NE[1]
178 if dim > 2:
179 nvals*=NE[2]
180 grid = np.array(range(nvals), dtype=self.dtype).reshape(tuple(reversed(NE)))
181 return grid
182
183 def write(self, data, filename):
184 self.domain.writeBinaryGrid(data, filename, self.byteorder, self.datatype)
185
186 def read(self, filename, ftype):
187 first = self.first[:self.domain.getDim()]
188 multiplier = self.multiplier[:self.domain.getDim()]
189 reverse = self.reverse[:self.domain.getDim()]
190 numValues=adjust(self.Ndata, ftype)
191 return readBinaryGrid(filename, ftype(self.domain),
192 shape=self.shape, fill=self.fill, byteOrder=self.byteorder,
193 dataType=self.datatype, first=first, numValues=numValues,
194 multiplier=multiplier, reverse=reverse)
195
196 def numpy2Data2Numpy(self, ref, ftype, fcode):
197 filename = os.path.join(SPECKLEY_WORKDIR, "_rgrid%dd%s"%(self.domain.getDim(),fcode))
198 filename = filename + self.dtype.replace('<','L').replace('>','B')
199 if mpiRank == 0:
200 ref.tofile(filename)
201 MPIBarrierWorld()
202 # step 2 - read
203 data = self.read(filename, ftype)
204 MPIBarrierWorld()
205 # step 3 - write
206 self.write(data, filename) # overwrite is ok
207 MPIBarrierWorld()
208 result = np.fromfile(filename, dtype=self.dtype)\
209 .reshape(tuple(reversed(adjust(self.NE,ftype))))
210 return result
211
212 def test_readGrid2D(self):
213 self.NE = [self.NX*mpiSize*self.multiplier[0],
214 self.NZ*self.multiplier[1]]
215 for order in range(2, 11):
216 self.domain = Rectangle(order, self.NE[0], self.NE[1], d0=mpiSize)
217 for ftype,fcode in self.fspaces:
218 self.Ndata = [(self.NX - 1)*mpiSize, self.NZ]
219 if ftype==ContinuousFunction:
220 self.Ndata[1] = self.NZ - 1
221 # step 1 - generate
222 ref = self.generateUniqueData(ftype)
223 # step 2 & 3
224 result = self.numpy2Data2Numpy(ref, ftype, fcode)
225 # apply transformations to be able to compare
226 if self.reverse[0]:
227 result = result[...,::-1]
228 if self.reverse[1]:
229 result = result[::-1,:]
230 for i in range(2):
231 ref = np.repeat(ref, self.multiplier[i], axis=1-i)
232
233 # if domain larger than data: add column(s)/row(s) with fill value
234 fill=np.array(self.fill, dtype=ref.dtype)
235 realNE = adjust(self.NE,ftype)
236 for d in range(2):
237 excess = realNE[d]-ref.shape[1-d]
238 if excess > 0:
239 shape = list(ref.shape)
240 shape[1-d] = excess
241 extra = fill * np.ones(shape)
242 if self.reverse[d]:
243 ref = np.append(extra, ref, axis=1-d)
244 else:
245 ref = np.append(ref, extra, axis=1-d)
246 # step 4 - compare
247 self.assertAlmostEqual(Lsup(ref-result), 0, delta=1e-9,
248 msg="Data doesn't match for "+str(ftype(self.domain))+\
249 "%d"%order)
250
251 def test_readGrid3D(self):
252 self.NE = [self.NX*mpiSize*self.multiplier[0],
253 self.NX*self.multiplier[1], self.NZ*self.multiplier[2]]
254 for order in range(2,11):
255 self.domain = Brick(order,self.NE[0], self.NE[1], self.NE[2],
256 d0=mpiSize)
257 for ftype,fcode in self.fspaces:
258 self.Ndata = [self.NX*mpiSize - 1, self.NX - 1, self.NZ - 1]
259 # step 1 - generate
260 ref = self.generateUniqueData(ftype)
261 # step 2 & 3
262 result = self.numpy2Data2Numpy(ref, ftype, fcode)
263 # apply transformations to be able to compare
264 if self.reverse[0]:
265 result = result[...,::-1]
266 if self.reverse[1]:
267 result = result[...,::-1,:]
268 if self.reverse[2]:
269 result = result[::-1,:,:]
270 for i in range(3):
271 ref = np.repeat(ref, self.multiplier[i], axis=2-i)
272
273 # if domain larger than data: add column(s)/row(s) with fill value
274 fill=np.array(self.fill, dtype=ref.dtype)
275 realNE = adjust(self.NE,ftype)
276 for d in range(3):
277 excess = realNE[d]-ref.shape[2-d]
278 if excess > 0:
279 shape = list(ref.shape)
280 shape[2-d] = excess
281 extra = fill * np.ones(shape)
282 if self.reverse[d]:
283 ref = np.append(extra, ref, axis=2-d)
284 else:
285 ref = np.append(ref, extra, axis=2-d)
286
287 # step 4 - compare
288 self.assertAlmostEqual(Lsup(ref-result), 0, delta=1e-9,
289 msg="Data doesn't match for "+str(ftype(self.domain))+\
290 "%d"%order)
291
292
293 # The following block tests the reader for different byte orders and data
294 # types with domain-filling data (i.e. multiplier=1, reverse=0 and N=NE)
295
296 class Test_readBinaryGridSpeckley_LITTLE_FLOAT32(ReadBinaryGridTestBase):
297 def setUp(self):
298 self.byteorder = BYTEORDER_LITTLE_ENDIAN
299 self.datatype = DATATYPE_FLOAT32
300 self.dtype = "<f4"
301
302 class Test_readBinaryGridSpeckley_LITTLE_FLOAT64(ReadBinaryGridTestBase):
303 def setUp(self):
304 self.byteorder = BYTEORDER_LITTLE_ENDIAN
305 self.datatype = DATATYPE_FLOAT64
306 self.dtype = "<f8"
307
308 #since using getX as the test, doubles required
309 def test_interpolationRead2D(self):
310 self.Ndata = [self.NX*mpiSize, self.NZ]
311 self.NE = [self.NX*mpiSize, self.NZ]
312 for order in range(2,11):
313 self.domain = Rectangle(order, self.NE[0], self.NE[1], d0=mpiSize)
314 X = self.domain.getX()
315 data = X[0] + X[1]
316 filename = os.path.join(SPECKLEY_WORKDIR, "_wgrid%dd%s"%(self.domain.getDim(),'CF'))
317 filename = filename + self.dtype.replace('<','L').replace('>','B')
318 self.write(data, filename)
319 result = self.read(filename, ContinuousFunction)
320 self.assertAlmostEqual(Lsup(data-result), 0, delta=1e-9,
321 msg="Data doesn't match for "+str(ContinuousFunction(self.domain)))
322
323 #since using getX as the test, doubles required
324 def test_interpolationRead3D(self):
325 self.Ndata = [self.NX*mpiSize, self.NX, self.NZ]
326 self.NE = [self.NX*mpiSize, self.NX, self.NZ]
327 for order in range(2,11):
328 self.domain = Brick(order, self.NE[0], self.NE[1], self.NE[2], d0=mpiSize)
329 X = self.domain.getX()
330 data = X[0] + X[1] + X[2]
331 filename = os.path.join(SPECKLEY_WORKDIR, "_wgrid%dd%s"%(self.domain.getDim(),'CF'))
332 filename = filename + self.dtype.replace('<','L').replace('>','B')
333 self.write(data, filename)
334 result = self.read(filename, ContinuousFunction)
335 self.assertAlmostEqual(Lsup(data-result), 0, delta=1e-9,
336 msg="Data doesn't match for "+str(ContinuousFunction(self.domain)))
337
338 class Test_readBinaryGridSpeckley_LITTLE_INT32(ReadBinaryGridTestBase):
339 def setUp(self):
340 self.byteorder = BYTEORDER_LITTLE_ENDIAN
341 self.datatype = DATATYPE_INT32
342 self.dtype = "<i4"
343
344 class Test_readBinaryGridSpeckley_BIG_FLOAT32(ReadBinaryGridTestBase):
345 def setUp(self):
346 self.byteorder = BYTEORDER_BIG_ENDIAN
347 self.datatype = DATATYPE_FLOAT32
348 self.dtype = ">f4"
349
350 class Test_readBinaryGridSpeckley_BIG_FLOAT64(ReadBinaryGridTestBase):
351 def setUp(self):
352 self.byteorder = BYTEORDER_BIG_ENDIAN
353 self.datatype = DATATYPE_FLOAT64
354 self.dtype = ">f8"
355
356 #since using getX as the test, doubles required
357 def test_interpolationRead2D(self):
358 self.Ndata = [self.NX*mpiSize, self.NZ]
359 self.NE = [self.NX*mpiSize, self.NZ]
360 for order in range(2,11):
361 self.domain = Rectangle(order, self.NE[0], self.NE[1], d0=mpiSize)
362 X = self.domain.getX()
363 data = X[0] + X[1]
364 filename = os.path.join(SPECKLEY_WORKDIR, "_wgrid%dd%s"%(self.domain.getDim(),'CF'))
365 filename = filename + self.dtype.replace('<','L').replace('>','B')
366 self.write(data, filename)
367 result = self.read(filename, ContinuousFunction)
368 self.assertAlmostEqual(Lsup(data-result), 0, delta=1e-9,
369 msg="Data doesn't match for "+str(ContinuousFunction(self.domain)))
370
371 #since using getX as the test, doubles required
372 def test_interpolationRead3D(self):
373 self.Ndata = [self.NX*mpiSize, self.NX, self.NZ]
374 self.NE = [self.NX*mpiSize, self.NX, self.NZ]
375 for order in range(2,11):
376 self.domain = Brick(order, self.NE[0], self.NE[1], self.NE[2], d0=mpiSize)
377 X = self.domain.getX()
378 data = X[0] + X[1] + X[2]
379 filename = os.path.join(SPECKLEY_WORKDIR, "_wgrid%dd%s"%(self.domain.getDim(),'CF'))
380 filename = filename + self.dtype.replace('<','L').replace('>','B')
381 self.write(data, filename)
382 result = self.read(filename, ContinuousFunction)
383 self.assertAlmostEqual(Lsup(data-result), 0, delta=1e-9,
384 msg="Data doesn't match for "+str(ContinuousFunction(self.domain)))
385
386
387 class Test_readBinaryGridSpeckley_BIG_INT32(ReadBinaryGridTestBase):
388 def setUp(self):
389 self.byteorder = BYTEORDER_BIG_ENDIAN
390 self.datatype = DATATYPE_INT32
391 self.dtype = ">i4"
392
393 @unittest.skip("reverseX not supported yet")
394 class Test_readBinaryGridSpeckley_reverseX(ReadBinaryGridTestBase):
395 def setUp(self):
396 self.reverse = [1,0,0]
397
398 @unittest.skip("reverseY not supported yet")
399 class Test_readBinaryGridSpeckley_reverseY(ReadBinaryGridTestBase):
400 def setUp(self):
401 self.reverse = [0,1,0]
402
403 class Test_readBinaryGridSpeckley_reverseZ(ReadBinaryGridTestBase):
404 def setUp(self):
405 self.reverse = [0,0,1]
406
407 class Test_readBinaryGridSpeckley_multiplierX(ReadBinaryGridTestBase):
408 def setUp(self):
409 self.multiplier = [2,1,1]
410
411 class Test_readBinaryGridSpeckley_multiplierY(ReadBinaryGridTestBase):
412 def setUp(self):
413 self.multiplier = [1,2,1]
414
415 class Test_readBinaryGridSpeckley_multiplierZ(ReadBinaryGridTestBase):
416 def setUp(self):
417 self.multiplier = [1,1,2]
418
419 class Test_readBinaryGridSpeckley_multiplierXYZ(ReadBinaryGridTestBase):
420 def setUp(self):
421 self.multiplier = [2,3,4]
422
423
424 @unittest.skipIf(getMPISizeWorld() > 1,
425 "Skipping compressed binary grid tests due to element stretching")
426 class Test_readBinaryGridZippedSpeckley(unittest.TestCase):
427 # constants
428 byteorder = BYTEORDER_NATIVE
429 datatype = DATATYPE_FLOAT64
430
431 def read(self, filename, FS, expected, zipped = False):
432 first = [0 for i in expected]
433 reverse = [0 for i in expected]
434 scale = [1 for i in expected]
435
436 if not zipped:
437 return readBinaryGrid(filename, FS, (), 50000,
438 self.byteorder, self.datatype, first, expected, scale, reverse)
439
440 if not HAVE_UNZIP:
441 raise unittest.SkipTest("unzip library not available (boost_iostreams)")
442 return speckleycpp._readBinaryGridFromZipped(filename, FS, (), 50000,
443 self.byteorder, self.datatype, first, expected, scale, reverse)
444
445 def test_readCompressed2D(self):
446 NE = [9, 10]
447 for order in range(2,11):
448 domain = Rectangle(order, NE[0], NE[1], d1=0)
449 for filename, ftype in [("RectConF%s.grid.gz", ContinuousFunction)]:
450 FS = ftype(domain)
451 filename = os.path.join("ref_data", filename%mpiSize)
452 unzipped = self.read(filename[:-3], FS, adjust(NE, ftype))
453 zipped = self.read(filename, FS, adjust(NE, ftype), True)
454 self.assertEqual(Lsup(zipped - unzipped), 0,
455 "Data objects don't match for "+str(FS))
456
457 def test_readCompressed3D(self):
458 NE = [9, 9, 10]
459 for order in range(2,11):
460 domain = Brick(order, NE[0], NE[1], NE[2], d1=0, d2=0)
461 for filename, ftype in [("BrickConF%s.grid.gz", ContinuousFunction)]:
462 FS = ftype(domain)
463 filename = os.path.join("ref_data", filename%mpiSize)
464 unzipped = self.read(filename[:-3], FS, adjust(NE, ftype))
465 zipped = self.read(filename, FS, adjust(NE, ftype), True)
466 self.assertEqual(Lsup(zipped - unzipped), 0,
467 "Data objects don't match for "+str(FS)+"%d"%order)
468
469
470 if __name__ == '__main__':
471 run_tests(__name__, exit_on_failure=True)
472

  ViewVC Help
Powered by ViewVC 1.1.26