From adeedf1bd25c1655abeaa1181db6763371de5f67 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 9 Nov 2017 20:39:38 -0600 Subject: [PATCH 1/9] f --- python/examples/cavity_arrayslice.py | 71 ++++++++-------------- python/simulation.py | 18 +++++- python/tests/cavity_arrayslice.py | 88 ++++++++++++++++++++++++++++ 3 files changed, 129 insertions(+), 48 deletions(-) create mode 100644 python/tests/cavity_arrayslice.py diff --git a/python/examples/cavity_arrayslice.py b/python/examples/cavity_arrayslice.py index aed6fb8be..b81cbf62c 100644 --- a/python/examples/cavity_arrayslice.py +++ b/python/examples/cavity_arrayslice.py @@ -1,11 +1,8 @@ -import unittest import meep as mp import numpy as np import matplotlib.pyplot as plt -################################################## # set up the geometry -################################################## eps = 13 w = 1.2 r = 0.36 @@ -33,60 +30,44 @@ geometry.append(mp.Cylinder(r, center=mp.Vector3(d / -2 - i))) sim = mp.Simulation(cell_size=cell, - geometry=geometry, - sources=[], - boundary_layers=[mp.PML(dpml)], - resolution=20) + geometry=geometry, + sources=[], + boundary_layers=[mp.PML(dpml)], + resolution=20) -################################################## -# add sources -################################################## +# add sources sim.sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), - mp.Hz, mp.Vector3())] + mp.Hz, mp.Vector3())] -#sim.symmetries = [mp.Mirror(mp.Y, phase=-1), mp.Mirror(mp.X, phase=-1)] - -################################################## # run until sources are finished (and no later) -################################################## -sim._run_sources_until(0,[]) +sim._run_sources_until(0, []) -################################################## # get 1D and 2D array slices -################################################## -xMin = -0.25*sx; -xMax = +0.25*sx; -yMin = -0.15*sy; -yMax = +0.15*sy; +xMin = -0.25 * sx +xMax = +0.25 * sx +yMin = -0.15 * sy +yMax = +0.15 * sy # 1D slice of Hz data -dims1d=np.zeros(3,dtype=np.int32) -v1d=mp.volume( mp.vec(xMin, 0.0), mp.vec(xMax, 0.0) ) -rank1d = sim.fields.get_array_slice_dimensions(v1d, dims1d); -NX1 = dims1d[0]; -slice1d = np.zeros(NX1, dtype=np.double); -sim.fields.get_array_slice(v1d, mp.Hz, slice1d); +size_1d = mp.Vector3(xMax - xMin) +center_1d = mp.Vector3((xMin + xMax) / 2) +slice1d = sim.get_array(center=center_1d, size=size_1d, component=mp.Hz) # 2D slice of Hz data -dims2d=np.zeros(3,dtype=np.int32) -v2d=mp.volume( mp.vec(xMin, yMin), mp.vec(xMax, yMax) ) -rank2d = sim.fields.get_array_slice_dimensions(v2d, dims2d); -NX2 = dims2d[0]; -NY2 = dims2d[1]; -N2 = NX2*NY2; -slice2d = np.zeros(N2, dtype=np.double); -sim.fields.get_array_slice(v2d, mp.Hz, slice2d); +size_2d = mp.Vector3(xMax - xMin, yMax - yMin) +center_2d = mp.Vector3((xMin + xMax) / 2, (yMin + yMax) / 2) +slice2d = sim.get_array(center=center_2d, size=size_2d, component=mp.Hz) # plot 1D slice -plt.subplot(1,2,1); -x1d=np.linspace(xMin, xMax, dims1d[0]); -plt.plot(x1d, slice1d); +plt.subplot(1, 2, 1) +x1d = np.linspace(xMin, xMax, len(slice1d)) +plt.plot(x1d, slice1d) # plot 2D slice -plt.subplot(1,2,2); -dy = (yMax-yMin) / dims2d[1]; -dx = (xMax-xMin) / dims2d[0]; -(x2d,y2d)=np.mgrid[ slice(xMin, xMax, dx), slice(yMin, yMax, dy)]; -plt.contourf( x2d, y2d, np.reshape(slice2d, (dims2d[0], dims2d[1])) ) -plt.colorbar(); +plt.subplot(1, 2, 2) +dy = (yMax - yMin) / sim.dim_sizes[1] +dx = (xMax - xMin) / sim.dim_sizes[0] +(x2d, y2d) = np.mgrid[slice(xMin, xMax, dx), slice(yMin, yMax, dy)] +plt.contourf(x2d, y2d, np.reshape(slice2d, (sim.dim_sizes[0], sim.dim_sizes[1]))) +plt.colorbar() plt.show() diff --git a/python/simulation.py b/python/simulation.py index a8452a381..d2d962631 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -652,6 +652,18 @@ def output_components(self, fname, *components): if self.output_append_h5 is None: self.output_h5_hook(self.fields.h5file_name(fname, self._get_filename_prefix(), True)) + def get_array(self, center, size, component=mp.Ez): + self.dim_sizes = np.zeros(3, dtype=np.int32) + vol = Volume(center, size=size, dims=self.dimensions) + rank = self.fields.get_array_slice_dimensions(vol.swigobj, self.dim_sizes) + sizes = self.dim_sizes.copy() + sizes[rank:] = 1 + n = np.prod(sizes) + arr = np.zeros(n, dtype=np.double) + self.fields.get_array_slice(vol.swigobj, component, arr) + + return arr + def change_k_point(self, k): self.k_point = k @@ -680,11 +692,11 @@ def run(self, *step_funcs, **kwargs): if kwargs: raise ValueError("Unrecognized keyword arguments: {}".format(kwargs.keys())) - if until_after_sources: + if until_after_sources is not None: self._run_sources_until(until_after_sources, step_funcs) - elif k_points: + elif k_points is not None: return self._run_k_points(k_points, *step_funcs) - elif until: + elif until is not None: self._run_until(until, step_funcs) else: raise ValueError("Invalid run configuration") diff --git a/python/tests/cavity_arrayslice.py b/python/tests/cavity_arrayslice.py new file mode 100644 index 000000000..41d6120e0 --- /dev/null +++ b/python/tests/cavity_arrayslice.py @@ -0,0 +1,88 @@ +from __future__ import division + +import unittest + +import meep as mp +import numpy as np + + +class TestCavityArraySlice(unittest.TestCase): + + def setUp(self): + + r = 0.36 + d = 1.4 + sy = 6 + pad = 2 + dpml = 1 + sx = (2 * (pad + dpml + 3)) + d - 1 + + cell = mp.Vector3(sx, sy, 0) + + blk = mp.Block(size=mp.Vector3(1e20, 1.2, 1e20), + material=mp.Medium(epsilon=13)) + + geometry = [blk] + + for i in range(3): + geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i))) + + for i in range(3): + geometry.append(mp.Cylinder(r, center=mp.Vector3(d / -2 - i))) + + sources = [mp.Source(mp.GaussianSource(0.25, fwidth=0.2), mp.Hz, mp.Vector3())] + + self.sim = mp.Simulation( + cell_size=cell, + geometry=geometry, + sources=sources, + boundary_layers=[mp.PML(dpml)], + resolution=20 + ) + + self.x_min = -0.25 * sx + self.x_max = +0.25 * sx + self.y_min = -0.15 * sy + self.y_max = +0.15 * sy + + def test_1d_slice(self): + self.sim.run(until_after_sources=0) + + # Low level 1D slice of Hz data + dims1d = np.zeros(3, dtype=np.int32) + v1d = mp.volume(mp.vec(self.x_min, 0.0), mp.vec(self.x_max, 0.0)) + self.sim.fields.get_array_slice_dimensions(v1d, dims1d) + NX1 = dims1d[0] + slice1d = np.zeros(NX1, dtype=np.double) + self.sim.fields.get_array_slice(v1d, mp.Hz, slice1d) + + # High level equivalent + size = mp.Vector3(self.x_max - self.x_min) + center = mp.Vector3((self.x_min + self.x_max) / 2) + hl_slice1d = self.sim.get_array(center=center, size=size, component=mp.Hz) + + np.testing.assert_allclose(slice1d, hl_slice1d) + + def test_2d_slice(self): + self.sim.run(until_after_sources=0) + + # Low level 2D slice of Hz data + dims2d = np.zeros(3, dtype=np.int32) + v2d = mp.volume(mp.vec(self.x_min, self.y_min), mp.vec(self.x_max, self.y_max)) + self.sim.fields.get_array_slice_dimensions(v2d, dims2d) + NX2 = dims2d[0] + NY2 = dims2d[1] + N2 = NX2 * NY2 + slice2d = np.zeros(N2, dtype=np.double) + self.sim.fields.get_array_slice(v2d, mp.Hz, slice2d) + + # High level equivalent + size = mp.Vector3(self.x_max - self.x_min, self.y_max - self.y_min) + center = mp.Vector3((self.x_min + self.x_max) / 2, (self.y_min + self.y_max) / 2) + hl_slice2d = self.sim.get_array(center=center, size=size, component=mp.Hz) + + np.testing.assert_allclose(slice2d, hl_slice2d) + + +if __name__ == '__main__': + unittest.main() From 035196660e742d62762d9a1236d183c4bf776bcb Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sat, 7 Oct 2017 12:28:30 -0500 Subject: [PATCH 2/9] Add cavity_arrayslice test to Makefile --- python/Makefile.am | 1 + 1 file changed, 1 insertion(+) diff --git a/python/Makefile.am b/python/Makefile.am index e0637e6ed..029ebf99e 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -33,6 +33,7 @@ TESTS = \ $(TEST_DIR)/3rd_harm_1d.py \ $(TEST_DIR)/antenna_radiation.py \ $(TEST_DIR)/bend_flux.py \ + $(TEST_DIR)/cavity_arrayslice.py \ $(TEST_DIR)/cyl_ellipsoid.py \ $(TEST_DIR)/geom.py \ $(TEST_DIR)/holey_wvg_bands.py \ From 3ea0ab79bcbcceaea14a635f3cbca4c902c221f9 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sun, 8 Oct 2017 09:09:25 -0500 Subject: [PATCH 3/9] Improve readability of if condition in change_k_point --- python/simulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index d2d962631..0518ba650 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -668,9 +668,9 @@ def change_k_point(self, k): self.k_point = k if self.fields: - cond1 = not (not self.k_point or self.k_point == mp.Vector3()) + needs_complex_fields = not (not self.k_point or self.k_point == mp.Vector3()) - if cond1 and self.fields.is_real != 0: + if needs_complex_fields and self.fields.is_real: self.fields = None self._init_fields() else: From 6955a965e210a99d72c5af36115a33ffe79de08e Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 12 Oct 2017 19:25:05 -0500 Subject: [PATCH 4/9] Return ndarray instead of storing shape --- python/examples/cavity_arrayslice.py | 6 +++--- python/simulation.py | 12 +++++++++--- python/tests/cavity_arrayslice.py | 3 ++- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/python/examples/cavity_arrayslice.py b/python/examples/cavity_arrayslice.py index b81cbf62c..dd98159e5 100644 --- a/python/examples/cavity_arrayslice.py +++ b/python/examples/cavity_arrayslice.py @@ -65,9 +65,9 @@ # plot 2D slice plt.subplot(1, 2, 2) -dy = (yMax - yMin) / sim.dim_sizes[1] -dx = (xMax - xMin) / sim.dim_sizes[0] +dy = (yMax - yMin) / slice2d.shape[1] +dx = (xMax - xMin) / slice2d.shape[0] (x2d, y2d) = np.mgrid[slice(xMin, xMax, dx), slice(yMin, yMax, dy)] -plt.contourf(x2d, y2d, np.reshape(slice2d, (sim.dim_sizes[0], sim.dim_sizes[1]))) +plt.contourf(x2d, y2d, slice2d) plt.colorbar() plt.show() diff --git a/python/simulation.py b/python/simulation.py index 0518ba650..4a77df8e5 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -652,17 +652,23 @@ def output_components(self, fname, *components): if self.output_append_h5 is None: self.output_h5_hook(self.fields.h5file_name(fname, self._get_filename_prefix(), True)) - def get_array(self, center, size, component=mp.Ez): + def get_array(self, center, size, component=mp.Ez, arr=None): self.dim_sizes = np.zeros(3, dtype=np.int32) vol = Volume(center, size=size, dims=self.dimensions) rank = self.fields.get_array_slice_dimensions(vol.swigobj, self.dim_sizes) sizes = self.dim_sizes.copy() sizes[rank:] = 1 n = np.prod(sizes) - arr = np.zeros(n, dtype=np.double) + + if arr is not None: + # TODO(chogan): Check size + pass + else: + arr = np.zeros(n, dtype=np.float64) + self.fields.get_array_slice(vol.swigobj, component, arr) - return arr + return arr.reshape([s for s in self.dim_sizes if s != 0]) def change_k_point(self, k): self.k_point = k diff --git a/python/tests/cavity_arrayslice.py b/python/tests/cavity_arrayslice.py index 41d6120e0..25a7f6236 100644 --- a/python/tests/cavity_arrayslice.py +++ b/python/tests/cavity_arrayslice.py @@ -81,7 +81,8 @@ def test_2d_slice(self): center = mp.Vector3((self.x_min + self.x_max) / 2, (self.y_min + self.y_max) / 2) hl_slice2d = self.sim.get_array(center=center, size=size, component=mp.Hz) - np.testing.assert_allclose(slice2d, hl_slice2d) + np.testing.assert_allclose(slice2d.reshape(NX2, NY2), hl_slice2d) + np.testing.assert_allclose(slice2d, hl_slice2d.flatten()) if __name__ == '__main__': From 4cf3b9dfd1594a6093b22a16e7a6d6d86a265bfe Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sat, 14 Oct 2017 10:28:47 -0500 Subject: [PATCH 5/9] Adjust array tests --- python/tests/cavity_arrayslice.py | 46 +++++++++++---------- python/tests/data/cavity_arrayslice_1d.npy | Bin 0 -> 1088 bytes python/tests/data/cavity_arrayslice_2d.npy | Bin 0 -> 38384 bytes 3 files changed, 25 insertions(+), 21 deletions(-) create mode 100644 python/tests/data/cavity_arrayslice_1d.npy create mode 100644 python/tests/data/cavity_arrayslice_2d.npy diff --git a/python/tests/cavity_arrayslice.py b/python/tests/cavity_arrayslice.py index 25a7f6236..a4e399f3b 100644 --- a/python/tests/cavity_arrayslice.py +++ b/python/tests/cavity_arrayslice.py @@ -1,5 +1,6 @@ from __future__ import division +import os import unittest import meep as mp @@ -8,6 +9,8 @@ class TestCavityArraySlice(unittest.TestCase): + data_dir = os.path.abspath(os.path.realpath(os.path.join(__file__, '..', 'data'))) + def setUp(self): r = 0.36 @@ -49,40 +52,41 @@ def test_1d_slice(self): self.sim.run(until_after_sources=0) # Low level 1D slice of Hz data - dims1d = np.zeros(3, dtype=np.int32) - v1d = mp.volume(mp.vec(self.x_min, 0.0), mp.vec(self.x_max, 0.0)) - self.sim.fields.get_array_slice_dimensions(v1d, dims1d) - NX1 = dims1d[0] - slice1d = np.zeros(NX1, dtype=np.double) - self.sim.fields.get_array_slice(v1d, mp.Hz, slice1d) - - # High level equivalent + # dims1d = np.zeros(3, dtype=np.int32) + # v1d = mp.volume(mp.vec(self.x_min, 0.0), mp.vec(self.x_max, 0.0)) + # self.sim.fields.get_array_slice_dimensions(v1d, dims1d) + # NX1 = dims1d[0] + # slice1d = np.zeros(NX1, dtype=np.double) + # self.sim.fields.get_array_slice(v1d, mp.Hz, slice1d) + + expected = np.load(os.path.join(self.data_dir, 'cavity_arrayslice_1d.npy')) + size = mp.Vector3(self.x_max - self.x_min) center = mp.Vector3((self.x_min + self.x_max) / 2) hl_slice1d = self.sim.get_array(center=center, size=size, component=mp.Hz) - np.testing.assert_allclose(slice1d, hl_slice1d) + np.testing.assert_allclose(expected, hl_slice1d) def test_2d_slice(self): self.sim.run(until_after_sources=0) # Low level 2D slice of Hz data - dims2d = np.zeros(3, dtype=np.int32) - v2d = mp.volume(mp.vec(self.x_min, self.y_min), mp.vec(self.x_max, self.y_max)) - self.sim.fields.get_array_slice_dimensions(v2d, dims2d) - NX2 = dims2d[0] - NY2 = dims2d[1] - N2 = NX2 * NY2 - slice2d = np.zeros(N2, dtype=np.double) - self.sim.fields.get_array_slice(v2d, mp.Hz, slice2d) - - # High level equivalent + # dims2d = np.zeros(3, dtype=np.int32) + # v2d = mp.volume(mp.vec(self.x_min, self.y_min), mp.vec(self.x_max, self.y_max)) + # self.sim.fields.get_array_slice_dimensions(v2d, dims2d) + # NX2 = dims2d[0] + # NY2 = dims2d[1] + # N2 = NX2 * NY2 + # slice2d = np.zeros(N2, dtype=np.double) + # self.sim.fields.get_array_slice(v2d, mp.Hz, slice2d) + + expected = np.load(os.path.join(self.data_dir, 'cavity_arrayslice_2d.npy')) + size = mp.Vector3(self.x_max - self.x_min, self.y_max - self.y_min) center = mp.Vector3((self.x_min + self.x_max) / 2, (self.y_min + self.y_max) / 2) hl_slice2d = self.sim.get_array(center=center, size=size, component=mp.Hz) - np.testing.assert_allclose(slice2d.reshape(NX2, NY2), hl_slice2d) - np.testing.assert_allclose(slice2d, hl_slice2d.flatten()) + np.testing.assert_allclose(expected, hl_slice2d) if __name__ == '__main__': diff --git a/python/tests/data/cavity_arrayslice_1d.npy b/python/tests/data/cavity_arrayslice_1d.npy new file mode 100644 index 0000000000000000000000000000000000000000..056caeb073db63a09b165c1991f6d32a7ea0180c GIT binary patch literal 1088 zcmX|=e=L-79LN2*Gs05LkJco+b1W-oA@0`qnZ_oaaj{VO(cQ_e{5o~o$gf25<3bHL zqQm@XoTytU-`jMh^#`jXRO9YRml?+$o0(dgPd)!Ud%k;~=ll7-->+9)L<}b~&eSB- zB$vS#Wb;IfAU6g(0NYwgT z9#z>+qURpD8!sM{Xgc#sciSBjy$#C0yu~F^)~45&Xd4n$b1B(LaxD`62=39bv}m+V z_cCxmgS3j-Fi!^!ij#`xLar{L2@U&vTJsEY_|&L>x7dLARSLSGOpUdYE^1%qckI<5 zE>DYJ!Hn_!CogaeSeVYKcMdY*u_bA7aD)+8R5kP$bBy?yzGLV_fDso6M;2+;MjYbP5AMo9fV0^Zb_5;?XiPP4>>$WY-;&B|}Y zUa+v8?%Y9q*-=jeU5})^?bhW2Bd(DIPWskZDQzlqs2|Jil{+lCUF3qZa^Rjix1sj1 zk}2-2SBhPgfpG7G^)+~p0^ha6e@fuLL9h=C_RWJiZoynlFy}nXt%Gwg;avW3&VT2g zf*eXA7Z=Eh4!Jc!j`onNE#zzpx$B{aF6c!KJyk<*ub{`L(CaYt>;k_v9oDDi0sm&{L@oV;j$#^C?VtDL%EVK6ptVRXvWVC>(^*2!*^kzTK` zTt@n)H1j`^E`#AgqntxMk0$mG{3p@ypJs{-C9`p zT(jj zN7JgJv5LH$diUy%#9Fexu}|g9-g@#V|IVj_3mQo5WwSqVkDAD-MGNi=PiQ8_7*jCs ziFF&WFA8yz5LX>}Y>;;aaBK#yMBuEHV{p5H$FMZR>tiOvQxm-Tpu-L5@(?)ZdfCrZ?iIDwPopJtg}B$10Lr*|6)7Lmlb!t!~}F z`4fqqG}}z=P!m~jquoMXt%Wq@3a9iYwvc&_D~e>=TgW!ulnpHdEo3J~Y0Pu5ZWi{% zAkGXwChmLWkwf0Sz`+At!;cu8BET)ZhvC7}XLwnG=OXa-f)2}|OA2&)^x0GJnTQ+J zS=h5^%MvHb-}E4$a;+Wxa&6o1+8LIV=ef3N_C-}H=^b#|b&Umi{OCsTsJA^iG2w^Z z{&|ta_)%D!^0sWUY)eP7#_tl+*I)j${Yn)XFi^ zMV@Ztoi@PWm3oHjZMq9|e#YP?h+c@Oh zyoGz>4|0h7)Ja#QPjQHKr-aAQc@7csJlMO=oI}Q8yovdZSa$*Y91zC_arYsQF7ieJ zhbeFg0q0-fW`l<>cga zPHhzy8wj6KAqtb?!Yy}QCytXPS3KppPkf!1__?2ZNb>8_7toS;!v0YmA1hQqOrPYQ znEJ7j%<@mXb=a()L~W~I{MD+7XfF94m)qDvE*}VgmY~KVqdSwut4?#s_{L+6iB~w} z$%OxOCG5xa`crXhog0TVr=(aqcymYx#%q{=igm*c%)Va4F-6>M$nzd~Wr5=?aG3#T zDsbNi4=yh^@MLXdcuPTt+t5W1ItlM}30~{$K#iJqxLh1MObQe5 z;gHlN*9YRBaY!u2ADA!0I$P{pjW{O|cL;eVAg?iS+yO3O;1mXK9`LXPFH`UoUd!;- zQDSsZ@@I6}dz#UyC?fYl&Kn0hIQ2!XKg))?6&znRI^Kvb6LvTCd2c}ST&L~TooGjn zzEs;G{?MH~D_ghMRq`o$_}3t<5xru+wyA-5Bos?^ z9Bm3`TgQIBN=LQpv^m7pHTa7ekS z{;L;u98!&O0Q1kW?jiQMBF=rpl|i1{$lC}U;=uI`ILm?C5jnwN*3o*PMix?fs zr!l%5i)3{AY%@dU^P78A@A*u;9vTge zIO<9UKGnMa5Pm|IcCIQhdX_|X?A>5KtW``L?d8rYAE+T3FVE*Jy4OGwCGSj--_SzT zOqIlVhFgfs=9c}Nr*lZ|pOe!rt>lo(=CvYeV}A!d+}-$Z?c$Kf$u$Bp8XWRY(fY)o z*4R3X-k5)abtTwWfH)qA%gqysya$0J0l1a`XF~^r`!#sH0WT}?bUeWD4ulR_&}FSH zqZ2th$@x+5eR}kPY}OafUD}wR;-V0BjVgco&r(NThkmJ)aNJR7PLxZ3%k^tIl0D7r zH*2POlf~Irt0%N3kn#?R*wV^EVv_B4sqJ7jnP1}Kvf{r_B=YU#2+53QqLlEe#yz8j zmFW4~lY#8ODv6 zkHb3gE6l!N#JP^R+&pQ>I}eg&KD)#!A+Ku&>DO>EdaYXex;`O}a#`N_!NVjgyQj;4l z#5k_V{LQR5@>@PceV$rAdDyi{VueQ~QQ3MebB{_rX}yp+TkuK~;cu5cy8KiNx${t~ z(#@fT@agT{mGrQMY~quMYYk{2T3^#heaP5o(h{BeriDDz;ohW=SzR$t~M3 zwe7Me4cI5XQFzpiN<9j{`@GJXPMCH@!++Bq8qDjIH|l+vsz!^@bGeG-o@W^Uab;`5 zwsDBIQj8)$tp~}D;9MekH+%ndp>krLE;MAIP)n{RMZLC7_(bykXVU!0W-{qD%l5<4 z7SblxnDuj23)!=5^uJ2k7Sf_{&GpOJ@2%n5143q-T8M7j_OSNNEyNh(CCq=tI$iAR zLY!%cTa7$>k=Ft^jsjOLa4xrHaBF!nJS4LiUS8l?G?n3f20C~`mjlo#g0-hpxjgrb0=u!=J@k3 zuZNOOx#cTLYchyx^{4r_wiJ<|2 zajT=~Z8JI2a#OBtZ1gHg`>!snnY?W?UicuqnWVjFPuvvUOwut%V%`_)w6U)baeg9h zAoA=&-Y(!+1zc}{Gen8O&ELTA*bH8};F${EEa+eYU5cR7_r9}%nrnRMo{0JnwrWr4 z3#T80DZ-EGl+iM^+#_x<~$C(haA%1-x)IQb%C$_Za###2EQTP*gC53eCwm!3{OZ}f@S-Cb>SBBqI) zAHK3scxp2-yp^f{m#3MuS}pDk9DC0G9ACG4+1Pvw#?P47!#XqU#bS7!1<&jehPN|xP=hX=(8+dHVusjfKRUs0&2-ToKPuL3qnVTI zN81R${I>&sbnQyhdHgqgskCpE^_|&I>72x668rg|(6y%LI-fNBNA-_hDOeN|NJiJ& z`Nr|bk|>8o*=YwT$u3p2UQ>}n{CJX2zpE%Dl4ow*s&XzTr>tALPg&Lwo2^f-ON{*; zwKA!HwM)2(%+ygabvn>Qj`~=B*lE~Aj$o9;{8g+whJCXUC;l`OcRljlMP4i5=$2z} zIUHkfu07A-mH-c7@R9}3bnuph4pX5^Aat^MP}b;O8bC8gMEAF7Kc{ym|Ed@K@SLtb zlan?1Kp<_uy^7a2F_8WV{XAQ8WbEI8jlNk^UeFme(;9jspVOiH-6J8zPifNu?NZ+i zH>%hFr{-_o*=JCS#vH;#4B84D!@8F?sbKGdPB|7+moW z7@RG@JqR8muNYp<;3*H@X3*gsbh!te>aFDcraTR#L&c9;cZ9s4ig|sn_lWU%S4#x4Pwj1)SQzj3aAObk6OK5P5- z=s+4&zu?dLRhP)?$|v97eSS=)Rxsn*56t{Utn(;f_7x(IHR5(5k2La{0LK>K(gV(+ zQw;7<@VEqC6Tx#7ytz7zFJW}yyUOUaSM`d?t4%>PsJLmrQ{hYcp!4!GL)Q@c+IS61 z`OYg^ZEx0b@6BskdM4htLpYp%oZ{dXus4!!*wt}!(O?o)Z5Mia@nQ}QGtDuQPNgo4fJxwmyb&=n(ncjPMN9oXuR%&wp zzRckz?etaptn!Z%9rR18hSC)4PP%c$@hqLYU+HFb(%yNZi>|Kl&3B&ojmkZ?y!YSH zZ*)0EY0S%G-7@SOK%5H1H9{UI(M|fT&72tFGv~+e1wcDzHTjE_D{k0>{ z?|5w@Re5ruvM4Q&{uK@`={Z_Ml_DyAeq=P!B{I)1?^gay*SN?pPMqFBcbuc?bt}5) zmsR?#|L%6vkPNy;M7oE53FkeeNqeZ?e8*V}4SJ~{5eZ;P_R;luljIJo_tB*oB`|*k z>%L-NIO235?gHdliM%Y}kO8i4;0$}r;Qj_4_rQw}JZFLTWaw}hy7WLNtA!!YYehro zXw4~ab%QY4a{k+zSj}+ipey?Ol35fDxN&cq^vqZ)=lob}oJ<1cInZYplbB4Kr6aHI zjVPe|HCH>xJo`X5lj}j+>p8SWr_t;DnJ@HO(1CV=lrAbEdehEIyoc`iydgAIvzIQv zW7{vdu8*pC+G}L>_0dMBis=q>`l-L7vgxMv{nXS%{gI?%KQ+eq4D-!cr-6Oah@*hG z{K#X4yzhZyBXAi3r{a7Dw*+`_c|8PAYw*s74ujC;ICT2BER6q@K` z);#!kS2%t5GfH~K%P6YmaKmAiM=Vt|llT?-JAujx1#R2hl|r|hDwQ7lUPLckS6F^X z?h_r{mpsX3!Do6qskbepyp!HkNH^5I{GEF5t_gg1xt9vxQ*GZB*hhV3g87$B?WZrt zo$x5y&`%G(UfFPKYd@`7Je8BK+D}tEd~R>w*H5!BaxkxfbrIO-hd2U=>w!G($mFo=DWpWFX>U68+TsTgwpEM>r)=-yrqvH ztQXvTD2mS7n(_O<?@Qo_2d+iH`C%J_ zI}JQOgO><+x`1~!bdZ5AM$n1n+U{ceE0ESXc7E(A3#Mu5>ZjITd`%l?XI-x)a0<;d+@P%}UP_zvS(`j(H&OlLR~CISZKst|+vYb1 zeWT~=wZy%aj{RKmA4@KK-AB!z*#2U5_0hEGn=jac{j}}Zlr1f@`YA_wFk#yAemX7Z zYMIrTep26GCasb0rzz3;OYdO|4uNGH-^`U)hUw*N(oWFLV12T-g;(C)^QJ z&0Zc)Cz<)DPfbatM#V|1ZTL&+h23sex(Ax5@$4B5_nO-2#(RfOUSxICeMyC$b$5EH zZrt}UfA2oZk-KrEF1?T5__O0!{f9o9Vlsct((ip#lf9o9|LvoDrn;Z-9n;anIOX%A z3H{Ut<2KBzVO<*bJwTlAh^vo0TaZ^BIJN+n0dRT(_dM{J0bXyxvlhG!p~Exi@(?=h z@bQk({Qj73DhklqHzj}`v^XU@ay^KO-Vc-GRKKDpJ8$(oUNrV|7O>scW+aBjiSKGH zaZRN8MQ413W|q+3d1cpvS2a;<_S1WX)$R1FMv=p&qun${#BW#n>t4F*;*Fibr~Bxm z#qGDF?fPgTjqeWg?W07zS0zI~wkeoQ}%`k3E@ zb*k7`ggBvyyBc}Ik+%pqP6O9Z;4}p8D)2Z3URL0#0p7=KYp$s~x5@TMIva0{c3eH*2DMj~m|ytgXWW^r`ktuL;2 z%0=`|@~ZL`I~%F_jPCoLD(!SYOYdVw-8ZUVGIgDkWiKs%w|slU>^?eVIi8dV z%$hE=r;n~wS(P&FSRc)Na#{by**?lY{pX{3CVljluaj1Z+1Sq^#^Gtqd?MD(!9H%B z--ug@Jlwp)z;Ob&b^zys+YD}gX@*A>cq#BPJb!|BHgvFsF51vZ^HN=qGv5^|I(y1u z`*eGneR+v$fwKqA-rXEk?EaKqh^`M+IrDCbf$Okt%hONxigl+xPXxw3pggguGZ$*-Q7o_g^yD z+DrfBvOj$4?WMoS+QsYt^wO-An$a4(ee}V{oUEh#ee@1S?)(p|8^pe&h+~boZOG$? zyi0*21-KG{^LZVEI|n?1z$*YerPLVSgV5mybWym>=rr)MFiNuH3Tg71d^msG5h7XU z>c4LK*m+0P96mSdKzZ_w4Gq1$>71oWb3DETP#rJXbAKH|>2e+S{ac>qQlo{fIa~LA zpy_#=-|gGkN~OgdWj+;j(&{GOXpQ|nG^WTbly6@zP5Jrd#lV#4gc@BakFhNjllQ;^J!Se!M=ltQ;4{+$dik_zeO1w zO~924oYxOCxMjd&4|th^r)mtt+q8nwfo;#|Vid&al(9TX#FwJ z!;dIuJBa@oeL_~*2q%jqdXY1)7uIj{aU!-e+>-bc&ynjJIZI!L$I}I`cb{dyE~UA< z=RdO!ucsPb=QATUIP{jF=@jv}FVuV0Ouhuo@6@xyu3~$>8MhEvuMwgS&Y5EnbhneLeq@n%op>x|q$(`4{GVN-i z#L$=bbWcGDSbeR<$H(Nmw07cumV&e^l?#+~FWD*N8^ ziP_e2I(@@Ri=6HC^mAU5)&%a?Z7N zQ@XQcuXJ5Et=)O2$Gx$e#%qPC1-6g9hPqz4;V%ee(&AgfmV@3(Yf4s;=Cv# zyJKX*+b2DCbxq&+TDcN0-(?*v}wr%GqbWy+6L18XR-{^Od zM~3Yp->9C8cG2MCZ}e4i#9-CtZ!`zvWz7G?x?t?vgE;nxtB*Y6gP6SSTNxb7W-z$c zhA=pn1NVLKcn4mKz;kpv!+RohsPkfUq0osAKkE~}9YZ?$&j}=ciY02j;sFD)@noIC z$wjBE5=e`wQHI0wL=qh6z-jhNB41o&!we56lb-fD@e#SHWYdG`fjTbfBuw@5-LB{q z;vp1O(Y8ODu+J3^-7h{#Y?mq4Hz%ji`ct=yy8o2Y(A13jqdz{;3ys$b9WQ;Mx5YQF zx_P^k-sczibxQRs^?Twv!{^3V>W0w_^Ojg=gMHzM(}cJ!$Rqul$?FUpI>5CRI0u2d z96ap6D-k>=gZD=0P!3&QL8o_dEY>5%1QM_3@JL@aiFg+*kQ z#$CTjNk*grl`UnHoO!bTFFs@tIqzp_Wiff=&nIQR3k_vt$>+(FrGHkD;f20ad1lrU z-!1DZ&sNlsDVB>)J}s{#mwzmfwc;rzxAy4u{;5bNq5PM>hV2NY!)b@LN>-H68vzey zu=h03rue$Een*>W1jb+f%)C)Lvu;F~*{8afiSvDciR*+sACR{aIMxAIByhR_cRF~) zg4b2>{13cuLx(reB?mf<-z<37Fe91Fi}`7?P&9lz-1l1f4`Dq;RE)-A-o*N7vBxVFfXjJyWG(F9zbzf+(aCf~bEB78J6$&|s*N6I&5le1ZV=XdVTCANRc^)#LH zNWIz*nOdF#(kI53P})>Ywx?g7#P;Hl9(BGs=d(IViaX z6GmmscVgY&*UY{>h+~YnYRHp?yy?Kfw~xW42b{6M9RnVxz)KrE$Afo2bdZEDlEsWp zv8N}7Tsh4q^P(5Krf$n556S|zZ$Fq#^hQ!@uLb0ikG#iZp2)l>b0gO-{WMxY2Kn_( zCNC@@?tb<%1_kxRa5YcaOW!Z#bm{fkc5`~k3X!%0+t&Ude!3NZe;EvtJG&GtB%_B& zxTNi{z`SAd`kI8vxZ+{b@g{+j(l<ppFs}w zc`2=VnMLj_zVTgDAeXe2+?+E#Jdfl!jO*8YRY1adUiuv4Eg`$|^Y1yHDkC}_jy=ti zO(f%vq4AD|-^h5)r6=1g2Z;TWyAyue{U(Cyj>7dGf6413U!^kCcv$=N9TV&K^0F=o z>1-)$eFPo#u_Js0uuiN&G)d$x9B; zCM!E){RTDPk%Wb|t-HtSZer8qzD^A;B2U|7K9?nyl1j%4!LtqJ#3xC)X52tC35jd9 zKYFu=NS&GFv$l;NeQi8nfHhI=W2Th7083nnbHhGJfF+Dk4D%bY&J6pE5oZqKa`Rk7 zUM`M_!1V|?!+^UGJiNi{Ie12aw>)%s3tg0;)0w|-pGlNu5{qxIr+wDUA%Ep{JEknl zBW1O1G5hojh*9I(v{um)GQTm9<2ZI7W0!PN-1grUMEY*GldIua9iX>zi!)yz`S9VY z%)zmrEBh}?&L@5zC7wIKieA0K$LidCBRe%8)USC>>B7GO10Ijy?XBEVYLo^7i? zPmqImW`3v0oxq%WX`Mi-)`g)uY)-` z7uyPm!WGBkFA7S?PRAz~+~$^(WJj~+?BGfwYTr~pY}`tYJ_;5RpWRQcm|fVlbJicy zspqvTj(2ta`UmkR?52vUt&ELDp)FwwND}bx*NxCE~n9+^@*90eSm@;{Lcfk7vbl3}Bu0bcymdQ3}?v35^8L#XyO)iHtb)Fk(SI!%&^IMyxoL@j3 zb7HM!D@%xg>I#kJisj^-eqTgOLnRUZyYNiUpH?C{Z0htPub&9lo~hNCF-&r;3m1yr z;AIu{In3VA&(Bh6KC&s^PJnf)Z_w&ng#as0=dF%Lw*aeOxa^1gI6;=_ESEDg#Q%?8 z&BLm*1X)s5uEw=<1z8I)-opF^tow<59}vd?akn8)9r7jvM^g=hOBOi$h8f&P!Q=Zb zhLK$1$ZyDb@`!i8+nITdjK0u=WY!Zoe9wryx$XWbIeBe1A zyxpL~eCTotI=ye5UcbwlO$5AGyjnLUizF}O9CDnJOWxE*JZW%!KUP24Ci={`h|He+ zqUG$8GSa+a`cUSF3KBlFsNSlul|0D09~C-0_B~4Kv8VIBVe&Y}=dWKjFYB0=t6vu3 zXRYANxHM_504rhHyiI#|39xuH3@Y~<2(YHdJXsZUO@Q@I{CCWSdjhNnrxT9Gc?ht| zryLXBFs5q(MikmnHcJ_C*h;Bo`bhrpc=9wOk?agyQr47^`K zhkh?cmsIHV)A~l~fOi_n)?Gf`P9cLB9M#DYoRUq1(#DOOc5dwFOrdtrbbkSvzA~h) zv#EsSjFas)dsI$p1<%yzPG}`ABERnpJM@#Qmjz32P8ud(|MnSda^Yn;+BB8#78_5Iyy$u@pgdf0|VqC5gD_w_+W!&3xUb0ivo^8k6Pfx`p1+<-F`xZ}ZNcm=~t6g>IB`yg~kgDyNa zj83wLvP5O}C6P#5ag|^3X=F$wPw)l2gvcyQzC=>*aQL#&<5aSAl4z8?XazmI`)3?W{=lsxY(9S0>no3{FAMRr)=ytSl)5H z`HAc((W&EU?k1WMeAa&^{UUG7FALa z4ma*%y&Z&a-(KdW6DYj54+m?^j zIJ@`V&#`}><4ssMZfyKDtH51iY-~v0{POPD*p0CR^9ESghJAApham1Q|}H@n&wK(>H&UUB?Cn3rT62$Cqrb?m+WnK#d%qeY-G3hMDVg$D!zv%2l2ALuUNm- z{V6YN>SkV#AunE*-SRp0hdp^&o*S3QXOE3O7_Bis0qeT4PX%$FAg&RKo>86f7S#H$JzG5ct`&+art6;ly)>-0rKCY|ZI-RgyIG#=suO>ko z{w_b#)w2M^1r+CY8kYaZ6za`Tci zzC5gwj_TUw4|rIzDytr*+~;ALi{+^Mx$>~CU{t_-5!TJbK10O$h`4)@XEO2%=`uKW z0oO<%gOksg!96*R;h_XxZ+#e^5epdJUyB$WO4l;Ftb|URzAT(7VDN^fo>#Jp^$n*z z^S61G9EzmA28*g><)f&drp7j6D_tP_ zF_;W!ruVsxy)UaQ>-}*izLK=FU)0D-w2UC_I&F|&wV^6Cj+)d&*BZRmx3ik`{Nd@{YS~E&-1cku z^$n5)jULx0SN{;5lUZLSjfcsKCPA_4kYS>P@hIk#u&x37J|fN>#N`qF?|*q@k+)Bk z!O2j5lEu=0N+e(8hYw14GKvmuepXqhU zInCF~p( zeTIME(-EEtj%8Qvof%AY4=_G~&2+f($>x2saT>~UwEPo9(KWxtvp$mz4>VJj61kDMDfLiH*t zKP0sNrZ+rtR_eL@pw~_h1$?#Xr#=`RFnC$;5C<*;VJcq;r#_Vb0Hl2A@ zJs-!h)om9`ZCfG4o<3LnWVy$9HYYzOO@7h@_PlTQ(SC0xu;)h)Ex&z#0{f|Ii}<=0 z5%!Er2VBL(MA#2qHNR^yvuIa7mHg56y|}!XR{g$~~#)8!s!wfJFkb%K&#q{8SI@oekozblI$Z&QcnfaB-uYV zM@-gil4K8KJdSzpIu7<}AdWfWrXbHrg%w&Ie{}AQ&dluU?{MEbOve|5XmFQzT4oR`!I#ynP>>}0|zM%{z&R4{Jj64>|JI;Z@aR|7^P8a@%GYGi3JT$;d8$6GJ_jl-^1zny( zCnXonth-?aG%~o_vAnE=_9i80R~#s(1$@62-CkNnV~luO!mifR!u5L>85(?~GKbB- z=-W5ZpC9P$+Oayf_^9((iv8SN=!K5^9QN|0|GM7F%w;QDhCE_h&SkH| zD2MqMSoaqDybxyuao-|OI`VRH+ypLT;A{i#)!>m0Uh&{r1KtA3j1IZbMHf0LV+_Xp zAoQxhzBt6mLtI(p=|Wz8;J6Q5_RuT(AH9BqmpgcdfOjBtSO8sE&}kFKMVR-;x?t>k zj5z&>%l+>XhrIWJ;|_4`15Unr=D#!dzdKj2Xz=9n{tX?tKS!$2$qheeviP~{#LuB7 zelAbo=X4)_Zo~2W=!BnZVf>sQ#m_x2eh;|6kB{(s@(sT?Bltb)!SB^I{GRFI_bvp# zhv67cV%`z!Iq$QZ4>c5ZV2C3rSUy)0N?wT;Co;dz8A{k zd!i1$H;UnVWIeuDO5uCv9DMI|!}rjc_+Dy<@2No;#W8<}yY8RwtL2D0^`GynyMQAf zxE2BDAK+dF9z6KIDh{5x;B5jO=0cYX(CHDr=jY>l|3^Frh~v3H7|#jbcy3sM=Llas zS198-gZtcZ6VD+McrJN>=adO}9ty{Ej3SGd3!#dvlTJ6 z;yIiyWuCXE;CcH3;t2OJ&+&=Kvmbd=f#VTyDFdglDf6B{l*PO^CA`qa5!&Qt%$65ASnK@Sdaw?@jjMeNMk5-si-e z|KEESJG^(7@)wOChzeMb?eVL$Vph{u|FZ`6*w zs=%=g@0p~4^AXKUmg$@?bWj%CSf%k=<@m_2M@5!?8-Yg&Q(PHsl zEe!A39Pz&JX%6$gFaz%kW$~U)CWCojn2Ps!+IX+G0`K|e;(cL0-UF`3dqHiyC%l1i zE9Q%_?l0aKu0kAh#BD&HPP{KP0FKMRHJHudoXla~7dC@OH+X#p&noa1T)@0Hb%ZV_ zp;Mnb^FF$L9rNBb4ez5V-beG}ee_-j=6!T2{r|p?cE)>NQ@r19;|ww;Xh^+RU7HBtoZXoTs(nJS`mONXGx1r)|Tz z(@&hIHQ+q05a($jIJbI+b1XK_wKn3MO9kg%691g1>EWD=XEt-5)`@epMvUKao@R)3 z)i_U^jX3g%tAIRjkv9%F{DJEea1H_YcktK`UgbDXy8_;yp+gsRc?g}daSnMH=aTw3 zr(A?{OB0-9is4*y2F^KaaqhVS=b-X97tO{wX%Eg#WpIv~j`LaZt+|?H6u%0-V zt-?8NJI2SDzl(Jy*yr=ldF@i<;pR;N4gy@Fz*z^}vf!~0ysm*Km-jj7un@W=K&Nn= zb64Qp`vcCwC*WM1drsbsdVn)HM_-L|b!VKjuf@5$KhEJ-;9R~M=k)t=Zl8j4{4Sj9 z%i^4$d+xs!bpY!yrlTG}0qX{^F9vaR5Vs0>mLe~=EeK zaC7Tse85W!JSU^h<_~m8g)VcT(+KLB5>dBv5_LR7sZ2eSB)857^-QU#1L{Cs&@JIX^H(gV~bNuf@u6m?6n|I{<(VtzW-Ibk2SF6tWMwj<9C zei%C#}jNCdpZoS-8 z?BmwmB_Zx26gKT=u{j1>Du(Ar4+Vs0-xs&V>%6(4`nUokpGEa?~9@M;+oM z)Fm!Moniy(7Js3RQ2}+0t5N58A9as;C* z#toQXfpwd)PYm^#&k@%Fd4!SoJa7a9*Cyc90&YIkm9~S|8t~-Sp}vI<8=#97bo#g6 zms{t026e9&PzO7Jx>z05$rhn*b_MEaf1$3n5Oua6PV4x- z@9T?tUlY_xA4T2tbJS6Lpsrd8b=LBzyN*B|w(LLkCo_OaKinJr+aUYf-^H7h&9cnTUHZ5*S4>uY`4Huk+$(a$J)>saJF37v zq{q0Iq>Ot?+PJsWf_qGpaj)qC?m0=|-cu6pLG8u8D1F?M`i^^3pK*`M3-_voanDK& zqbufHu&xLDt{~2F#C?N2Qpn52F%h_`fpax*6Y%f{FAARR;Oz_@&O(;^4({ch`{zEN z81C`$V${cdzOz_Yihc2j^Bi&Qk>?TessV=yaB*>_0{0K_5dY^spAmRJf(}a1r5HN> z!hO+YxEI-pdy?~UZ}KqiQJUdiWiIYnuAj)<7oCfHn166zGz<4M*|@iP1ot?P;$G(w z-18j8z0c3M2f73ILf_+_=y!}8G0(v|e%v$NfjC--yB~RCkvA4Ns(@<_aMt0zXexM| z1Fr+%IT5^jp~F<@!i)Q&;<%5z4EJu2;2!Q(+{-Yco63}H(k-4|cy~oYH*X@jZ-WIs`ZGn5>iMSUY zjeFwdUd(-W#RlfSy8`Z&=ir|C2HZP;g?s1^aNm77?z=A%V(zW0;J*7f+;^|TJ@-nC z^D%FWb@>9!zHT+$Xf~=yMc=U=UUu%4*-uK@Hz>eC&1ebIy}aG_jqn! z0Q3o{MBjkj=p!J6z5*0|1}>xTKq~qWRC4<%+-LeJ$nIwPDU1p-{S*QoGkp!xmNNYm zH1wH%3W8}&KZQ%^i{OMl31`ta!3lj7Trgh4yd&0K!#-E^QxHMiTj-}?TfpSCM?ZyL z^iw#DehPPi`x$uXftN3M_R2GTB$h&lbm-y_ojB-2(T=_p6VRt(BKmdALm!J=^tG@> zp9}wDre8;5Akzn9J^FS0M8A$k^y|2TJ{n2rt8o~8HUiLhqZ@rVZlW(o8v1mkVmy!e z4y=2EeR7Cnfw=a_lZm_+SWLf;1mIE%U~oP!VsIY;4`1*ak!Sjn^q^nIcPFNANiKBp zhE9FxbJC5zC+_HjvKM_()}l{JFZ!nJLLZd?HKwnMGWx;Tqwh+S9Mcb`AN^p$(5Iyv zeOnUI$3-1|T}sjCr4oH#s?Z1KHO3CiCt}@x>`OwNO2o}Wo@(S}14jsOO$JW(ZU*-f z@F)eZA3aRp8!7N!4ITWTODA-SMjxGU^wlXupPe}L-O)uK9#Qn=VNGWG^gKr2o?YnU z6N`RB%hBhD+xN!{eSiYc7ic^B1WBN8kR1959sQ@@P%`=r9l-nwtXqP8s)+L#ak+U8 zq2JI(;7A3ogTOg?1JfsIvIN6}2fXHi=L+yPg$@s)%MIuhi9SsS(a)+7eVo`U(a-7< z)8|QeH`C856@8#Y(HBY_{j8$VH!2c+q-@bw>InKwji8^^RrH~{hrU#+(5LDm`c^GN zAFG)duVQ{Z*4@Xx%ZMX_epcMRS;@%j2prr#Tf^vQwHvtiqA!<}8pEr}n&DZvg5ezm z9kxT4>(I$Igy{no;KB3-vnpl!gvmuSeZ!>CM@$KQ#jc~z*dg>w3qT*Tujot0kA7*H z=v$VIK4y97Yjy#B&YIEp?EOFe(q5xaS`5Y-^hTaQL zT{-&Ltw%pTBlNkeMc+Fm^ucpMU%UqN$+Jb@JRS7W<9)~UCi9 zQ3(AW7o%@uGWs}*px@(d^m*it+&+*+SSO8rBZ^Gl$PvV?MxJBHD*_zcK9q96SqCE&K)j*%k)#$q!hCZBI(3ex> zecAu|iC#qC&H)XkpJ)X7dM^E^pJ+b%fa;(xr~>+g9!K9$4fGK`gwX)=@>s|1Lt2G6 ze&{DEk30#;+YTItz}1OyB{~=GCz-3H$62=M&9 z!sv%A1zrolGZVaTLWf}JG7tTbmC;9f)SBtH9gIHHm(h3n=a!8B_1msOU+OsYsjfiZ zY8CXczJb2hq3Cmc3w^IYqYw6Z^u^}($-afY**55-9f7{ucIdNx3u7qeFJYY}_FYB4 zZ9~N6=CMUyU*LEHT*|PAm*``jj=tvH zKIc2o_k0cdpl?E7^r_dGe(qu#Oh0#b^ielQU-if6v(A5y>F4f>KI}2*%dU$)?XBqB zo{T>3F&MvKUJvUMu=RRZIbAa&+FhSq_gXqJ*8-4ln(Wl=Uefviz zFn$3_=$eg9K<1rD+a%S3&>jq92&r73Y?RH zTNFGL!Rr}#x`6jr=&%mDJcVC?7-jq@{NTH=5T8VeiY{LqiFia zk3tbX5bNLzp$eafA^1kb!bgJZE3pqg6XNilm;@h+=@R}JQsuaEa>0^T_mAX4SYO8;p;)*^KlNo9~}6A{Dm*b zGx&t$!S96YBk};gA{p=*ae?p1F8GiH!k45TJ|*(-E!hMglj-m^nFgPeMHuIEeNeD& z9`*I`}ntLZ|icx!MNbD<}AYEkDKhfvts4mK}Vv!r`N} zAAVp%@Yynb$M}H-z=vx+e7U~Br^^SvUF-(N56lR@UP4nCKd?6yjPDnJ72^YTsD>H! zN}2g@lbCf|ux}7?jv{U>@*F|leBj6cu5-Y-8Mt47hdFrhg6DYfHh>O`pvx@i6#szn zn;SRE_{}B3XYC|>*ZA!izd0WGvZcbO&9#d0o4W(Qxz7EJuNzMhn{LbGS%bWKFBlxz z@WC^3X8h)S;WzgYJpRH*uM2*2%HaI~IyA#?&J#Mt!KY8rkMT1+4j;c-`1&n^&tKtF z#?R0tf$=l!(P4Z6Pr%PG5PpW*@DZF@#`qa}z-O=%zJn{^Lud$JLL>MTio>^14?c$b zFfPD6cin01(?Xny@JVz>o{7kN!V$iTz$F8mEL{e7G_Uw=oI4%fB(cj^BM5 zUG$*S+fN#J=EybnW%jnE|(w22)?*N-i%M~H2CI<*f4&-tKh5q2tK=t@ZGI}5ARI)^0vdLci%sLzXt#K{hor~ zuMy^(;rDwO`@SO1WW=rb$M2VmV*_wS!tYnzhVlDl!MC_6hVlE2hTks%@5|7E>w|0z zokZc2oDbjRF8E0wgs-xx1mm;Z3EyQle3*ye%ghg-<{R*Bo&z7}Pw;i-fzPuGe4oeb zF@Dk;;3v)XiPnK{v=@A&Z809h{2us8%VQtcmzwJ*Jr#KxkarqzcmUT@;4B611o&ca zgHLvZVF!Ji{#{0KWCP@Ubt3ue~>X?tS2U zU#Mt(Vo-}i#%z24H*)8Fu#UTrjW`bxR~vbfkyi^iVu6c`Qx&-N!D9h<=~)L`eO+lo zYPA^NJRXb=eb6NmI{8QDUdValK&Q@lQR~mLp*+RMSB;K0qNe=rratct$U)a>dvzz; z5s{Z_JH#Kl6P>bkdtD`;lF;uPlLrqallZIAqQZ5>M9M_rx!KDaA|iIher;0&u}LVF z>Nwg$mMH2Sh?>tK?i(^p2e)#FN9m|`oi>NubqznPa*jiShpu)!x8#r^bDQPbHXKr9 zs{iVR-I!h&dodq~bspH~gg8!!D~mkN$lEQ-;1B?=aNtw~?!Dk42wvvkX(hn$4k%=F z*gF0Hsp;IKq0YlFt}VN&ZA+>4uo#N0E=#d#nckcglSZf^sa!HstX9~<8mU}{BpD=> zOgFX8cGC<|)>0uxZ^X0twa0qICo40fLaxM(l&B+Ot7K6N@TD1A89R8iW zI_J{jhmh235XccIA&`52lB3UKa58@4=ph`0wy8mz_b+6?=tuntCqD+V&yD}y8^ZvX zb+^Ak#6X6dY4_FR3|ws6@q2C=1EsgG`evVJpbXz)oQrX9Bi^N;jv?wQ(8mFNS0RTG zxw?>ZLnGzp+@U#MV=fy$%_-SNbNiUm9@NRSm(Q>#p4uTIjdz&TtzI=TzvmEX-dZjj z){iA?9IXTfe0Q*$DZ3KqA%w*_xwb$h^)H^bOnq{q4R(wce3O<%_Aey@Z2TOtfsFd$IO z-&p0$Krp^aoM+))0^X&f&OOu>p-(6JnrJD933Ba6&Sjrc?o&lH2d9wc(krAnLlSB3 z0u}9H-U8Z76!uhixNw}aJ&$lR%?_UGJVxq1?{Dcd&mu0n&qbOgrjqCy&k=TVJPF}% zZEu+2PbBLnwGAnwK;sZ!CEsunI7KdwqkIBZ!;!Cgy?%vR(Md8Rm*0RL;v%}A*$G{W z0S~{79xy3YKJj&a43#rYVm$rSV45`agO{fUUN;%kzLjg>%S_vcONKPS4ytw3Xf?17 zUt65N!aXy*lcUa`r>JfN`f$j zfjw0`z3rPeRzOl7>}9?%%_jzmfM+AxToSO$zuNLmHc{7ex2YU5$fp*AZILAh$@A3@ zdWzr2kZuPhCn}-@W+bg4#`-c?;ibBsbM6lKPkiMeDQt$aQ||f29NqtYtovfcxUwDO z>wOYBRGlC;O8@a&c{i+7jkxo0E2}3a((_64?UK3w-V_ssiGApi@MFX!%yqH6CyPYq zb-<0}STM3l?sYzV4){ghkB1Le!nXGfH+BhXA(;7^{lT&x%A=e`F-M!Ab;jhJWzSnc zbX~bD)T$M}DYdNB>i*y2*3R$pD@!|Iy?^Phm%*LjIWp_pvvx}0;>*VQ8QgQkyUVE4 shPo2;DMa5}? Date: Sat, 14 Oct 2017 14:23:57 -0500 Subject: [PATCH 6/9] Fix get_array_slice typemaps * Allow user arrays * Allow complex arrays --- python/meep.i | 35 +++++++++++++- python/simulation.py | 33 ++++++++----- python/tests/cavity_arrayslice.py | 77 ++++++++++++++++++++----------- src/array_slice.cpp | 9 ++-- src/meep.hpp | 10 ++-- 5 files changed, 111 insertions(+), 53 deletions(-) diff --git a/python/meep.i b/python/meep.i index 3abf1630b..a7dbbb814 100644 --- a/python/meep.i +++ b/python/meep.i @@ -332,8 +332,39 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, } // Typemap suite for array_slice -// TODO: add (cdouble *, int) version -%apply (double* INPLACE_ARRAY1, int DIM1) {(double *slice, int slice_length)}; + +%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* slice { + $1 = is_array($input); +} + +%typemap(in, fragment="NumPy_Macros") double* slice { + $1 = (double *)array_data($input); +} + +%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex* slice { + $1 = is_array($input); +} + +%typemap(in) std::complex* slice { + $1 = (std::complex *)array_data($input); +} + +%typecheck(SWIG_TYPECHECK_POINTER) meep::component { + $1 = PyLong_Check($input) && PyLong_AsLong($input) < 100; +} + +%typemap(in) meep::component { + $1 = static_cast(PyLong_AsLong($input)); +} + +%typecheck(SWIG_TYPECHECK_POINTER) meep::derived_component { + $1 = PyLong_Check($input) && PyLong_AsLong($input) >= 100; +} + +%typemap(in) meep::derived_component { + $1 = static_cast(PyLong_AsLong($input)); +} + %apply int INPLACE_ARRAY1[ANY] { int [3] }; // Rename python builtins diff --git a/python/simulation.py b/python/simulation.py index 4a77df8e5..abfef2d59 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -652,23 +652,34 @@ def output_components(self, fname, *components): if self.output_append_h5 is None: self.output_h5_hook(self.fields.h5file_name(fname, self._get_filename_prefix(), True)) - def get_array(self, center, size, component=mp.Ez, arr=None): - self.dim_sizes = np.zeros(3, dtype=np.int32) + def get_array(self, center, size, component=mp.Ez, cmplx=False, arr=None): + dim_sizes = np.zeros(3, dtype=np.int32) + # TODO(chogan): Account for cylindrical Volume after PR 106 is merged vol = Volume(center, size=size, dims=self.dimensions) - rank = self.fields.get_array_slice_dimensions(vol.swigobj, self.dim_sizes) - sizes = self.dim_sizes.copy() - sizes[rank:] = 1 - n = np.prod(sizes) + self.fields.get_array_slice_dimensions(vol.swigobj, dim_sizes) + + dims = [s for s in dim_sizes if s != 0] if arr is not None: - # TODO(chogan): Check size - pass + if cmplx and not np.iscomplexobj(arr): + raise ValueError("Requested a complex slice, but provided array of type {}.".format(arr.dtype)) + + for a, b in zip(arr.shape, dims): + if a != b: + fmt = "Expected dimensions {}, but got {}" + raise ValueError(fmt.format(dims, arr.shape)) + + arr = np.require(arr, requirements=['C', 'W']) + else: - arr = np.zeros(n, dtype=np.float64) + arr = np.zeros(dims, dtype=np.complex128 if cmplx else np.float64) - self.fields.get_array_slice(vol.swigobj, component, arr) + if np.iscomplexobj(arr): + self.fields.get_complex_array_slice(vol.swigobj, component, arr) + else: + self.fields.get_array_slice(vol.swigobj, component, arr) - return arr.reshape([s for s in self.dim_sizes if s != 0]) + return arr def change_k_point(self, k): self.k_point = k diff --git a/python/tests/cavity_arrayslice.py b/python/tests/cavity_arrayslice.py index a4e399f3b..93dfd09c0 100644 --- a/python/tests/cavity_arrayslice.py +++ b/python/tests/cavity_arrayslice.py @@ -10,6 +10,8 @@ class TestCavityArraySlice(unittest.TestCase): data_dir = os.path.abspath(os.path.realpath(os.path.join(__file__, '..', 'data'))) + expected_1d = np.load(os.path.join(data_dir, 'cavity_arrayslice_1d.npy')) + expected_2d = np.load(os.path.join(data_dir, 'cavity_arrayslice_2d.npy')) def setUp(self): @@ -48,46 +50,65 @@ def setUp(self): self.y_min = -0.15 * sy self.y_max = +0.15 * sy + self.size_1d = mp.Vector3(self.x_max - self.x_min) + self.center_1d = mp.Vector3((self.x_min + self.x_max) / 2) + + self.size_2d = mp.Vector3(self.x_max - self.x_min, self.y_max - self.y_min) + self.center_2d = mp.Vector3((self.x_min + self.x_max) / 2, (self.y_min + self.y_max) / 2) + def test_1d_slice(self): self.sim.run(until_after_sources=0) + hl_slice1d = self.sim.get_array(center=self.center_1d, size=self.size_1d, component=mp.Hz) + np.testing.assert_allclose(self.expected_1d, hl_slice1d) - # Low level 1D slice of Hz data - # dims1d = np.zeros(3, dtype=np.int32) - # v1d = mp.volume(mp.vec(self.x_min, 0.0), mp.vec(self.x_max, 0.0)) - # self.sim.fields.get_array_slice_dimensions(v1d, dims1d) - # NX1 = dims1d[0] - # slice1d = np.zeros(NX1, dtype=np.double) - # self.sim.fields.get_array_slice(v1d, mp.Hz, slice1d) - - expected = np.load(os.path.join(self.data_dir, 'cavity_arrayslice_1d.npy')) + def test_2d_slice(self): + self.sim.run(until_after_sources=0) + hl_slice2d = self.sim.get_array(center=self.center_2d, size=self.size_2d, component=mp.Hz) + np.testing.assert_allclose(self.expected_2d, hl_slice2d) - size = mp.Vector3(self.x_max - self.x_min) - center = mp.Vector3((self.x_min + self.x_max) / 2) - hl_slice1d = self.sim.get_array(center=center, size=size, component=mp.Hz) + def test_1d_slice_user_array(self): + self.sim.run(until_after_sources=0) + arr = np.zeros(126, dtype=np.float64) + self.sim.get_array(center=self.center_1d, size=self.size_1d, component=mp.Hz, arr=arr) + np.testing.assert_allclose(self.expected_1d, arr) - np.testing.assert_allclose(expected, hl_slice1d) + def test_2d_slice_user_array(self): + self.sim.run(until_after_sources=0) + arr = np.zeros((126, 38), dtype=np.float64) + self.sim.get_array(center=self.center_2d, size=self.size_2d, component=mp.Hz, arr=arr) + np.testing.assert_allclose(self.expected_2d, arr) - def test_2d_slice(self): + def test_illegal_user_array(self): self.sim.run(until_after_sources=0) - # Low level 2D slice of Hz data - # dims2d = np.zeros(3, dtype=np.int32) - # v2d = mp.volume(mp.vec(self.x_min, self.y_min), mp.vec(self.x_max, self.y_max)) - # self.sim.fields.get_array_slice_dimensions(v2d, dims2d) - # NX2 = dims2d[0] - # NY2 = dims2d[1] - # N2 = NX2 * NY2 - # slice2d = np.zeros(N2, dtype=np.double) - # self.sim.fields.get_array_slice(v2d, mp.Hz, slice2d) + with self.assertRaises(ValueError): + arr = np.zeros(128) + self.sim.get_array(center=self.center_1d, size=self.size_1d, + component=mp.Hz, arr=arr) - expected = np.load(os.path.join(self.data_dir, 'cavity_arrayslice_2d.npy')) + with self.assertRaises(ValueError): + arr = np.zeros((126, 39)) + self.sim.get_array(center=self.center_2d, size=self.size_2d, + component=mp.Hz, arr=arr) - size = mp.Vector3(self.x_max - self.x_min, self.y_max - self.y_min) - center = mp.Vector3((self.x_min + self.x_max) / 2, (self.y_min + self.y_max) / 2) - hl_slice2d = self.sim.get_array(center=center, size=size, component=mp.Hz) + with self.assertRaises(ValueError): + arr = np.zeros((126, 38)) + self.sim.get_array(center=self.center_2d, size=self.size_2d, + component=mp.Hz, cmplx=True, arr=arr) - np.testing.assert_allclose(expected, hl_slice2d) + def test_1d_complex_slice(self): + self.sim.run(until_after_sources=0) + hl_slice1d = self.sim.get_array(center=self.center_1d, size=self.size_1d, + component=mp.Hz, cmplx=True) + self.assertTrue(hl_slice1d.dtype == np.complex128) + self.assertTrue(hl_slice1d.shape[0] == 126) + def test_2d_complex_slice(self): + self.sim.run(until_after_sources=0) + hl_slice2d = self.sim.get_array(center=self.center_2d, size=self.size_2d, + component=mp.Hz, cmplx=True) + self.assertTrue(hl_slice2d.dtype == np.complex128) + self.assertTrue(hl_slice2d.shape[0] == 126 and hl_slice2d.shape[1] == 38) if __name__ == '__main__': unittest.main() diff --git a/src/array_slice.cpp b/src/array_slice.cpp index e6d9d452b..934ea15c2 100644 --- a/src/array_slice.cpp +++ b/src/array_slice.cpp @@ -461,9 +461,8 @@ cdouble *fields::get_complex_array_slice(const volume &where, } double *fields::get_array_slice(const volume &where, component c, - double *slice, int slice_length) + double *slice) { - (void) slice_length; std::vector components(1); components[0]=c; return (double *)do_get_array_slice(where, components, @@ -473,9 +472,8 @@ double *fields::get_array_slice(const volume &where, component c, double *fields::get_array_slice(const volume &where, derived_component c, - double *slice, int slice_length) + double *slice) { - (void) slice_length; int nfields; component carray[12]; field_rfunction rfun = derived_component_func(c, gv, nfields, carray); @@ -486,9 +484,8 @@ double *fields::get_array_slice(const volume &where, } cdouble *fields::get_complex_array_slice(const volume &where, component c, - cdouble *slice, int slice_length) + cdouble *slice) { - (void) slice_length; std::vector components(1); components[0]=c; return (cdouble *)do_get_array_slice(where, components, diff --git a/src/meep.hpp b/src/meep.hpp index 782a8c589..3b75c95a3 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1314,14 +1314,12 @@ class fields { // alternative entry points for when you have no field // function, i.e. you want just a single component or - // derived component.) (The slice_length parameter is only - // present to facilitate SWIG-generated python wrappers; - // it is ignored by the C code). - double *get_array_slice(const volume &where, component c, double *slice=0, int slice_length=0); - double *get_array_slice(const volume &where, derived_component c, double *slice=0, int slice_length=0); + // derived component.) + double *get_array_slice(const volume &where, component c, double *slice=0); + double *get_array_slice(const volume &where, derived_component c, double *slice=0); std::complex *get_complex_array_slice(const volume &where, component c, - std::complex *slice=0, int slice_length=0); + std::complex *slice=0); // master routine for all above entry points void *do_get_array_slice(const volume &where, From c5181becea451fe12240e354515bc3790f3e1f13 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sat, 14 Oct 2017 15:27:57 -0500 Subject: [PATCH 7/9] Py2 compatability --- python/meep.i | 8 ++++---- python/typemap_utils.cpp | 4 ++++ 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/python/meep.i b/python/meep.i index a7dbbb814..9dc32ebda 100644 --- a/python/meep.i +++ b/python/meep.i @@ -350,19 +350,19 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, } %typecheck(SWIG_TYPECHECK_POINTER) meep::component { - $1 = PyLong_Check($input) && PyLong_AsLong($input) < 100; + $1 = PyInteger_Check($input) && PyInteger_AsLong($input) < 100; } %typemap(in) meep::component { - $1 = static_cast(PyLong_AsLong($input)); + $1 = static_cast(PyInteger_AsLong($input)); } %typecheck(SWIG_TYPECHECK_POINTER) meep::derived_component { - $1 = PyLong_Check($input) && PyLong_AsLong($input) >= 100; + $1 = PyInteger_Check($input) && PyInteger_AsLong($input) >= 100; } %typemap(in) meep::derived_component { - $1 = static_cast(PyLong_AsLong($input)); + $1 = static_cast(PyInteger_AsLong($input)); } %apply int INPLACE_ARRAY1[ANY] { int [3] }; diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index d65b9f7e8..4184a112b 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -19,8 +19,12 @@ #if PY_MAJOR_VERSION >= 3 #define PyObject_ToCharPtr(n) PyUnicode_AsUTF8(n) + #define PyInteger_Check(n) PyLong_Check(n) + #define PyInteger_AsLong(n) PyLong_AsLong(n) #else #define PyObject_ToCharPtr(n) PyString_AsString(n) + #define PyInteger_Check(n) PyInt_Check(n) + #define PyInteger_AsLong(n) PyInt_AsLong(n) #endif static PyObject *py_geometric_object() { From b18022d2ef0fcb4571e1013404762aa07b3035ff Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 19 Oct 2017 20:18:25 -0500 Subject: [PATCH 8/9] Default complex in get_array --- python/simulation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python/simulation.py b/python/simulation.py index abfef2d59..8c81bf565 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -652,7 +652,7 @@ def output_components(self, fname, *components): if self.output_append_h5 is None: self.output_h5_hook(self.fields.h5file_name(fname, self._get_filename_prefix(), True)) - def get_array(self, center, size, component=mp.Ez, cmplx=False, arr=None): + def get_array(self, center, size, component=mp.Ez, cmplx=None, arr=None): dim_sizes = np.zeros(3, dtype=np.int32) # TODO(chogan): Account for cylindrical Volume after PR 106 is merged vol = Volume(center, size=size, dims=self.dimensions) @@ -660,6 +660,9 @@ def get_array(self, center, size, component=mp.Ez, cmplx=False, arr=None): dims = [s for s in dim_sizes if s != 0] + if cmplx is None: + cmplx = component < mp.Dielectric and not self.fields.is_real + if arr is not None: if cmplx and not np.iscomplexobj(arr): raise ValueError("Requested a complex slice, but provided array of type {}.".format(arr.dtype)) From 3798ab07ff68223c61ed376c0d444713c0647071 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 26 Oct 2017 18:20:34 -0500 Subject: [PATCH 9/9] Account for cylindrical Volume --- python/simulation.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 8c81bf565..c3f82bfe6 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -654,8 +654,7 @@ def output_components(self, fname, *components): def get_array(self, center, size, component=mp.Ez, cmplx=None, arr=None): dim_sizes = np.zeros(3, dtype=np.int32) - # TODO(chogan): Account for cylindrical Volume after PR 106 is merged - vol = Volume(center, size=size, dims=self.dimensions) + vol = Volume(center, size=size, dims=self.dimensions, is_cylindrical=self.is_cylindrical) self.fields.get_array_slice_dimensions(vol.swigobj, dim_sizes) dims = [s for s in dim_sizes if s != 0]