1import functools 

2import hashlib 

3import operator 

4import pickle 

5import re 

6from enum import Enum 

7from itertools import chain 

8from typing import List, Optional, Sequence, Set, Tuple 

9 

10import numpy as np 

11import sympy as sp 

12from sympy.core.cache import cacheit 

13 

14import pystencils 

15from pystencils.alignedarray import aligned_empty 

16from pystencils.data_types import StructType, TypedSymbol, create_type 

17from pystencils.kernelparameters import FieldShapeSymbol, FieldStrideSymbol 

18from pystencils.stencil import ( 

19 direction_string_to_offset, inverse_direction, offset_to_direction_string) 

20from pystencils.sympyextensions import is_integer_sequence 

21 

22__all__ = ['Field', 'fields', 'FieldType', 'AbstractField'] 

23 

24 

25class FieldType(Enum): 

26 # generic fields 

27 GENERIC = 0 

28 # index fields are currently only used for boundary handling 

29 # the coordinates are not the loop counters in that case, but are read from this index field 

30 INDEXED = 1 

31 # communication buffer, used for (un)packing data in communication. 

32 BUFFER = 2 

33 # unsafe fields may be accessed in an absolute fashion - the index depends on the data 

34 # and thus may lead to out-of-bounds accesses 

35 CUSTOM = 3 

36 # staggered field 

37 STAGGERED = 4 

38 # staggered field that reverses sign when accessed via opposite direction 

39 STAGGERED_FLUX = 5 

40 

41 @staticmethod 

42 def is_generic(field): 

43 assert isinstance(field, Field) 

44 return field.field_type == FieldType.GENERIC 

45 

46 @staticmethod 

47 def is_indexed(field): 

48 assert isinstance(field, Field) 

49 return field.field_type == FieldType.INDEXED 

50 

51 @staticmethod 

52 def is_buffer(field): 

53 assert isinstance(field, Field) 

54 return field.field_type == FieldType.BUFFER 

55 

56 @staticmethod 

57 def is_custom(field): 

58 assert isinstance(field, Field) 

59 return field.field_type == FieldType.CUSTOM 

60 

61 @staticmethod 

62 def is_staggered(field): 

63 assert isinstance(field, Field) 

64 return field.field_type == FieldType.STAGGERED or field.field_type == FieldType.STAGGERED_FLUX 

65 

66 @staticmethod 

67 def is_staggered_flux(field): 

68 assert isinstance(field, Field) 

69 return field.field_type == FieldType.STAGGERED_FLUX 

70 

71 

72def fields(description=None, index_dimensions=0, layout=None, field_type=FieldType.GENERIC, **kwargs): 

73 """Creates pystencils fields from a string description. 

74 

75 Examples: 

76 Create a 2D scalar and vector field: 

77 >>> s, v = fields("s, v(2): double[2D]") 

78 >>> assert s.spatial_dimensions == 2 and s.index_dimensions == 0 

79 >>> assert (v.spatial_dimensions, v.index_dimensions, v.index_shape) == (2, 1, (2,)) 

80 

81 Create an integer field of shape (10, 20): 

82 >>> f = fields("f : int32[10, 20]") 

83 >>> f.has_fixed_shape, f.shape 

84 (True, (10, 20)) 

85 

86 Numpy arrays can be used as template for shape and data type of field: 

87 >>> arr_s, arr_v = np.zeros([20, 20]), np.zeros([20, 20, 2]) 

88 >>> s, v = fields("s, v(2)", s=arr_s, v=arr_v) 

89 >>> assert s.index_dimensions == 0 and s.dtype.numpy_dtype == arr_s.dtype 

90 >>> assert v.index_shape == (2,) 

91 

92 Format string can be left out, field names are taken from keyword arguments. 

93 >>> fields(f1=arr_s, f2=arr_s) 

94 [f1: double[20,20], f2: double[20,20]] 

95 

96 The keyword names ``index_dimension`` and ``layout`` have special meaning, don't use them for field names 

97 >>> f = fields(f=arr_v, index_dimensions=1) 

98 >>> assert f.index_dimensions == 1 

99 >>> f = fields("pdfs(19) : float32[3D]", layout='fzyx') 

100 >>> f.layout 

101 (2, 1, 0) 

102 """ 

103 result = [] 

104 if description: 

105 field_descriptions, dtype, shape = _parse_description(description) 

106 layout = 'numpy' if layout is None else layout 

107 for field_name, idx_shape in field_descriptions: 

108 if field_name in kwargs: 

109 arr = kwargs[field_name] 

110 idx_shape_of_arr = () if not len(idx_shape) else arr.shape[-len(idx_shape):] 

111 assert idx_shape_of_arr == idx_shape 

112 f = Field.create_from_numpy_array(field_name, kwargs[field_name], index_dimensions=len(idx_shape), 

113 field_type=field_type) 

114 elif isinstance(shape, tuple): 114 ↛ 115line 114 didn't jump to line 115, because the condition on line 114 was never true

115 f = Field.create_fixed_size(field_name, shape + idx_shape, dtype=dtype, 

116 index_dimensions=len(idx_shape), layout=layout, field_type=field_type) 

117 elif isinstance(shape, int): 117 ↛ 120line 117 didn't jump to line 120, because the condition on line 117 was never false

118 f = Field.create_generic(field_name, spatial_dimensions=shape, dtype=dtype, 

119 index_shape=idx_shape, layout=layout, field_type=field_type) 

120 elif shape is None: 

121 f = Field.create_generic(field_name, spatial_dimensions=2, dtype=dtype, 

122 index_shape=idx_shape, layout=layout, field_type=field_type) 

123 else: 

124 assert False 

125 result.append(f) 

126 else: 

127 assert layout is None, "Layout can not be specified when creating Field from numpy array" 

128 for field_name, arr in kwargs.items(): 

129 result.append(Field.create_from_numpy_array(field_name, arr, index_dimensions=index_dimensions, 

130 field_type=field_type)) 

131 

132 if len(result) == 0: 132 ↛ 133line 132 didn't jump to line 133, because the condition on line 132 was never true

133 return None 

134 elif len(result) == 1: 

135 return result[0] 

136 else: 

137 return result 

138 

139 

140class AbstractField: 

141 

142 class AbstractAccess: 

143 pass 

144 

145 

146class Field(AbstractField): 

