1111import numpy as np
1212
1313
14+ def chunk_to_slice (chunk ):
15+ """
16+ Convert an openPMD_api.ChunkInfo to np.s_
17+ """
18+ #print("chunk={} {}".format(chunk, type(chunk)))
19+ stops = [a + b for a , b in zip (chunk .offset , chunk .extent )]
20+ indices_per_dim = zip (chunk .offset , stops )
21+ index_tuple = map (lambda s : slice (s [0 ], s [1 ], None ), indices_per_dim )
22+ return tuple (index_tuple )
23+
1424def get_data (series , record_component , i_slice = None , pos_slice = None ,
1525 output_type = np .float64 ):
1626 """
@@ -46,8 +56,17 @@ def get_data(series, record_component, i_slice=None, pos_slice=None,
4656 if i_slice is not None and not isinstance (i_slice , list ):
4757 i_slice = [i_slice ]
4858
59+ chunks = record_component .available_chunks ()
60+
4961 if pos_slice is None :
50- data = record_component [()]
62+ # mask invalid regions with zero
63+ data = np .zeros_like (record_component )
64+ for chunk in chunks :
65+ chunk_slice = chunk_to_slice (chunk )
66+ # read only valid region
67+ x = record_component [chunk_slice ]
68+ series .flush ()
69+ data [chunk_slice ] = x
5170 else :
5271 # Get largest element of pos_slice
5372 max_pos = max (pos_slice )
@@ -60,8 +79,33 @@ def get_data(series, record_component, i_slice=None, pos_slice=None,
6079 list_index [dir_index ] = i_slice [count ]
6180 # Convert list_index into a tuple
6281 tuple_index = tuple (list_index )
63- # Slice dset according to tuple_index
64- data = record_component [tuple_index ]
82+ print ("tuple_index={}" .format (tuple_index ))
83+
84+ # potentially a better approach as below, since we only slice
85+ # out hyperplanes, planes & lines:
86+ # - allocate zero array for result, which is a hyperplane/plane/line
87+ # - iterate over slices in tuple_index
88+ # - reduce selected read range to "valid" range
89+
90+ # initial experiment:
91+ # full_indices can be HUGE, avoid!!
92+ full_indices = np .indices (record_component .shape )[0 ]
93+ #full_shape = full_indices.shape
94+ #print("full_shape.shape={}".format(full_shape))
95+ #print("full_shape={}".format(full_shape))
96+
97+ # prepare sliced data according to tuple_index
98+ slice_indices = full_indices [tuple_index ]
99+ slice_shape = slice_indices .shape
100+ data = np .zeros (slice_shape , dtype = output_type )
101+ # write now in index space between intersection of slice_indices and chunk indices
102+ for chunk in chunks :
103+ chunk_slice = chunk_to_slice (chunk )
104+ chunk_indices = full_indices [chunk_slice ]
105+ intersect_indices = np .intersect1d (chunk_indices , slice_indices )
106+ print (intersect_indices )
107+ data [slice_indices ] = record_component [intersect_indices ]
108+ #data = np.zeros_like(record_component)[tuple_index] # just avoid invalid reads for now
65109
66110 series .flush ()
67111
0 commit comments