collapse empty dimensions of time-domain grid slices by default and remove snap argument from get_array_metadata #1456

oskooi · 2020-12-20T03:45:57Z

Fixes #759.

Fixes #1431 based on the source of the bug described in #1431 (comment).

The snap parameter of get_array_metadata has been removed since it does not really serve a practical purpose compared with the argument collapse which does interpolation. Also, documentation has been added for the collapse argument which was missing.

The tests described in #1431 (comment) now produce the correct results:

import meep as mp
import numpy as np

sim = mp.Simulation(cell_size=mp.Vector3(5,5),
                    resolution=10)

sim.init_sim()

ds = np.linspace(0,0.1,11)

print("sweep grid slice along y")
for d in ds:
    (x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(0,d), size=mp.Vector3(0.5,0), collapse=False)
    if (w.ndim == 2):
	print('{:.3f}, {:.3f}, {:.3f}'.format(d,np.sum(w[:,0]),np.sum(w[:,1])))
    elif (w.ndim == 1):
	print('{:.3f}, {:.3f}'.format(d,np.sum(w)))

print("sweep grid slice along x")
for d in ds:
    (x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(d,0), size=mp.Vector3(0,0.5), collapse=False)
    if (w.ndim == 2):
	print('{:.3f}, {:.3f}, {:.3f}'.format(d,np.sum(w[0,:]),np.sum(w[1,:])))
    elif (w.ndim == 1):
	print('{:.3f}, {:.3f}'.format(d,np.sum(w)))

sweep grid slice along y
0.000, 0.250, 0.250
0.010, 0.200, 0.300
0.020, 0.150, 0.350
0.030, 0.100, 0.400
0.040, 0.050, 0.450
0.050, 0.500
0.060, 0.450, 0.050
0.070, 0.400, 0.100
0.080, 0.350, 0.150
0.090, 0.300, 0.200
0.100, 0.250, 0.250
sweep grid slice along x
0.000, 0.250, 0.250
0.010, 0.200, 0.300
0.020, 0.150, 0.350
0.030, 0.100, 0.400
0.040, 0.050, 0.450
0.050, 0.500
0.060, 0.450, 0.050
0.070, 0.400, 0.100
0.080, 0.350, 0.150
0.090, 0.300, 0.200
0.100, 0.250, 0.250

Note how the sum of the weights from each of the uncollapsed dimensions (second and third columns) changes as the grid slice is swept across a distance of one pixel (first column).

import meep as mp
import numpy as np

sim = mp.Simulation(cell_size=mp.Vector3(5,5,5),
                    resolution=10)

sim.init_sim()

ds = np.linspace(0,0.1,11)+0.058712

print("grid slice along z")
for d in ds:
    (x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(0,0,d), size=mp.Vector3(0.5,0.5,0), collapse=False)
    if (w.ndim == 3):
        print('{:.3f}, {:.3f}, {:.3f}'.format(d,np.sum(w[:,:,0]),np.sum(w[:,:,1])))
    elif (w.ndim == 2):
        print('{:.3f}, {:.3f}'.format(d,np.sum(w)))

print("grid slice along y")
for d in ds:
    (x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(0,d,0), size=mp.Vector3(0.5,0,0.5), collapse=False)
    if (w.ndim == 3):
        print('{:.3f}, {:.3f}, {:.3f}'.format(d,np.sum(w[:,0,:]),np.sum(w[:,1,:])))
    elif (w.ndim == 2):
        print('{:.3f}, {:.3f}'.format(d,np.sum(w)))

print("grid slice along x")
for d in ds:
    (x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(d,0,0), size=mp.Vector3(0,0.5,0.5), collapse=False)
    if (w.ndim == 3):
	print('{:.3f}, {:.3f}, {:.3f}'.format(d,np.sum(w[0,:,:]),np.sum(w[1,:,:])))
    elif (w.ndim == 2):
	print('{:.3f}, {:.3f}'.format(d,np.sum(w)))

grid slice along z
0.059, 0.228, 0.022
0.069, 0.203, 0.047
0.079, 0.178, 0.072
0.089, 0.153, 0.097
0.099, 0.128, 0.122
0.109, 0.103, 0.147
0.119, 0.078, 0.172
0.129, 0.053, 0.197
0.139, 0.028, 0.222
0.149, 0.003, 0.247
0.159, 0.228, 0.022
grid slice along y
0.059, 0.228, 0.022
0.069, 0.203, 0.047
0.079, 0.178, 0.072
0.089, 0.153, 0.097
0.099, 0.128, 0.122
0.109, 0.103, 0.147
0.119, 0.078, 0.172
0.129, 0.053, 0.197
0.139, 0.028, 0.222
0.149, 0.003, 0.247
0.159, 0.228, 0.022
grid slice along x
0.059, 0.228, 0.022
0.069, 0.203, 0.047
0.079, 0.178, 0.072
0.089, 0.153, 0.097
0.099, 0.128, 0.122
0.109, 0.103, 0.147
0.119, 0.078, 0.172
0.129, 0.053, 0.197
0.139, 0.028, 0.222
0.149, 0.003, 0.247
0.159, 0.228, 0.022

