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