147 """ 

148 With fields one can formulate stencil-like update rules on structured grids. 

149 This Field class knows about the dimension, memory layout (strides) and optionally about the size of an array. 

150 

151 Creating Fields: 

152 The preferred method to create fields is the `fields` function. 

153 Alternatively one can use one of the static functions `Field.create_generic`, `Field.create_from_numpy_array` 

154 and `Field.create_fixed_size`. Don't instantiate the Field directly! 

155 Fields can be created with known or unknown shapes: 

156 

157 1. If you want to create a kernel with fixed loop sizes i.e. the shape of the array is already known. 

158 This is usually the case if just-in-time compilation directly from Python is done. 

159 (see `Field.create_from_numpy_array` 

160 2. create a more general kernel that works for variable array sizes. This can be used to create kernels 

161 beforehand for a library. (see `Field.create_generic`) 

162 

163 Dimensions and Indexing: 

164 A field has spatial and index dimensions, where the spatial dimensions come first. 

165 The interpretation is that the field has multiple cells in (usually) two or three dimensional space which are 

166 looped over. Additionally N values are stored per cell. In this case spatial_dimensions is two or three, 

167 and index_dimensions equals N. If you want to store a matrix on each point in a two dimensional grid, there 

168 are four dimensions, two spatial and two index dimensions: ``len(arr.shape) == spatial_dims + index_dims`` 

169 

170 The shape of the index dimension does not have to be specified. Just use the 'index_dimensions' parameter. 

171 However, it is good practice to define the shape, since out of bounds accesses can be directly detected in this 

172 case. The shape can be passed with the 'index_shape' parameter of the field creation functions. 

173 

174 When accessing (indexing) a field the result is a `Field.Access` which is derived from sympy Symbol. 

175 First specify the spatial offsets in [], then in case index_dimension>0 the indices in () 

176 e.g. ``f[-1,0,0](7)`` 

177 

178 Staggered Fields: 

179 Staggered fields are used to store a value on a second grid shifted by half a cell with respect to the usual 

180 grid. 

181 

182 The first index dimension is used to specify the position on the staggered grid (e.g. 0 means half-way to the 

183 eastern neighbor, 1 is half-way to the northern neighbor, etc.), while additional indices can be used to store 

184 multiple values at each position. 

185 

186 Example using no index dimensions: 

187 >>> a = np.zeros([10, 10]) 

188 >>> f = Field.create_from_numpy_array("f", a, index_dimensions=0) 

189 >>> jacobi = (f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) / 4 

190 

191 Examples for index dimensions to create LB field and implement stream pull: 

192 >>> from pystencils import Assignment 

193 >>> stencil = np.array([[0,0], [0,1], [0,-1]]) 

194 >>> src, dst = fields("src(3), dst(3) : double[2D]") 

195 >>> assignments = [Assignment(dst[0,0](i), src[-offset](i)) for i, offset in enumerate(stencil)]; 

196 """ 

197 

198 @staticmethod 

199 def create_generic(field_name, spatial_dimensions, dtype=np.float64, index_dimensions=0, layout='numpy', 

200 index_shape=None, field_type=FieldType.GENERIC) -> 'Field': 

201 """ 

202 Creates a generic field where the field size is not fixed i.e. can be called with arrays of different sizes 

203 

204 Args: 

205 field_name: symbolic name for the field 

206 dtype: numpy data type of the array the kernel is called with later 

207 spatial_dimensions: see documentation of Field 

208 index_dimensions: see documentation of Field 

209 layout: tuple specifying the loop ordering of the spatial dimensions e.g. (2, 1, 0 ) means that 

210 the outer loop loops over dimension 2, the second outer over dimension 1, and the inner loop 

211 over dimension 0. Also allowed: the strings 'numpy' (0,1,..d) or 'reverse_numpy' (d, ..., 1, 0) 

212 index_shape: optional shape of the index dimensions i.e. maximum values allowed for each index dimension, 

213 has to be a list or tuple 

214 field_type: besides the normal GENERIC fields, there are INDEXED fields that store indices of the domain 

215 that should be iterated over, BUFFER fields that are used to generate communication 

216 packing/unpacking kernels, and STAGGERED fields, which store values half-way to the next 

217 cell 

218 """ 

219 if index_shape is not None: 219 ↛ 222line 219 didn't jump to line 222, because the condition on line 219 was never false

220 assert index_dimensions == 0 or index_dimensions == len(index_shape) 

221 index_dimensions = len(index_shape) 

222 if isinstance(layout, str): 222 ↛ 225line 222 didn't jump to line 225, because the condition on line 222 was never false

223 layout = spatial_layout_string_to_tuple(layout, dim=spatial_dimensions) 

224 

225 total_dimensions = spatial_dimensions + index_dimensions 

226 if index_shape is None or len(index_shape) == 0: 226 ↛ 229line 226 didn't jump to line 229, because the condition on line 226 was never false

227 shape = tuple([FieldShapeSymbol([field_name], i) for i in range(total_dimensions)]) 

228 else: 

229 shape = tuple([FieldShapeSymbol([field_name], i) for i in range(spatial_dimensions)] + list(index_shape)) 

230 

231 strides = tuple([FieldStrideSymbol(field_name, i) for i in range(total_dimensions)]) 

232 

233 np_data_type = np.dtype(dtype) 

234 if np_data_type.fields is not None: 234 ↛ 235line 234 didn't jump to line 235, because the condition on line 234 was never true

235 if index_dimensions != 0: 

236 raise ValueError("Structured arrays/fields are not allowed to have an index dimension") 

237 shape += (1,) 

238 strides += (1,) 

239 if field_type == FieldType.STAGGERED and index_dimensions == 0: 239 ↛ 240line 239 didn't jump to line 240, because the condition on line 239 was never true

240 raise ValueError("A staggered field needs at least one index dimension") 

241 

242 return Field(field_name, field_type, dtype, layout, shape, strides) 

243 

244 @staticmethod 

245 def create_from_numpy_array(field_name: str, array: np.ndarray, index_dimensions: int = 0, 

246 field_type=FieldType.GENERIC) -> 'Field': 

247 """Creates a field based on the layout, data type, and shape of a given numpy array. 

248 

249 Kernels created for these kind of fields can only be called with arrays of the same layout, shape and type. 

250 

251 Args: 

252 field_name: symbolic name for the field 

253 array: numpy array 

254 index_dimensions: see documentation of Field 

255 field_type: kind of field 

256 """ 

257 spatial_dimensions = len(array.shape) - index_dimensions 

258 if spatial_dimensions < 1: 258 ↛ 259line 258 didn't jump to line 259, because the condition on line 258 was never true

259 raise ValueError("Too many index dimensions. At least one spatial dimension required") 

260 