Note: the documentation in Python_User_Interface.md has been regenerated using meep/doc/generate_py_api.py which is why an unrelated change due to plot2D (from #1446) appears.

…of grid slice

python/simulation.py

ianwilliamson · 2020-12-20T05:31:49Z

python/simulation.py

        """
        This routine provides geometric information useful for interpreting the arrays
        returned by `get_array` or `get_dft_array` for the spatial region defined by `vol`
        or `center`/`size`. In both cases, the return value is a tuple `(x,y,z,w)`, where:

-        + `x,y,z` are 1d NumPy arrays storing the $x,y,z$ coordinates of the points in the
+        + `x,y,z` are tuples storing the $x,y,z$ coordinates of the points in the


Just another drive by comment following this from the reference to #1431... x, y, and z look like they're NumPy arrays, not tuples.

Thanks for calling this out. The "tics" return values (x,y,z) should be NumPy arrays but in fact they are set up as tuples:

meep/python/simulation.py

Lines 3278 to 3289 in f5fc4cd

xyzw_vector = self.fields.get_array_metadata(v, collapse, snap)

offset, tics = 0, []

for n in range(3):

N = int(xyzw_vector[offset])

tics.append( xyzw_vector[offset+1:offset+1+N] )

offset += 1+N

wshape=[len(t) for t in tics if len(t)>1]

weights=np.reshape(xyzw_vector[offset:],wshape)

if return_pw:

points=[ mp.Vector3(x,y,z) for x in tics[0] for y in tics[1] for z in tics[2] ]

return points,weights

return tics + [weights]

I verified this using type(x), type(y), type(z) and the results was a tuple. (w though is a NumPy array.) If we really want these return values as 1d NumPy arrays then we will need to modify get_array_metadata in a separate PR. I just updated the documentation here to make it consistent with what is actually happening.

Ah, good to know. If you wanted, you could also call np.asarray() on them to do the conversion in Python.

This looks like a SWIG problem — self.fields.get_array_metadata(v, collapse, snap) is nowadays returning a tuple, whereas the intention (and I think the original behavior) was to return a NumPy array.

We definitely want a numpy array here — a tuple is a terrible data structure for a large array. And we want to return a numpy array directly, not return a tuple and then convert to numpy.

fields::get_array_metadata returns a std::vector<double>, which we want SWIG to convert to a NumPy array. I'm not sure why that isn't happening.

According to the documentation for SWIG's std_vector.i (which we are using), it's supposed to be using a "list object in the target language" for output typemaps.

I was able to fix it by adding:

// use NumPy arrays for returning common std::vector types: %typemap(out) std::vector<double> { npy_intp vec_len = (npy_intp) $1.size(); $result = PyArray_SimpleNew(1, &vec_len, NPY_DOUBLE); memcpy(PyArray_DATA((PyArrayObject*) $result), &$1[0], vec_len * sizeof(double)); }

to meep.i (after %template(DoubleVector) std::vector<double>;).

We can do that in a separate PR.

stevengj · 2020-12-22T18:32:02Z

Note that until we get proper interpolation of the time-domain slices (#759), we may want to have a snap option here if you want the metadata for a time-domain slice.

stevengj · 2020-12-22T18:41:22Z

Seems like it should be easy to fix #759 nowadays, so that we can remove the snap option completely? For the time-domain get_slice, just set snap=false and call collapse_array on the result before returning it.

ianwilliamson · 2020-12-22T22:12:51Z

Incorporating these changes does seem to eliminate the ASAN finding that I reported in #1431. ASAN now only returns the result reported in #1432, which is an issue in a separate function.

oskooi · 2020-12-27T04:26:46Z

The snap option has been removed and the time-domain slices are now all collapsed. The docs has also been updated with these changes.

All the tests in the make check suite are passing on my local machine. However, it seems that the Travis CI tests are disabled?

oskooi · 2020-12-27T04:31:42Z

The latest changes also affect the call to get_array_slice_dimensions in material_grids_add_gradient of meepgeom.cpp. Thus, this PR may also fix #1432?

…hdf5 only snaps to grid

…avity_arrayslice.py test

python/simulation.py

stevengj · 2020-12-29T18:58:00Z

Would be good to re-instate the snap option for get_array_slice (but not get_array_metadata), if only for comparison the the output_h5 feature (which snaps) in order to re-instate the tests (and also to ease migration if someone is migrating from an HDF5-based workflow to a Python get_array_slice workflow.

In the long run, I'm not sure that the current output_h5 implementation will support collapsing interpolated slices, because output_h5 is designed for saving huge datasets that may not fit into the memory of a single process — since it never loads the uncollapsed array into a single process's memory, implementing collapse_array could be difficult. However, we could still add a snap=false option to output_h5 — in that case, we could just call get_array_slice and then h5file::write, circumventing the whole parallel-I/O implementation entirely.

oskooi · 2020-12-29T19:02:52Z

Verifying types of return values of get_array_metadata:

import meep as mp

sim = mp.Simulation(cell_size=mp.Vector3(5,5,5),
                    resolution=20)

sim.init_sim()

(x,y,z,w) = sim.get_array_metadata(center=mp.Vector3(0,0,0.149), size=mp.Vector3(4.5,2.4,0))
print(len(x),len(y),len(z),w.shape)

print(type(x),type(y),type(z),type(w))

output

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
time for choose_chunkdivision = 2.7325e-05 s
Working in 3D dimensions.
Computational cell is 5 x 5 x 5 with resolution 20
time for set_epsilon = 3.19013 s
-----------
92 50 1 (92, 50)
<class 'tuple'> <class 'tuple'> <class 'tuple'> <class 'numpy.ndarray'>

stevengj · 2020-12-29T19:53:43Z

Note that if snap==true you should omit the IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0) weights in get_array_slice.

oskooi · 2020-12-30T05:26:14Z

The snap argument to get_array has been re-instated as well as the tests which involve comparing the output from get_array with output_hdf5 (tests/array-slice-ll.cpp, python/tests/simulation.py). The make check suite is passing on my local machine.

oskooi · 2020-12-30T06:05:28Z

All tests from the most recent commit (33dfce9) are passing on Travis (https://travis-ci.com/github/oskooi/meep/builds/211270487).

python/simulation.py

stevengj · 2020-12-31T18:30:51Z

Note that the weights returned by get_array_metadata are somewhat inconsistent. Realize that these weights include three things:

the volume element dV (≠ constant in cylindrical coords)
an interpolation weight for slices
a fractional weight at the edges of the region, corresponding to the fractional dV that those pixels contribute to an integral.

Currently, the get_array_metadata weights include 1 and 3 — weight (2) goes away when you collapse empty dimensions (because interpolation weights sum to 1). With this PR, the get_array_slice fields include weights (3). So, if you use sum(fields*weights) to compute an integral, you are including weights (3) twice, making a first-order error.

Note that the dft_fields already include weights 3 as well, so there is the same first-order error in the current release for computing integrals of DFT fields with these integration weights. And the current release makes a first-order error in get_array_slice when it snaps the field — so in some sense you are swapping one first-order error for another.

We can probably fix this in another PR. The fix is simple: this line should be changed to fields[i] = dV0 + dV1 * loop_i2; (computing weights 1 only). But some of the test results might change slightly.

Correction the dft_fields do not include weight 3, so there is no problem there. The only problem is with array slices?

python/simulation.py

src/array_slice.cpp

…emove snap argument from get_array_metadata (NanoComp#1456) * fix for get_array_metadata related to collapsing of empty dimensions of grid slice * remove snap_empty_dimensions and collapse arrays in get_array_ * fix get_source_slice and get_array_metadata * update hard-coded validation binary files for cavity_arrayslice.py * remove collapse_empty_dimensions argument from tests/array_metadata.cpp * add array_to_all for parallel simulations and update docs * remove tests/array-slice-ll.cpp from make check suite because output_hdf5 only snaps to grid * add interpolation weights to slices and update hard-coded values of cavity_arrayslice.py test * re-instate snap argument for do_get_array_slice * fix bug in do_get_array_slice related to snapping * re-instate snap_empty_dimensions argument in get_array_slice_dimensions * include snap_empty_dimensions in meep.i * fixes * condense/simplify Co-authored-by: Steven G. Johnson <stevenj@alum.mit.edu>

fix for get_array_metadata related to collapsing of empty dimensions …

9c663f7

…of grid slice

oskooi added the enhancement label Dec 20, 2020

ianwilliamson reviewed Dec 20, 2020

View reviewed changes

python/simulation.py Show resolved Hide resolved

ianwilliamson reviewed Dec 20, 2020

View reviewed changes

oskooi added 5 commits December 26, 2020 15:22

remove snap_empty_dimensions and collapse arrays in get_array_

5cd145c

fix get_source_slice and get_array_metadata

823d89f

update hard-coded validation binary files for cavity_arrayslice.py

cc20805

remove collapse_empty_dimensions argument from tests/array_metadata.cpp

f468cfb

add array_to_all for parallel simulations and update docs

c4347c3

oskooi changed the title ~~fix for get_array_metadata related to collapsing of empty dimensions of grid slice~~ collapse empty dimensions of grid slices by default and remove snap argument Dec 27, 2020

oskooi added 2 commits December 29, 2020 09:05

remove tests/array-slice-ll.cpp from make check suite because output_…

cd043e1

…hdf5 only snaps to grid

add interpolation weights to slices and update hard-coded values of c…

0614b0f

…avity_arrayslice.py test

stevengj reviewed Dec 29, 2020

View reviewed changes

python/simulation.py Outdated Show resolved Hide resolved

stevengj mentioned this pull request Dec 29, 2020

return numpy arrays from std::vectors #1458

Merged

oskooi added 2 commits December 29, 2020 20:38

re-instate snap argument for do_get_array_slice

e6f7123

fix bug in do_get_array_slice related to snapping

33dfce9

oskooi changed the title ~~collapse empty dimensions of grid slices by default and remove snap argument~~ collapse empty dimensions of time-domain grid slices by default and remove snap argument from get_array_metadata Dec 30, 2020