261 full_layout = get_layout_of_array(array) 

262 spatial_layout = tuple([i for i in full_layout if i < spatial_dimensions]) 

263 assert len(spatial_layout) == spatial_dimensions 

264 

265 strides = tuple([s // np.dtype(array.dtype).itemsize for s in array.strides]) 

266 shape = tuple(int(s) for s in array.shape) 

267 

268 numpy_dtype = np.dtype(array.dtype) 

269 if numpy_dtype.fields is not None: 269 ↛ 270line 269 didn't jump to line 270, because the condition on line 269 was never true

270 if index_dimensions != 0: 

271 raise ValueError("Structured arrays/fields are not allowed to have an index dimension") 

272 shape += (1,) 

273 strides += (1,) 

274 if field_type == FieldType.STAGGERED and index_dimensions == 0: 274 ↛ 275line 274 didn't jump to line 275, because the condition on line 274 was never true

275 raise ValueError("A staggered field needs at least one index dimension") 

276 

277 return Field(field_name, field_type, array.dtype, spatial_layout, shape, strides) 

278 

279 @staticmethod 

280 def create_fixed_size(field_name: str, shape: Tuple[int, ...], index_dimensions: int = 0, 

281 dtype=np.float64, layout: str = 'numpy', strides: Optional[Sequence[int]] = None, 

282 field_type=FieldType.GENERIC) -> 'Field': 

283 """ 

284 Creates a field with fixed sizes i.e. can be called only with arrays of the same size and layout 

285 

286 Args: 

287 field_name: symbolic name for the field 

288 shape: overall shape of the array 

289 index_dimensions: how many of the trailing dimensions are interpreted as index (as opposed to spatial) 

290 dtype: numpy data type of the array the kernel is called with later 

291 layout: full layout of array, not only spatial dimensions 

292 strides: strides in bytes or None to automatically compute them from shape (assuming no padding) 

293 field_type: kind of field 

294 """ 

295 spatial_dimensions = len(shape) - index_dimensions 

296 assert spatial_dimensions >= 1 

297 

298 if isinstance(layout, str): 

299 layout = layout_string_to_tuple(layout, spatial_dimensions + index_dimensions) 

300 

301 shape = tuple(int(s) for s in shape) 

302 if strides is None: 

303 strides = compute_strides(shape, layout) 

304 else: 

305 assert len(strides) == len(shape) 

306 strides = tuple([s // np.dtype(dtype).itemsize for s in strides]) 

307 

308 numpy_dtype = np.dtype(dtype) 

309 if numpy_dtype.fields is not None: 

310 if index_dimensions != 0: 

311 raise ValueError("Structured arrays/fields are not allowed to have an index dimension") 

312 shape += (1,) 

313 strides += (1,) 

314 if field_type == FieldType.STAGGERED and index_dimensions == 0: 

315 raise ValueError("A staggered field needs at least one index dimension") 

316 

317 spatial_layout = list(layout) 

318 for i in range(spatial_dimensions, len(layout)): 

319 spatial_layout.remove(i) 

320 return Field(field_name, field_type, dtype, tuple(spatial_layout), shape, strides) 

321 

322 def __init__(self, field_name, field_type, dtype, layout, shape, strides): 

323 """Do not use directly. Use static create* methods""" 

324 self._field_name = field_name 

325 assert isinstance(field_type, FieldType) 

326 assert len(shape) == len(strides) 

327 self.field_type = field_type 

328 self._dtype = create_type(dtype) 

329 self._layout = normalize_layout(layout) 

330 self.shape = shape 

331 self.strides = strides 

332 self.latex_name: Optional[str] = None 

333 self.coordinate_origin: tuple[float, sp.Symbol] = sp.Matrix(tuple( 

334 0 for _ in range(self.spatial_dimensions) 

335 )) 

336 self.coordinate_transform = sp.eye(self.spatial_dimensions) 

337 if field_type == FieldType.STAGGERED: 337 ↛ 338line 337 didn't jump to line 338, because the condition on line 337 was never true

338 assert self.staggered_stencil 

339 

340 def new_field_with_different_name(self, new_name): 

341 if self.has_fixed_shape: 

342 return Field(new_name, self.field_type, self._dtype, self._layout, self.shape, self.strides) 

343 else: 

344 return Field.create_generic(new_name, self.spatial_dimensions, self.dtype.numpy_dtype, 

345 self.index_dimensions, self._layout, self.index_shape, self.field_type) 

346 

347 @property 

348 def spatial_dimensions(self) -> int: 

349 return len(self._layout) 

350 

351 @property 

352 def index_dimensions(self) -> int: 

353 return len(self.shape) - len(self._layout) 

354 

355 @property 

356 def ndim(self) -> int: 

357 return len(self.shape) 

358 

359 def values_per_cell(self) -> int: 

360 return functools.reduce(operator.mul, self.index_shape, 1) 

361 

362 @property 

363 def layout(self): 

364 return self._layout 

365 

366 @property 

367 def name(self) -> str: 

368 return self._field_name 

369 

370 @property 

371 def spatial_shape(self) -> Tuple[int, ...]: 

372 return self.shape[:self.spatial_dimensions] 

373 

374 @property 

375 def has_fixed_shape(self): 

376 return is_integer_sequence(self.shape) 

377 

378 @property 

379 def index_shape(self): 

380 return self.shape[self.spatial_dimensions:] 

381 

382 @property 

383 def has_fixed_index_shape(self): 

384 return is_integer_sequence(self.index_shape) 

385 

386 @property 

387 def spatial_strides(self): 

388 return self.strides[:self.spatial_dimensions] 

389 

390 @property 

391 def index_strides(self): 

392 return self.strides[self.spatial_dimensions:] 

393 

394 @property 

395 def dtype(self): 

396 return self._dtype 

397 

398 @property 

399 def itemsize(self): 

400 return self.dtype.numpy_dtype.itemsize 

401 

402 def __repr__(self): 

403 if any(isinstance(s, sp.Symbol) for s in self.spatial_shape): 

404 spatial_shape_str = f'{self.spatial_dimensions}d' 

405 else: 

406 spatial_shape_str = ','.join(str(i) for i in self.spatial_shape) 

407 index_shape_str = ','.join(str(i) for i in self.index_shape) 

408 

409 if self.index_shape: 

410 return f'{self._field_name}({index_shape_str}): {self.dtype}[{spatial_shape_str}]' 

411 else: 

412 return f'{self._field_name}: {self.dtype}[{spatial_shape_str}]' 

413 

414 def __str__(self): 

415 return self.name 

416 

417 def neighbor(self, coord_id, offset): 

418 offset_list = [0] * self.spatial_dimensions 

419 offset_list[coord_id] = offset 

420 return Field.Access(self, tuple(offset_list)) 

421 

422 def neighbors(self, stencil): 

423 return [self.__getitem__(s) for s in stencil] 

424 

425 @property 

426 def center_vector(self): 

427 index_shape = self.index_shape 

428 if len(index_shape) == 0: 

429 return sp.Matrix([self.center]) 

430 elif len(index_shape) == 1: 

431 return sp.Matrix([self(i) for i in range(index_shape[0])]) 

432 elif len(index_shape) == 2: 

433 return sp.Matrix([[self(i, j) for j in range(index_shape[1])] for i in range(index_shape[0])]) 

434 elif len(index_shape) == 3: 

435 return sp.Matrix([[[self(i, j, k) for k in range(index_shape[2])] 

436 for j in range(index_shape[1])] for i in range(index_shape[0])]) 

437 else: 

438 raise NotImplementedError("center_vector is not implemented for more than 3 index dimensions") 

439 

440 @property 

441 def center(self): 

442 center = tuple([0] * self.spatial_dimensions) 

443 return Field.Access(self, center) 

444 

445 def neighbor_vector(self, offset): 

446 """Like neighbor, but returns the entire vector/tensor stored at offset.""" 

447 if self.spatial_dimensions == 2 and len(offset) == 3: 

448 assert offset[2] == 0 

449 offset = offset[:2] 

450 

451 if self.index_dimensions == 0: 

452 return sp.Matrix([self.__getitem__(offset)]) 

453 elif self.index_dimensions == 1: 

454 return sp.Matrix([self.__getitem__(offset)(i) for i in range(self.index_shape[0])]) 

455 elif self.index_dimensions == 2: 

456 return sp.Matrix([[self.__getitem__(offset)(i, k) for k in range(self.index_shape[1])] 

457 for i in range(self.index_shape[0])]) 

458 else: 

459 raise NotImplementedError("neighbor_vector is not implemented for more than 2 index dimensions") 

460 

461 def __getitem__(self, offset): 

462 if type(offset) is np.ndarray: 462 ↛ 463line 462 didn't jump to line 463, because the condition on line 462 was never true

463 offset = tuple(offset) 

464 if type(offset) is str: 464 ↛ 465line 464 didn't jump to line 465, because the condition on line 464 was never true

465 offset = tuple(direction_string_to_offset(offset, self.spatial_dimensions)) 

466 if type(offset) is not tuple: 466 ↛ 467line 466 didn't jump to line 467, because the condition on line 466 was never true

467 offset = (offset,) 

468 if len(offset) != self.spatial_dimensions: 468 ↛ 469line 468 didn't jump to line 469, because the condition on line 468 was never true

469 raise ValueError("Wrong number of spatial indices: " 

470 "Got %d, expected %d" % (len(offset), self.spatial_dimensions)) 

471 return Field.Access(self, offset) 

472 

473 def absolute_access(self, offset, index): 

474 assert FieldType.is_custom(self) 

475 return Field.Access(self, offset, index, is_absolute_access=True) 

476 

477 def interpolated_access(self, 

478 offset: Tuple, 

479 interpolation_mode='linear', 

480 address_mode='BORDER', 

481 allow_textures=True): 

482 """Provides access to field values at non-integer positions 

483 

484 ``interpolated_access`` is similar to :func:`Field.absolute_access` except that 

485 it allows non-integer offsets and automatic handling of out-of-bound accesses. 

486 

487 :param offset: Tuple of spatial coordinates (can be floats) 

488 :param interpolation_mode: One of :class:`pystencils.interpolation_astnodes.InterpolationMode` 

489 :param address_mode: How boundaries are handled can be 'border', 'wrap', 'mirror', 'clamp' 

490 :param allow_textures: Allow implementation by texture accesses on GPUs 

491 """ 

492 from pystencils.interpolation_astnodes import Interpolator 

493 return Interpolator(self, 

494 interpolation_mode, 

495 address_mode, 

496 allow_textures=allow_textures).at(offset) 

497 

498 def staggered_access(self, offset, index=None): 

499 """If this field is a staggered field, it can be accessed using half-integer offsets. 

500 For example, an offset of ``(0, sp.Rational(1,2))`` or ``"E"`` corresponds to the staggered point to the east 

501 of the cell center, i.e. half-way to the eastern-next cell. 

502 If the field stores more than one value per staggered point (e.g. a vector or a tensor), the index (integer or 

503 tuple of integers) refers to which of these values to access. 

504 """ 

505 assert FieldType.is_staggered(self) 

506 

507 offset_orig = offset 

508 if type(offset) is np.ndarray: 508 ↛ 509line 508 didn't jump to line 509, because the condition on line 508 was never true

509 offset = tuple(offset) 

510 if type(offset) is str: 510 ↛ 513line 510 didn't jump to line 513, because the condition on line 510 was never false

511 offset = tuple(direction_string_to_offset(offset, self.spatial_dimensions)) 

512 offset = tuple([o * sp.Rational(1, 2) for o in offset]) 

513 if len(offset) != self.spatial_dimensions: 513 ↛ 514line 513 didn't jump to line 514, because the condition on line 513 was never true

514 raise ValueError("Wrong number of spatial indices: " 

515 "Got %d, expected %d" % (len(offset), self.spatial_dimensions)) 

516 

517 prefactor = 1 

518 neighbor_vec = [0] * len(offset) 

519 for i in range(self.spatial_dimensions): 

520 if (offset[i] + sp.Rational(1, 2)).is_Integer: 

521 neighbor_vec[i] = sp.sign(offset[i]) 

522 neighbor = offset_to_direction_string(neighbor_vec) 

523 if neighbor not in self.staggered_stencil: 523 ↛ 524line 523 didn't jump to line 524, because the condition on line 523 was never true

524 neighbor_vec = inverse_direction(neighbor_vec) 

525 neighbor = offset_to_direction_string(neighbor_vec) 

526 if FieldType.is_staggered_flux(self): 

527 prefactor = -1 

528 if neighbor not in self.staggered_stencil: 528 ↛ 529line 528 didn't jump to line 529, because the condition on line 528 was never true

529 raise ValueError("{} is not a valid neighbor for the {} stencil".format(offset_orig, 

530 self.staggered_stencil_name)) 

531 

532 offset = tuple(sp.Matrix(offset) - sp.Rational(1, 2) * sp.Matrix(neighbor_vec)) 

533 

534 idx = self.staggered_stencil.index(neighbor) 

535 

536 if self.index_dimensions == 1: # this field stores a scalar value at each staggered position 536 ↛ 541line 536 didn't jump to line 541, because the condition on line 536 was never false

537 if index is not None: 537 ↛ 538line 537 didn't jump to line 538, because the condition on line 537 was never true

538 raise ValueError("Cannot specify an index for a scalar staggered field") 

539 return prefactor * Field.Access(self, offset, (idx,)) 

540 else: # this field stores a vector or tensor at each staggered position 

541 if index is None: 

542 raise ValueError("Wrong number of indices: " 

543 "Got %d, expected %d" % (0, self.index_dimensions - 1)) 

544 if type(index) is np.ndarray: 

545 index = tuple(index) 

546 if type(index) is not tuple: 

547 index = (index,) 

548 if self.index_dimensions != len(index) + 1: 

549 raise ValueError("Wrong number of indices: " 

550 "Got %d, expected %d" % (len(index), self.index_dimensions - 1)) 

551 

552 return prefactor * Field.Access(self, offset, (idx, *index)) 

553 

554 def staggered_vector_access(self, offset): 

555 """Like staggered_access, but returns the entire vector/tensor stored at offset.""" 

556 assert FieldType.is_staggered(self) 

557 

558 if self.index_dimensions == 1: 

559 return sp.Matrix([self.staggered_access(offset)]) 

560 elif self.index_dimensions == 2: 

561 return sp.Matrix([self.staggered_access(offset, i) for i in range(self.index_shape[1])]) 

562 elif self.index_dimensions == 3: 

563 return sp.Matrix([[self.staggered_access(offset, (i, k)) for k in range(self.index_shape[2])] 

564 for i in range(self.index_shape[1])]) 

565 else: 

566 raise NotImplementedError("staggered_vector_access is not implemented for more than 3 index dimensions") 

567 

568 @property 

569 def staggered_stencil(self): 

570 assert FieldType.is_staggered(self) 

571 stencils = { 

572 2: { 

573 2: ["W", "S"], # D2Q5 

574 4: ["W", "S", "SW", "NW"] # D2Q9 

575 }, 

576 3: { 

577 3: ["W", "S", "B"], # D3Q7 

578 7: ["W", "S", "B", "BSW", "TSW", "BNW", "TNW"], # D3Q15 

579 9: ["W", "S", "B", "SW", "NW", "BW", "TW", "BS", "TS"], # D3Q19 

580 13: ["W", "S", "B", "SW", "NW", "BW", "TW", "BS", "TS", "BSW", "TSW", "BNW", "TNW"] # D3Q27 

581 } 

582 } 

583 if not self.index_shape[0] in stencils[self.spatial_dimensions]: 583 ↛ 584line 583 didn't jump to line 584, because the condition on line 583 was never true

584 raise ValueError(f"No known stencil has {self.index_shape[0]} staggered points") 

585 return stencils[self.spatial_dimensions][self.index_shape[0]] 

586 

587 @property 

588 def staggered_stencil_name(self): 

589 assert FieldType.is_staggered(self) 

590 return "D%dQ%d" % (self.spatial_dimensions, self.index_shape[0] * 2 + 1) 

591 

592 def __call__(self, *args, **kwargs): 

593 center = tuple([0] * self.spatial_dimensions) 

594 return Field.Access(self, center)(*args, **kwargs) 

595 

596 def hashable_contents(self): 

597 return (self._layout, 

598 self.shape, 

599 self.strides, 

600 self.field_type, 

601 self._field_name, 

602 self.latex_name, 

603 self._dtype) 

604 

605 def __hash__(self): 

606 return hash(self.hashable_contents()) 

607 

608 def __eq__(self, other): 

609 if not isinstance(other, Field): 609 ↛ 610line 609 didn't jump to line 610, because the condition on line 609 was never true

610 return False 

611 return self.hashable_contents() == other.hashable_contents() 

612 

613 @property 

614 def physical_coordinates(self): 

615 if hasattr(self.coordinate_transform, '__call__'): 

616 return self.coordinate_transform(self.coordinate_origin + pystencils.x_vector(self.spatial_dimensions)) 

617 else: 

618 return self.coordinate_transform @ (self.coordinate_origin + pystencils.x_vector(self.spatial_dimensions)) 

619 

620 @property 

621 def physical_coordinates_staggered(self): 

622 return self.coordinate_transform @ \ 

623 (self.coordinate_origin + pystencils.x_staggered_vector(self.spatial_dimensions)) 

624 

625 def index_to_physical(self, index_coordinates: sp.Matrix, staggered=False): 

626 if staggered: 

627 index_coordinates = sp.Matrix([0.5] * len(self.coordinate_origin)) + index_coordinates 

628 if hasattr(self.coordinate_transform, '__call__'): 

629 return self.coordinate_transform(self.coordinate_origin + index_coordinates) 

630 else: 

631 return self.coordinate_transform @ (self.coordinate_origin + index_coordinates) 

632 

633 def physical_to_index(self, physical_coordinates: sp.Matrix, staggered=False): 

634 if hasattr(self.coordinate_transform, '__call__'): 

635 if hasattr(self.coordinate_transform, 'inv'): 

636 return self.coordinate_transform.inv()(physical_coordinates) - self.coordinate_origin 

637 else: 

638 idx = sp.Matrix(sp.symbols(f'index_coordinates:{self.ndim}', real=True)) 

639 rtn = sp.solve(self.index_to_physical(idx) - physical_coordinates, idx) 

640 assert rtn, f'Could not find inverese of coordinate_transform: {self.index_to_physical(idx)}' 

641 return rtn 

642 

643 else: 

644 rtn = self.coordinate_transform.inv() @ physical_coordinates - self.coordinate_origin 

645 if staggered: 

646 rtn = sp.Matrix([i - 0.5 for i in rtn]) 

647 

648 return rtn 

649 

650 def set_coordinate_origin_to_field_center(self): 

651 self.coordinate_origin = -sp.Matrix([i / 2 for i in self.spatial_shape]) 

652 

653 # noinspection PyAttributeOutsideInit,PyUnresolvedReferences 

654 class Access(TypedSymbol, AbstractField.AbstractAccess): 

655 """Class representing a relative access into a `Field`. 

656 

657 This class behaves like a normal sympy Symbol, it is actually derived from it. One can built up 

658 sympy expressions using field accesses, solve for them, etc. 

659 

660 Examples: 

661 >>> vector_field_2d = fields("v(2): double[2D]") # create a 2D vector field 

662 >>> northern_neighbor_y_component = vector_field_2d[0, 1](1) 

663 >>> northern_neighbor_y_component 

664 v_N^1 

665 >>> central_y_component = vector_field_2d(1) 

666 >>> central_y_component 

667 v_C^1 

668 >>> central_y_component.get_shifted(1, 0) # move the existing access 

669 v_E^1 

670 >>> central_y_component.at_index(0) # change component 

671 v_C^0 

672 """ 

673 _iterable = False # see https://i10git.cs.fau.de/pycodegen/pystencils/-/merge_requests/166#note_10680 

674 

675 def __new__(cls, name, *args, **kwargs): 

676 obj = Field.Access.__xnew_cached_(cls, name, *args, **kwargs) 

677 return obj 

678 

679 def __new_stage2__(self, field, offsets=(0, 0, 0), idx=None, is_absolute_access=False, dtype=None): 

680 field_name = field.name 

681 offsets_and_index = (*offsets, *idx) if idx is not None else offsets 

682 constant_offsets = not any([isinstance(o, sp.Basic) and not o.is_Integer for o in offsets_and_index]) 

683 

684 if not idx: 

685 idx = tuple([0] * field.index_dimensions) 

686 

687 if constant_offsets: 687 ↛ 701line 687 didn't jump to line 701, because the condition on line 687 was never false

688 offset_name = offset_to_direction_string(offsets) 

689 if field.index_dimensions == 0: 

690 superscript = None 

691 elif field.index_dimensions == 1: 691 ↛ 694line 691 didn't jump to line 694, because the condition on line 691 was never false

692 superscript = str(idx[0]) 

693 else: 

694 idx_str = ",".join([str(e) for e in idx]) 

695 superscript = idx_str 

696 if field.has_fixed_index_shape and not isinstance(field.dtype, StructType): 696 ↛ 704line 696 didn't jump to line 704, because the condition on line 696 was never false

697 for i, bound in zip(idx, field.index_shape): 

698 if i >= bound: 698 ↛ 699line 698 didn't jump to line 699, because the condition on line 698 was never true

699 raise ValueError("Field index out of bounds") 

700 else: 

701 offset_name = hashlib.md5(pickle.dumps(offsets_and_index)).hexdigest()[:12] 

702 superscript = None 

703 

704 symbol_name = f"{field_name}_{offset_name}" 

705 if superscript is not None: 

706 symbol_name += "^" + superscript 

707 

708 obj = super(Field.Access, self).__xnew__(self, symbol_name, field.dtype) 

709 obj._field = field 

710 obj._offsets = [] 

711 for o in offsets: 

712 if isinstance(o, sp.Basic): 

713 obj._offsets.append(o) 

714 else: 

715 obj._offsets.append(int(o)) 

716 obj._offsets = tuple(sp.sympify(obj._offsets)) 

717 obj._offsetName = offset_name 

718 obj._superscript = superscript 

719 obj._index = idx 

720 

721 obj._indirect_addressing_fields = set() 

722 for e in chain(obj._offsets, obj._index): 

723 if isinstance(e, sp.Basic): 

724 obj._indirect_addressing_fields.update(a.field for a in e.atoms(Field.Access)) 

725 

726 obj._is_absolute_access = is_absolute_access 

727 return obj 

728 

729 def __getnewargs__(self): 

730 return self.field, self.offsets, self.index, self.is_absolute_access, self.dtype 

731 

732 def __getnewargs_ex__(self): 

733 return (self.field, self.offsets, self.index, self.is_absolute_access, self.dtype), {} 

734 

735 # noinspection SpellCheckingInspection 

736 __xnew__ = staticmethod(__new_stage2__) 

737 # noinspection SpellCheckingInspection 

738 __xnew_cached_ = staticmethod(cacheit(__new_stage2__)) 

739 

740 def __call__(self, *idx): 

741 if self._index != tuple([0] * self.field.index_dimensions): 741 ↛ 742line 741 didn't jump to line 742, because the condition on line 741 was never true

742 raise ValueError("Indexing an already indexed Field.Access") 

743 

744 idx = tuple(idx) 

745 

746 if self.field.index_dimensions == 0 and idx == (0,): 746 ↛ 747line 746 didn't jump to line 747, because the condition on line 746 was never true

747 idx = () 

748 

749 if len(idx) != self.field.index_dimensions: 749 ↛ 750line 749 didn't jump to line 750, because the condition on line 749 was never true

750 raise ValueError("Wrong number of indices: " 

751 "Got %d, expected %d" % (len(idx), self.field.index_dimensions)) 

752 return Field.Access(self.field, self._offsets, idx, dtype=self.dtype) 

753 

754 def __getitem__(self, *idx): 

755 return self.__call__(*idx) 

756 

757 @property 

758 def field(self) -> 'Field': 

759 """Field that the Access points to""" 

760 return self._field 

761 

762 @property 

763 def offsets(self) -> Tuple: 

764 """Spatial offset as tuple""" 

765 return self._offsets 

766 

767 @property 

768 def required_ghost_layers(self) -> int: 

769 """Largest spatial distance that is accessed.""" 

770 return int(np.max(np.abs(self._offsets))) 

771 

772 @property 

773 def nr_of_coordinates(self): 

774 return len(self._offsets) 

775 

776 @property 

777 def offset_name(self) -> str: 

778 """Spatial offset as string, East-West for x, North-South for y and Top-Bottom for z coordinate. 

779 

780 Example: 

781 >>> f = fields("f: double[2D]") 

782 >>> f[1, 1].offset_name # north-east 

783 'NE' 

784 """ 

785 return self._offsetName 

786 

787 @property 

788 def index(self): 

789 """Value of index coordinates as tuple.""" 

790 return self._index 

791 

792 def neighbor(self, coord_id: int, offset: int) -> 'Field.Access': 

793 """Returns a new Access with changed spatial coordinates. 

794 

795 Args: 

796 coord_id: index of the coordinate to change (0 for x, 1 for y,...) 

797 offset: incremental change of this coordinate 

798 

799 Example: 

800 >>> f = fields('f: [2D]') 

801 >>> f[0,0].neighbor(coord_id=1, offset=-1) 

802 f_S 

803 """ 

804 offset_list = list(self.offsets) 

805 offset_list[coord_id] += offset 

806 return Field.Access(self.field, tuple(offset_list), self.index, dtype=self.dtype) 

807 

808 def get_shifted(self, *shift) -> 'Field.Access': 

809 """Returns a new Access with changed spatial coordinates 

810 

811 Example: 

812 >>> f = fields("f: [2D]") 

813 >>> f[0,0].get_shifted(1, 1) 

814 f_NE 

815 """ 

816 return Field.Access(self.field, 

817 tuple(a + b for a, b in zip(shift, self.offsets)), 

818 self.index, 

819 dtype=self.dtype) 

820 

821 def at_index(self, *idx_tuple) -> 'Field.Access': 

822 """Returns new Access with changed index. 

823 

824 Example: 

825 >>> f = fields("f(9): [2D]") 

826 >>> f(0).at_index(8) 

827 f_C^8 

828 """ 

829 return Field.Access(self.field, self.offsets, idx_tuple, dtype=self.dtype) 

830 

831 def _eval_subs(self, old, new): 

832 return Field.Access(self.field, 

833 tuple(sp.sympify(a).subs(old, new) for a in self.offsets), 

834 tuple(sp.sympify(a).subs(old, new) for a in self.index), 

835 dtype=self.dtype) 

836 

837 @property 

838 def is_absolute_access(self) -> bool: 

839 """Indicates if a field access is relative to the loop counters (this is the default) or absolute""" 

840 return self._is_absolute_access 

841 

842 @property 

843 def indirect_addressing_fields(self) -> Set['Field']: 

844 """Returns a set of fields that the access depends on. 

845 

846 e.g. f[index_field[1, 0]], the outer access to f depends on index_field 

847 """ 

848 return self._indirect_addressing_fields 

849 

850 def _hashable_content(self): 

851 super_class_contents = super(Field.Access, self)._hashable_content() 

852 return (super_class_contents, self._field.hashable_contents(), *self._index, *self._offsets) 

853 

854 def _staggered_offset(self, offsets, index): 

855 assert FieldType.is_staggered(self._field) 

856 neighbor = self._field.staggered_stencil[index] 

857 neighbor = direction_string_to_offset(neighbor, self._field.spatial_dimensions) 

858 return [(o + sp.Rational(int(neighbor[i]), 2)) for i, o in enumerate(offsets)] 

859 

860 def _latex(self, _): 

861 n = self._field.latex_name if self._field.latex_name else self._field.name 

862 offset_str = ",".join([sp.latex(o) for o in self.offsets]) 

863 if FieldType.is_staggered(self._field): 

864 offset_str = ",".join([sp.latex(self._staggered_offset(self.offsets, self.index[0])[i]) 

865 for i in range(len(self.offsets))]) 

866 if self.is_absolute_access: 

867 offset_str = f"\\mathbf{offset_str}" 

868 elif self.field.spatial_dimensions > 1: 

869 offset_str = f"({offset_str})" 

870 

871 if FieldType.is_staggered(self._field): 

872 if self.index and self.field.index_dimensions > 1: 

873 return "{{%s}_{%s}^{%s}}" % (n, offset_str, self.index[1:] 

874 if len(self.index) > 2 else self.index[1]) 

875 else: 

876 return "{{%s}_{%s}}" % (n, offset_str) 

877 else: 

878 if self.index and self.field.index_dimensions > 0: 

879 return "{{%s}_{%s}^{%s}}" % (n, offset_str, self.index if len(self.index) > 1 else self.index[0]) 

880 else: 

881 return "{{%s}_{%s}}" % (n, offset_str) 

882 

883 def __str__(self): 

884 n = self._field.latex_name if self._field.latex_name else self._field.name 

885 offset_str = ",".join([sp.latex(o) for o in self.offsets]) 

886 if FieldType.is_staggered(self._field): 

887 offset_str = ",".join([sp.latex(self._staggered_offset(self.offsets, self.index[0])[i]) 

888 for i in range(len(self.offsets))]) 

889 if self.is_absolute_access: 

890 offset_str = f"[abs]{offset_str}" 

891 

892 if FieldType.is_staggered(self._field): 

893 if self.index and self.field.index_dimensions > 1: 

894 return f"{n}[{offset_str}]({self.index[1:] if len(self.index) > 2 else self.index[1]})" 

895 else: 

896 return f"{n}[{offset_str}]" 

897 else: 

898 if self.index and self.field.index_dimensions > 0: 

899 return f"{n}[{offset_str}]({self.index if len(self.index) > 1 else self.index[0]})" 

900 else: 

901 return f"{n}[{offset_str}]" 

902 

903 

904def get_layout_from_strides(strides: Sequence[int], index_dimension_ids: Optional[List[int]] = None): 

905 index_dimension_ids = [] if index_dimension_ids is None else index_dimension_ids 

906 coordinates = list(range(len(strides))) 

907 relevant_strides = [stride for i, stride in enumerate(strides) if i not in index_dimension_ids] 

908 result = [x for (y, x) in sorted(zip(relevant_strides, coordinates), key=lambda pair: pair[0], reverse=True)] 

909 return normalize_layout(result) 

910 

911 

912def get_layout_of_array(arr: np.ndarray, index_dimension_ids: Optional[List[int]] = None): 

913 """ Returns a list indicating the memory layout (linearization order) of the numpy array. 

914 

915 Examples: 

916 >>> get_layout_of_array(np.zeros([3,3,3])) 

917 (0, 1, 2) 

918 

919 In this example the loop over the zeroth coordinate should be the outermost loop, 

920 followed by the first and second. Elements arr[x,y,0] and arr[x,y,1] are adjacent in memory. 

921 Normally constructed numpy arrays have this order, however by stride tricks or other frameworks, arrays 

922 with different memory layout can be created. 

923 

924 The index_dimension_ids parameter leaves specifies which coordinates should not be 

925 """ 

926 index_dimension_ids = [] if index_dimension_ids is None else index_dimension_ids 

927 return get_layout_from_strides(arr.strides, index_dimension_ids) 

928 

929 

930def create_numpy_array_with_layout(shape, layout, alignment=False, byte_offset=0, **kwargs): 

931 """Creates numpy array with given memory layout. 

932 

933 Args: 

934 shape: shape of the resulting array 

935 layout: layout as tuple, where the coordinates are ordered from slow to fast 

936 alignment: number of bytes to align the beginning and the innermost coordinate to, or False for no alignment 

937 byte_offset: only used when alignment is specified, align not beginning but address at this offset 

938 mostly used to align first inner cell, not ghost cells 

939 

940 Example: 

941 >>> res = create_numpy_array_with_layout(shape=(2, 3, 4, 5), layout=(3, 2, 0, 1)) 

942 >>> res.shape 

943 (2, 3, 4, 5) 

944 >>> get_layout_of_array(res) 

945 (3, 2, 0, 1) 

946 """ 

947 assert set(layout) == set(range(len(shape))), "Wrong layout descriptor" 

948 cur_layout = list(range(len(shape))) 

949 swaps = [] 

950 for i in range(len(layout)): 

951 if cur_layout[i] != layout[i]: 

952 index_to_swap_with = cur_layout.index(layout[i]) 

953 swaps.append((i, index_to_swap_with)) 

954 cur_layout[i], cur_layout[index_to_swap_with] = cur_layout[index_to_swap_with], cur_layout[i] 

955 assert tuple(cur_layout) == tuple(layout) 

956 

957 shape = list(shape) 

958 for a, b in swaps: 

959 shape[a], shape[b] = shape[b], shape[a] 

960 

961 if not alignment: 

962 res = np.empty(shape, order='c', **kwargs) 

963 else: 

964 res = aligned_empty(shape, alignment, byte_offset=byte_offset, **kwargs) 

965 

966 for a, b in reversed(swaps): 

967 res = res.swapaxes(a, b) 

968 return res 

969 

970 

971def spatial_layout_string_to_tuple(layout_str: str, dim: int) -> Tuple[int, ...]: 

972 if layout_str in ('fzyx', 'zyxf'): 972 ↛ 973line 972 didn't jump to line 973, because the condition on line 972 was never true

973 assert dim <= 3 

974 return tuple(reversed(range(dim))) 

975 

976 if layout_str in ('fzyx', 'f', 'reverse_numpy', 'SoA'): 

977 return tuple(reversed(range(dim))) 

978 elif layout_str in ('c', 'numpy', 'AoS'): 978 ↛ 980line 978 didn't jump to line 980, because the condition on line 978 was never false

979 return tuple(range(dim)) 

980 raise ValueError("Unknown layout descriptor " + layout_str) 

981 

982 

983def layout_string_to_tuple(layout_str, dim): 

984 layout_str = layout_str.lower() 

985 if layout_str == 'fzyx' or layout_str == 'soa': 985 ↛ 988line 985 didn't jump to line 988, because the condition on line 985 was never false

986 assert dim <= 4 

987 return tuple(reversed(range(dim))) 

988 elif layout_str == 'zyxf' or layout_str == 'aos': 

989 assert dim <= 4 

990 return tuple(reversed(range(dim - 1))) + (dim - 1,) 

991 elif layout_str == 'f' or layout_str == 'reverse_numpy': 

992 return tuple(reversed(range(dim))) 

993 elif layout_str == 'c' or layout_str == 'numpy': 

994 return tuple(range(dim)) 

995 raise ValueError("Unknown layout descriptor " + layout_str) 

996 

997 

998def normalize_layout(layout): 

999 """Takes a layout tuple and subtracts the minimum from all entries""" 

1000 min_entry = min(layout) 

1001 return tuple(i - min_entry for i in layout) 

1002 

1003 

1004def compute_strides(shape, layout): 

1005 """ 

1006 Computes strides assuming no padding exists 

1007 

1008 Args: 

1009 shape: shape (size) of array 

1010 layout: layout specification as tuple 

1011 

1012 Returns: 

1013 strides in elements, not in bytes 

1014 """ 

1015 dim = len(shape) 

1016 assert len(layout) == dim 

1017 assert len(set(layout)) == dim 

1018 strides = [0] * dim 

1019 product = 1 

1020 for j in reversed(layout): 

1021 strides[j] = product 

1022 product *= shape[j] 

1023 return tuple(strides) 

1024 

1025 

1026# ---------------------------------------- Parsing of string in fields() function -------------------------------------- 

1027 

1028field_description_regex = re.compile(r""" 

1029 \s* # ignore leading white spaces 

1030 (\w+) # identifier is a sequence of alphanumeric characters, is stored in first group 

1031 (?: # optional index specification e.g. (1, 4, 2) 

1032 \s* 

1033 \( 

1034 ([^\)]+) # read everything up to closing bracket 

1035 \) 

1036 \s* 

1037 )? 

1038 \s*,?\s* # ignore trailing white spaces and comma 

1039""", re.VERBOSE) 

1040 

1041type_description_regex = re.compile(r""" 

1042 \s* 

1043 (\w+)? # optional dtype 

1044 \s* 

1045 \[ 

1046 ([^\]]+) 

1047 \] 

1048 \s* 

1049""", re.VERBOSE | re.IGNORECASE) 

1050 

1051 

1052def _parse_part1(d): 

1053 result = field_description_regex.match(d) 

1054 while result: 

1055 name, index_str = result.group(1), result.group(2) 

1056 index = tuple(int(e) for e in index_str.split(",")) if index_str else () 

1057 yield name, index 

1058 d = d[result.end():] 

1059 result = field_description_regex.match(d) 

1060 

1061 

1062def _parse_description(description): 

1063 

1064 def parse_part2(d): 

1065 result = type_description_regex.match(d) 

1066 if result: 1066 ↛ 1080line 1066 didn't jump to line 1080, because the condition on line 1066 was never false

1067 data_type_str, size_info = result.group(1), result.group(2).strip().lower() 

1068 if data_type_str is None: 

1069 data_type_str = 'float64' 

1070 data_type_str = data_type_str.lower().strip() 

1071 

1072 if not data_type_str: 1072 ↛ 1073line 1072 didn't jump to line 1073, because the condition on line 1072 was never true

1073 data_type_str = 'float64' 

1074 if size_info.endswith('d'): 1074 ↛ 1077line 1074 didn't jump to line 1077, because the condition on line 1074 was never false

1075 size_info = int(size_info[:-1]) 

1076 else: 

1077 size_info = tuple(int(e) for e in size_info.split(",")) 

1078 return data_type_str, size_info 

1079 else: 

1080 raise ValueError("Could not parse field description") 

1081 

1082 if ':' in description: 1082 ↛ 1085line 1082 didn't jump to line 1085, because the condition on line 1082 was never false

1083 field_description, field_info = description.split(':') 

1084 else: 

1085 field_description, field_info = description, 'float64[2D]' 

1086 

1087 fields_info = [e for e in _parse_part1(field_description)] 

1088 if not field_info: 1088 ↛ 1089line 1088 didn't jump to line 1089, because the condition on line 1088 was never true

1089 raise ValueError("Could not parse field description") 

1090 

1091 data_type, size = parse_part2(field_info) 

1092 return fields_info, data_type, size