From 66bdb9c4b65138e13b6ba45e998a2a7eca6f2a3b Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 28 Nov 2017 11:44:37 -0800 Subject: [PATCH 01/15] added B factor reference implementation --- test/scatter/test_scatter.py | 120 +++++++++++++++++++++++++++++++++-- 1 file changed, 114 insertions(+), 6 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 3f7bc47..03149fb 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -167,7 +167,7 @@ def form_factor_reference(qvector, atomz): def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, - dont_rotate=False): + dont_rotate=False, Bfactors=None): """ Simulate a single x-ray scattering shot off an ensemble of identical molecules. @@ -191,6 +191,10 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, rfloats : ndarray, float, n x 3 A bunch of random floats, uniform on [0,1] to be used to seed the quaternion computation. + + B_factors : ndarray, float, 1d or 2d + An n x 1 or n x 3 x 3 array of B-factors for isotropic and anisotropic + cases, respectively. Returns ------- @@ -206,6 +210,12 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, else: rs = np.random.RandomState(RANDOM_SEED) rfs = rs.rand(3, num_molecules) + + # convert B factors in PDB file to V matrix + if Bfactors is not None: + V = Bfactors/(8*np.square(np.pi)) + else: + V = np.zeros(xyzlist.shape[0]) for i,qvector in enumerate(q_grid): for n in range(num_molecules): @@ -218,8 +228,14 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, q_mag = np.linalg.norm(qvector) fi = scatter.atomic_formfactor(atomic_numbers[j], q_mag) r = rotated_xyzlist[j,:] - A[i] += fi * np.sin( np.dot(qvector, r) ) - A[i] += 1j * fi * np.cos( np.dot(qvector, r) ) + + if len(V.shape) == 1: + qVq = np.square(q_mag)*V[j] + else: + qVq = np.dot(np.dot(qvector, V[j]), qvector) + + A[i] += fi * np.sin( np.dot(qvector, r) ) * np.exp(- 0.5 * qVq) + A[i] += 1j * fi * np.cos( np.dot(qvector, r) ) * np.exp(- 0.5 * qVq) return A @@ -347,9 +363,20 @@ def setup(self): self.q_grid = np.loadtxt(ref_file('512_q.xyz'))[:self.nq] - self.num_molecules = 512 + self.num_molecules = 512 + self.iso_B = 10.0*np.random.rand(self.nr) + self.aniso_B = 10.0*np.abs(np.random.randn(self.nr, 3, 3)) + for i in range(self.nr): + self.aniso_B[i] += self.aniso_B[i].T + self.ref_A = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid) + self.ref_A_isoB = ref_simulate_shot(self.xyzlist, self.atomic_numbers, + self.num_molecules, self.q_grid, + Bfactors=self.iso_B) + self.ref_A_anisoB = ref_simulate_shot(self.xyzlist, self.atomic_numbers, + self.num_molecules, self.q_grid, + Bfactors=self.aniso_B) def test_cpu_scatter(self): @@ -371,6 +398,44 @@ def test_cpu_scatter(self): assert not np.all( cpu_A == 0.0 ) assert not np.sum( cpu_A == np.nan ) + def test_cpu_scatter_isoB(self): + + cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + + print 'num_molecules:', self.num_molecules + cpu_A_isoB = _cppscatter.cpp_scatter(self.num_molecules, + self.xyzlist, + self.q_grid, + atom_types, + cromermann_parameters, + device_id='CPU', + random_state=np.random.RandomState(RANDOM_SEED), + Bfactors=self.iso_B) + + assert_allclose(cpu_A_isoB, self.ref_A_isoB, rtol=1e-3, atol=1.0, + err_msg='scatter: c-cpu/cpu reference mismatch with isotropic B factors') + assert not np.all( cpu_A_isoB == 0.0 ) + assert not np.sum( cpu_A_isoB == np.nan ) + + def test_cpu_scatter_anisoB(self): + + cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + + print 'num_molecules:', self.num_molecules + cpu_A_anisoB = _cppscatter.cpp_scatter(self.num_molecules, + self.xyzlist, + self.q_grid, + atom_types, + cromermann_parameters, + device_id='CPU', + random_state=np.random.RandomState(RANDOM_SEED), + Bfactors=self.aniso_B) + + assert_allclose(cpu_A_anisoB, self.ref_A_anisoB, rtol=1e-3, atol=1.0, + err_msg='scatter: c-cpu/cpu reference mismatch with anisotropic B factors') + assert not np.all( cpu_A_anisoB == 0.0 ) + assert not np.sum( cpu_A_anisoB == np.nan ) + def test_gpu_scatter(self): if not GPU: raise SkipTest @@ -389,8 +454,51 @@ def test_gpu_scatter(self): assert_allclose(gpu_A, self.ref_A, rtol=1e-3, atol=1.0, err_msg='scatter: cuda-gpu/cpu reference mismatch') assert not np.all( gpu_A == 0.0 ) - assert not np.sum( gpu_A == np.nan ) - + assert not np.sum( gpu_A == np.nan ) + + def test_gpu_scatter_isoB(self): + + if not GPU: raise SkipTest + + cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + + print 'num_molecules:', self.num_molecules + gpu_A_isoB = _cppscatter.cpp_scatter(self.num_molecules, + self.xyzlist, + self.q_grid, + atom_types, + cromermann_parameters, + device_id=0, + random_state=np.random.RandomState(RANDOM_SEED), + Bfactors=self.iso_B) + + assert_allclose(gpu_A_isoB, self.ref_A_isoB, rtol=1e-3, atol=1.0, + err_msg='scatter: cuda-gpu/cpu reference mismatch with isotropic B factors') + assert not np.all( gpu_A_isoB == 0.0 ) + assert not np.sum( gpu_A_isoB == np.nan ) + + def test_gpu_scatter_anisoB(self): + + if not GPU: raise SkipTest + + cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + + print 'num_molecules:', self.num_molecules + gpu_A_anisoB = _cppscatter.cpp_scatter(self.num_molecules, + self.xyzlist, + self.q_grid, + atom_types, + cromermann_parameters, + device_id=0, + random_state=np.random.RandomState(RANDOM_SEED), + Bfactors=self.aniso_B) + + assert_allclose(gpu_A_anisoB, self.ref_A_anisoB, rtol=1e-3, atol=1.0, + err_msg='scatter: cuda-gpu/cpu reference mismatch with anisotropic B factors') + assert not np.all( gpu_A_anisoB == 0.0 ) + assert not np.sum( gpu_A_anisoB == np.nan ) + + def test_parallel_interface(self): cp, at = get_cromermann_parameters(self.atomic_numbers) From b416342e92c42399ebc29258bb66db7af2b4c8d0 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 28 Nov 2017 15:09:10 -0800 Subject: [PATCH 02/15] added reference implementation for B factors --- test/scatter/test_scatter.py | 88 ++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 03149fb..448872c 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -167,7 +167,7 @@ def form_factor_reference(qvector, atomz): def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, - dont_rotate=False, Bfactors=None): + dont_rotate=False, U=None): """ Simulate a single x-ray scattering shot off an ensemble of identical molecules. @@ -192,9 +192,11 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, A bunch of random floats, uniform on [0,1] to be used to seed the quaternion computation. - B_factors : ndarray, float, 1d or 2d - An n x 1 or n x 3 x 3 array of B-factors for isotropic and anisotropic - cases, respectively. + U : ndarray, float, 1d or 2d + An n x 1 or n x 3 x 3 array of atomic displacement parameters for + isotropic and anisotropic cases, respectively. Related to B factors + listed in PDBs by factors of 8*pi^2 in the isotropic case, and 10^4 + in the anisotropic case. Returns ------- @@ -212,10 +214,8 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, rfs = rs.rand(3, num_molecules) # convert B factors in PDB file to V matrix - if Bfactors is not None: - V = Bfactors/(8*np.square(np.pi)) - else: - V = np.zeros(xyzlist.shape[0]) + if U is None: + U = np.zeros(xyzlist.shape[0]) for i,qvector in enumerate(q_grid): for n in range(num_molecules): @@ -229,13 +229,13 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, fi = scatter.atomic_formfactor(atomic_numbers[j], q_mag) r = rotated_xyzlist[j,:] - if len(V.shape) == 1: - qVq = np.square(q_mag)*V[j] + if len(U.shape) == 1: + qUq = np.square(q_mag)*U[j] else: - qVq = np.dot(np.dot(qvector, V[j]), qvector) + qUq = np.dot(np.dot(qvector, U[j]), qvector) - A[i] += fi * np.sin( np.dot(qvector, r) ) * np.exp(- 0.5 * qVq) - A[i] += 1j * fi * np.cos( np.dot(qvector, r) ) * np.exp(- 0.5 * qVq) + A[i] += fi * np.sin( np.dot(qvector, r) ) * np.exp(- 0.5 * qUq) + A[i] += 1j * fi * np.cos( np.dot(qvector, r) ) * np.exp(- 0.5 * qUq) return A @@ -364,19 +364,19 @@ def setup(self): self.q_grid = np.loadtxt(ref_file('512_q.xyz'))[:self.nq] self.num_molecules = 512 - self.iso_B = 10.0*np.random.rand(self.nr) - self.aniso_B = 10.0*np.abs(np.random.randn(self.nr, 3, 3)) + self.iso_U = np.random.rand(self.nr) + self.aniso_U = np.abs(np.random.randn(self.nr, 3, 3)) for i in range(self.nr): - self.aniso_B[i] += self.aniso_B[i].T + self.aniso_U[i] += self.aniso_U[i].T self.ref_A = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid) - self.ref_A_isoB = ref_simulate_shot(self.xyzlist, self.atomic_numbers, + self.ref_A_isoU = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid, - Bfactors=self.iso_B) - self.ref_A_anisoB = ref_simulate_shot(self.xyzlist, self.atomic_numbers, + U = self.iso_U) + self.ref_A_anisoU = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid, - Bfactors=self.aniso_B) + U = self.aniso_U) def test_cpu_scatter(self): @@ -398,43 +398,43 @@ def test_cpu_scatter(self): assert not np.all( cpu_A == 0.0 ) assert not np.sum( cpu_A == np.nan ) - def test_cpu_scatter_isoB(self): + def test_cpu_scatter_isoU(self): cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) print 'num_molecules:', self.num_molecules - cpu_A_isoB = _cppscatter.cpp_scatter(self.num_molecules, + cpu_A_isoU = _cppscatter.cpp_scatter(self.num_molecules, self.xyzlist, self.q_grid, atom_types, cromermann_parameters, device_id='CPU', random_state=np.random.RandomState(RANDOM_SEED), - Bfactors=self.iso_B) + U=self.iso_U) - assert_allclose(cpu_A_isoB, self.ref_A_isoB, rtol=1e-3, atol=1.0, + assert_allclose(cpu_A_isoU, self.ref_A_isoU, rtol=1e-3, atol=1.0, err_msg='scatter: c-cpu/cpu reference mismatch with isotropic B factors') - assert not np.all( cpu_A_isoB == 0.0 ) - assert not np.sum( cpu_A_isoB == np.nan ) + assert not np.all( cpu_A_isoU == 0.0 ) + assert not np.sum( cpu_A_isoU == np.nan ) - def test_cpu_scatter_anisoB(self): + def test_cpu_scatter_anisoU(self): cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) print 'num_molecules:', self.num_molecules - cpu_A_anisoB = _cppscatter.cpp_scatter(self.num_molecules, + cpu_A_anisoU = _cppscatter.cpp_scatter(self.num_molecules, self.xyzlist, self.q_grid, atom_types, cromermann_parameters, device_id='CPU', random_state=np.random.RandomState(RANDOM_SEED), - Bfactors=self.aniso_B) + U=self.aniso_U) - assert_allclose(cpu_A_anisoB, self.ref_A_anisoB, rtol=1e-3, atol=1.0, + assert_allclose(cpu_A_anisoU, self.ref_A_anisoU, rtol=1e-3, atol=1.0, err_msg='scatter: c-cpu/cpu reference mismatch with anisotropic B factors') - assert not np.all( cpu_A_anisoB == 0.0 ) - assert not np.sum( cpu_A_anisoB == np.nan ) + assert not np.all( cpu_A_anisoU == 0.0 ) + assert not np.sum( cpu_A_anisoU == np.nan ) def test_gpu_scatter(self): @@ -456,47 +456,47 @@ def test_gpu_scatter(self): assert not np.all( gpu_A == 0.0 ) assert not np.sum( gpu_A == np.nan ) - def test_gpu_scatter_isoB(self): + def test_gpu_scatter_isoU(self): if not GPU: raise SkipTest cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) print 'num_molecules:', self.num_molecules - gpu_A_isoB = _cppscatter.cpp_scatter(self.num_molecules, + gpu_A_isoU = _cppscatter.cpp_scatter(self.num_molecules, self.xyzlist, self.q_grid, atom_types, cromermann_parameters, device_id=0, random_state=np.random.RandomState(RANDOM_SEED), - Bfactors=self.iso_B) + U=self.iso_U) - assert_allclose(gpu_A_isoB, self.ref_A_isoB, rtol=1e-3, atol=1.0, + assert_allclose(gpu_A_isoU, self.ref_A_isoU, rtol=1e-3, atol=1.0, err_msg='scatter: cuda-gpu/cpu reference mismatch with isotropic B factors') - assert not np.all( gpu_A_isoB == 0.0 ) - assert not np.sum( gpu_A_isoB == np.nan ) + assert not np.all( gpu_A_isoU == 0.0 ) + assert not np.sum( gpu_A_isoU == np.nan ) - def test_gpu_scatter_anisoB(self): + def test_gpu_scatter_anisoU(self): if not GPU: raise SkipTest cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) print 'num_molecules:', self.num_molecules - gpu_A_anisoB = _cppscatter.cpp_scatter(self.num_molecules, + gpu_A_anisoU = _cppscatter.cpp_scatter(self.num_molecules, self.xyzlist, self.q_grid, atom_types, cromermann_parameters, device_id=0, random_state=np.random.RandomState(RANDOM_SEED), - Bfactors=self.aniso_B) + U=self.aniso_U) - assert_allclose(gpu_A_anisoB, self.ref_A_anisoB, rtol=1e-3, atol=1.0, + assert_allclose(gpu_A_anisoU, self.ref_A_anisoU, rtol=1e-3, atol=1.0, err_msg='scatter: cuda-gpu/cpu reference mismatch with anisotropic B factors') - assert not np.all( gpu_A_anisoB == 0.0 ) - assert not np.sum( gpu_A_anisoB == np.nan ) + assert not np.all( gpu_A_anisoU == 0.0 ) + assert not np.sum( gpu_A_anisoU == np.nan ) def test_parallel_interface(self): From 25aadcec37856c3e274b199598a596bc111bd2b1 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 28 Nov 2017 15:18:50 -0800 Subject: [PATCH 03/15] added python reference implementation for b-factors --- test/scatter/test_scatter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 448872c..a219168 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -364,10 +364,11 @@ def setup(self): self.q_grid = np.loadtxt(ref_file('512_q.xyz'))[:self.nq] self.num_molecules = 512 - self.iso_U = np.random.rand(self.nr) + self.iso_U = np.random.rand(self.nr) / 8.0 self.aniso_U = np.abs(np.random.randn(self.nr, 3, 3)) for i in range(self.nr): self.aniso_U[i] += self.aniso_U[i].T + self.aniso_U /= 30.0 self.ref_A = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid) From 54ebdcb90549db95684120e5f840bc348da8e58d Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 28 Nov 2017 16:20:29 -0800 Subject: [PATCH 04/15] added python reference implementation for b-factors --- test/scatter/test_scatter.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index a219168..b1ea86c 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -364,20 +364,20 @@ def setup(self): self.q_grid = np.loadtxt(ref_file('512_q.xyz'))[:self.nq] self.num_molecules = 512 - self.iso_U = np.random.rand(self.nr) / 8.0 - self.aniso_U = np.abs(np.random.randn(self.nr, 3, 3)) + self.iso_U = np.random.rand(self.nr) / 10.0 + self.aniso_U = np.random.rand(self.nr, 3, 3) for i in range(self.nr): self.aniso_U[i] += self.aniso_U[i].T - self.aniso_U /= 30.0 + self.aniso_U /= 20.0 self.ref_A = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid) self.ref_A_isoU = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid, - U = self.iso_U) + U=self.iso_U) self.ref_A_anisoU = ref_simulate_shot(self.xyzlist, self.atomic_numbers, self.num_molecules, self.q_grid, - U = self.aniso_U) + U=self.aniso_U) def test_cpu_scatter(self): From 4a59e7aa30939cd28f102a2becd59c074cf6b222 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Wed, 29 Nov 2017 11:25:24 -0800 Subject: [PATCH 05/15] added python reference implementation for b-factors --- test/scatter/test_scatter.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index b1ea86c..1d378cd 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -195,8 +195,8 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, U : ndarray, float, 1d or 2d An n x 1 or n x 3 x 3 array of atomic displacement parameters for isotropic and anisotropic cases, respectively. Related to B factors - listed in PDBs by factors of 8*pi^2 in the isotropic case, and 10^4 - in the anisotropic case. + listed in PDB files by factors of 8*pi^2 in the isotropic case, and + 10^4 in the anisotropic case. Returns ------- @@ -213,7 +213,6 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, rs = np.random.RandomState(RANDOM_SEED) rfs = rs.rand(3, num_molecules) - # convert B factors in PDB file to V matrix if U is None: U = np.zeros(xyzlist.shape[0]) @@ -229,6 +228,7 @@ def ref_simulate_shot(xyzlist, atomic_numbers, num_molecules, q_grid, fi = scatter.atomic_formfactor(atomic_numbers[j], q_mag) r = rotated_xyzlist[j,:] + # compute the Debye-Waller factor if len(U.shape) == 1: qUq = np.square(q_mag)*U[j] else: From 6b57f4c0eb0281569f63ced63db2830a52aa2365 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Thu, 30 Nov 2017 14:41:29 -0800 Subject: [PATCH 06/15] B-factors implemented and functioning for CPU --- src/python/scatter.py | 27 ++++++++++++++++-- src/scatter/cpp_scatter.cpp | 49 ++++++++++++++++++++++++++++---- src/scatter/cpp_scatter.cu | 23 ++++++++++----- src/scatter/cpp_scatter.hh | 7 +++++ src/scatter/cpp_scatter_wrap.pyx | 31 +++++++++++++++----- test/scatter/test_scatter.py | 1 + 6 files changed, 116 insertions(+), 22 deletions(-) diff --git a/src/python/scatter.py b/src/python/scatter.py index def7b60..ce9001a 100644 --- a/src/python/scatter.py +++ b/src/python/scatter.py @@ -69,7 +69,7 @@ def rand(self, *args): def simulate_atomic(traj, num_molecules, detector, traj_weights=None, finite_photon=False, ignore_hydrogens=False, dont_rotate=False, procs_per_node=1, - nodes=[], devices=[], random_state=None): + nodes=[], devices=[], random_state=None, U=None): """ Simulate a scattering 'shot', i.e. one exposure of x-rays to a sample. @@ -117,6 +117,10 @@ def simulate_atomic(traj, num_molecules, detector, traj_weights=None, dont_rotate : bool Don't apply a random rotation to the molecule before computing the scattering. Good for debugging and the like. + + U : ndarray, float + Atomic displacement parameters, either a 1d n_atoms array [isotropic case] + or 3d (n_atoms, 3, 3,) tensor. Returns ------- @@ -180,6 +184,21 @@ def simulate_atomic(traj, num_molecules, detector, traj_weights=None, atoms_to_keep = np.ones(n_atoms, dtype=np.bool) + # if a 1d U matrix is passed, expand to 3d; if no U matrix is input, + # set all elements to zero + n_atoms = len(atoms_to_keep) + if U is None: + U = np.zeros((n_atoms, 3, 3)) + elif U.shape == (n_atoms,): + e = np.eye(3) + U = U[:,None,None] * e[None,None,:] + elif U.shape == (n_atoms, 3, 3): + # correctly sized for lower level interface + pass + else: + raise TypeError("U matrix has invalide shape.") + + qxyz = _qxyz_from_detector(detector) amplitudes = np.zeros(qxyz.shape[0], dtype=np.complex128) @@ -196,7 +215,8 @@ def simulate_atomic(traj, num_molecules, detector, traj_weights=None, procs_per_node=procs_per_node, nodes=nodes, devices=devices, - random_state=random_state) + random_state=random_state, + U=U) if finite_photon is not False: intensities = np.square(np.abs(amplitudes)) @@ -254,7 +274,8 @@ def simulate_density(grid, grid_spacing, num_molecules, detector, procs_per_node=procs_per_node, nodes=nodes, devices=devices, - random_state=random_state) + random_state=random_state, + U=None) # put the output back in sq grid form if requested if reshape_output: diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index 6667bab..f2ad636 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -142,6 +142,35 @@ void qVq_product(float const * const V, } +#ifdef __CUDACC__ + __host__ __device__ +#endif +void qUq_product(float const * const U, + int i, + float qx, + float qy, + float qz, + + float &o_qUq + + ){ + + // compute < q | U_ii | q > + + float qUq; + + qUq = qx * qx * U[9*i + 0]; + qUq += qy * qy * U[9*i + 4]; + qUq += qz * qz * U[9*i + 8]; + qUq += 2 * qx * qy * U[9*i + 1]; + qUq += 2 * qx * qz * U[9*i + 2]; + qUq += 2 * qy * qz * U[9*i + 5]; + + o_qUq = qUq; + +} + + /****************************************************************************** * Hybrid CPU/GPU Code * decides whether to try and call GPU code or raise an exception, depending @@ -168,6 +197,9 @@ void gpuscatter (int device_id_, int * h_atom_types, float * h_cromermann, + // atomic displacement parameters + float * h_U, + // random numbers for rotations int n_rotations, float * rand1, @@ -183,7 +215,7 @@ void gpuscatter (int device_id_, _gpuscatter( device_id_, n_q, h_qx, h_qy, h_qz, n_atoms, h_rx, h_ry, h_rz, - n_atom_types, h_atom_types, h_cromermann, + n_atom_types, h_atom_types, h_cromermann, h_U, n_rotations, rand1, rand2, rand3, h_q_out_real, h_q_out_imag); #else @@ -252,6 +284,8 @@ void cpu_kernel( int const n_q, int const * const __restrict__ atom_types, float const * const __restrict__ cromermann, + float const * const __restrict__ U, + int const n_rotations, float const * const __restrict__ randN1, float const * const __restrict__ randN2, @@ -273,6 +307,7 @@ void cpu_kernel( int const n_q, * n_atom_types : the number of unique atom types (formfactors) * atom_types : the atom "type", which is an arbitrary index * cromermann : 9 params specifying formfactor for each atom type + * U : 3d array of the atomic displacement parameters * n_rotations/ : The number of molecules to independently rotate and * randN{1,2,3} the random numbers used to perform those rotations * @@ -289,6 +324,7 @@ void cpu_kernel( int const n_q, float mq, qo, fi; // mag of q, formfactor for atom i float q_sum_real, q_sum_imag; // partial sum of real and imaginary amplitude float qr; // dot product of q and r + float qUq; // Debye Waller factor, qT * U_ii * q // we will use a small array to store form factors float * formfactors = (float *) malloc(n_atom_types * sizeof(float)); @@ -352,10 +388,11 @@ void cpu_kernel( int const n_q, ax, ay, az); qr = ax*qx + ay*qy + az*qz; - + qUq_product(U, a, qx, qy, qz, qUq); + // FIXME :: swap cos and sin???? e^i*t = cos(t) + i sin(t) - q_sum_real += fi * sinf(qr); - q_sum_imag += fi * cosf(qr); + q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq); + q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq); } // finished one atom (3rd loop) } // finished one molecule (2nd loop) @@ -388,6 +425,8 @@ void cpuscatter( int n_q, int * atom_types, float * cromermann, + float * U, + int n_rotations, float * randN1, float * randN2, @@ -399,7 +438,7 @@ void cpuscatter( int n_q, cpu_kernel( n_q, q_x, q_y, q_z, n_atoms, r_x, r_y, r_z, - n_atom_types, atom_types, cromermann, + n_atom_types, atom_types, cromermann, U, n_rotations, randN1, randN2, randN3, q_out_real, q_out_imag ); } diff --git a/src/scatter/cpp_scatter.cu b/src/scatter/cpp_scatter.cu index e2c2e4f..ccd09f9 100644 --- a/src/scatter/cpp_scatter.cu +++ b/src/scatter/cpp_scatter.cu @@ -64,6 +64,7 @@ void __global__ gpu_kernel(int const n_q, int const n_atom_types, int const * const __restrict__ atom_types, float const * const __restrict__ cromermann, + float const * const __restrict__ U, int const n_rotations, float const * const __restrict__ q0, @@ -91,7 +92,7 @@ void __global__ gpu_kernel(int const n_q, float ax, ay, az; // rotated r vector float mq, qo, fi; // mag of q, formfactor for atom i float qr; // dot product of q and r - + float qUq; // matrix product of qT * U_ii * q while(gid < n_q) { @@ -143,8 +144,10 @@ void __global__ gpu_kernel(int const n_q, qr = ax*qx + ay*qy + az*qz; - q_sum.x += fi*__sinf(qr); - q_sum.y += fi*__cosf(qr); + qUq_product(U, a, qx, qy, qz, qUq); + + q_sum.x += fi*__sinf(qr) * exp(- 0.5 * qUq); + q_sum.y += fi*__cosf(qr) * exp(- 0.5 * qUq); } // finished one atom (3rd loop) } // finished one molecule (2nd loop) @@ -324,6 +327,9 @@ void _gpuscatter(int device_id, int * h_atom_types, float * h_cromermann, + // atomic displacement parameters + float * h_U, + // random numbers for rotations int n_rotations, float * rand1, @@ -371,7 +377,7 @@ void _gpuscatter(int device_id, const unsigned int id_size = n_atoms * sizeof(int); const unsigned int cm_size = 9 * n_atom_types * sizeof(float); const unsigned int quat_size = n_rotations * sizeof(float); - + const unsigned int U_size = n_atoms * 9 * sizeof(float); err = cudaGetLastError(); if (err != cudaSuccess) { @@ -390,7 +396,8 @@ void _gpuscatter(int device_id, int *d_id; deviceMalloc( (void **) &d_id, id_size); float *d_cm; deviceMalloc( (void **) &d_cm, cm_size); - + float *d_U; deviceMalloc( (void **) &d_U, U_size); + float *d_q0; deviceMalloc( (void **) &d_q0, quat_size); float *d_q1; deviceMalloc( (void **) &d_q1, quat_size); float *d_q2; deviceMalloc( (void **) &d_q2, quat_size); @@ -428,7 +435,8 @@ void _gpuscatter(int device_id, cudaMemcpy(d_id, &h_atom_types[0], id_size, cudaMemcpyHostToDevice); cudaMemcpy(d_cm, &h_cromermann[0], cm_size, cudaMemcpyHostToDevice); - + cudaMemcpy(d_U, &h_U[0], U_size, cudaMemcpyHostToDevice); + cudaMemcpy(d_q0, &h_q0[0], quat_size, cudaMemcpyHostToDevice); cudaMemcpy(d_q1, &h_q1[0], quat_size, cudaMemcpyHostToDevice); cudaMemcpy(d_q2, &h_q2[0], quat_size, cudaMemcpyHostToDevice); @@ -444,7 +452,7 @@ void _gpuscatter(int device_id, // execute the kernel gpu_kernel <<>> (n_q, d_qx, d_qy, d_qz, n_atoms, d_rx, d_ry, d_rz, - n_atom_types, d_id, d_cm, + n_atom_types, d_id, d_cm, d_U, n_rotations, d_q0, d_q1, d_q2, d_q3, d_q_out_real, d_q_out_imag); cudaThreadSynchronize(); @@ -479,6 +487,7 @@ void _gpuscatter(int device_id, cudaFree(d_id); cudaFree(d_cm); + cudaFree(d_U); cudaFree(d_q0); cudaFree(d_q1); diff --git a/src/scatter/cpp_scatter.hh b/src/scatter/cpp_scatter.hh index 904d515..ada8771 100644 --- a/src/scatter/cpp_scatter.hh +++ b/src/scatter/cpp_scatter.hh @@ -20,6 +20,8 @@ void gpuscatter(int device_id, int * atom_types, float * cromermann, + float * U, + int n_rotations, float * randN1, float * randN2, @@ -44,6 +46,8 @@ void cpuscatter( int * atom_types, float * cromermann, + float * U, + int n_rotations, float * randN1, float * randN2, @@ -117,6 +121,9 @@ void _gpuscatter(int device_id, int * h_atom_types, float * h_cromermann, + // atomic displacement parameters + float * h_U, + // random numbers for rotations int n_rotations, float * rand1, diff --git a/src/scatter/cpp_scatter_wrap.pyx b/src/scatter/cpp_scatter_wrap.pyx index 4b4d0a0..73d5088 100644 --- a/src/scatter/cpp_scatter_wrap.pyx +++ b/src/scatter/cpp_scatter_wrap.pyx @@ -54,6 +54,8 @@ cdef extern from "cpp_scatter.hh": int * atom_types, float * cromermann, + float * U, + int n_rotations, float * randN1, float * randN2, @@ -76,6 +78,8 @@ cdef extern from "cpp_scatter.hh": int * atom_types, float * cromermann, + float * U, + int n_rotations, float * randN1, float * randN2, @@ -172,6 +176,7 @@ def parallel_cpp_scatter(n_molecules, nodes=[], devices=[], random_state=None, + U=None, ignore_gpu_check=False): """ Multi-threaded interface to `cpp_scatter`. Specify the devices the @@ -245,7 +250,7 @@ def parallel_cpp_scatter(n_molecules, continue print('CPU Thread %d :: %d shots' % (cpu_thread, num)) cpu_args = (num, rxyz, qxyz, atom_types, cromermann_parameters, - 'CPU', random_state) + 'CPU', random_state, U) t_cpu = Thread(target=t_fxn, args=cpu_args) t_cpu.start() threads.append(t_cpu) @@ -256,7 +261,7 @@ def parallel_cpp_scatter(n_molecules, continue print('GPU Device %d :: %d shots' % (gpu_device, num)) gpu_args = (num, rxyz, qxyz, atom_types, cromermann_parameters, - gpu_device, random_state) + gpu_device, random_state, U) t_gpu = Thread(target=t_fxn, args=gpu_args) t_gpu.start() threads.append(t_gpu) @@ -274,7 +279,8 @@ def cpp_scatter(n_molecules, np.ndarray atom_types, np.ndarray cromermann_parameters, device_id='CPU', - random_state=None): + random_state=None, + U=None): """ A python interface to the C++ and CUDA scattering code. The idea here is to mirror that interface closely, but in a pythonic fashion. @@ -307,7 +313,11 @@ def cpp_scatter(n_molecules, random_state : np.random.RandomState Seed the random state. For testing only. - + + U : ndarray, float + An n x 3 x 3 array of the mean-square atomic displacement parameters of + each atom. + Returns ------- amplitudes : ndarray, complex128 @@ -339,7 +349,7 @@ def cpp_scatter(n_molecules, cdef np.ndarray[ndim=2, dtype=np.float32_t, mode="c"] c_qxyz cdef np.ndarray[ndim=2, dtype=np.float32_t, mode="c"] c_rxyz cdef np.ndarray[ndim=1, dtype=np.float32_t] c_cromermann - + c_qxyz = np.ascontiguousarray(qxyz.T, dtype=np.float32) c_rxyz = np.ascontiguousarray(rxyz.T, dtype=np.float32) c_cromermann = np.ascontiguousarray(cromermann_parameters, dtype=np.float32) @@ -358,6 +368,13 @@ def cpp_scatter(n_molecules, cdef np.ndarray[ndim=2, dtype=np.float32_t, mode="c"] c_rfloats c_rfloats = np.ascontiguousarray( random_state.rand(3, n_molecules), dtype=np.float32 ) + + # generate U matrix of zeros if no ADPs input + num_atoms = rxyz.shape[0] + if U is None: + U = np.zeros((num_atoms, 3, 3)) + cdef np.ndarray[ndim=1, dtype=np.float32_t, mode="c"] c_U + c_U = np.ascontiguousarray(U.flatten(), dtype=np.float32) # initialize output arrays cdef np.ndarray[ndim=1, dtype=np.float32_t] real_amplitudes @@ -371,7 +388,7 @@ def cpp_scatter(n_molecules, if device_id == 'CPU': cpuscatter(qxyz.shape[0], &c_qxyz[0,0], &c_qxyz[1,0], &c_qxyz[2,0], rxyz.shape[0], &c_rxyz[0,0], &c_rxyz[1,0], &c_rxyz[2,0], - num_atom_types, &c_atom_types[0], &c_cromermann[0], + num_atom_types, &c_atom_types[0], &c_cromermann[0], &c_U[0], n_molecules, &c_rfloats[0,0], &c_rfloats[1,0], &c_rfloats[2,0], &real_amplitudes[0], &imag_amplitudes[0]) @@ -383,7 +400,7 @@ def cpp_scatter(n_molecules, gpuscatter(device_id, qxyz.shape[0], &c_qxyz[0,0], &c_qxyz[1,0], &c_qxyz[2,0], rxyz.shape[0], &c_rxyz[0,0], &c_rxyz[1,0], &c_rxyz[2,0], - num_atom_types, &c_atom_types[0], &c_cromermann[0], + num_atom_types, &c_atom_types[0], &c_cromermann[0], &c_U[0], n_molecules, &c_rfloats[0,0], &c_rfloats[1,0], &c_rfloats[2,0], &real_amplitudes[0], &imag_amplitudes[0]) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 1d378cd..7af1e45 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -402,6 +402,7 @@ def test_cpu_scatter(self): def test_cpu_scatter_isoU(self): cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + self.iso_U = self.iso_U[:,None,None] * np.eye(3)[None,None,:] print 'num_molecules:', self.num_molecules cpu_A_isoU = _cppscatter.cpp_scatter(self.num_molecules, From 835e9d595a2e2bc941ba44049999e73500b9c402 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Thu, 30 Nov 2017 16:14:23 -0800 Subject: [PATCH 07/15] added tests of updated CppScatter Python interface --- test/scatter/test_scatter.py | 93 ++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 7af1e45..6415fdf 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -536,6 +536,99 @@ def test_parallel_interface(self): # err_msg='error in 2 thread cpu') +class TestCppScatterPython(object): + + def setup(self): + self.device = 'CPU' + + def test_python_interface_noU(self): + + nq = 100 # number of detector vectors to do + q_grid = np.loadtxt(ref_file('512_q.xyz'))[:nq] + + traj = mdtraj.load(ref_file('ala2.pdb')) + atomic_numbers = np.array([ a.element.atomic_number for a in traj.topology.atoms ]) + rxyz = traj.xyz[0] * 10.0 + + ref_A = ref_simulate_shot(rxyz, + atomic_numbers, + 1, + q_grid, + dont_rotate=True) + + cpu_A = scatter.simulate_atomic(traj, + 1, + q_grid, + ignore_hydrogens=False, + dont_rotate=True) + + assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, + err_msg='scatter: c-cpu/cpu reference mismatch') + assert not np.all( cpu_A == 0.0 ) + assert not np.sum( cpu_A == np.nan ) + + def test_python_interface_isoU(self): + + nq = 100 # number of detector vectors to do + q_grid = np.loadtxt(ref_file('512_q.xyz'))[:nq] + + traj = mdtraj.load(ref_file('ala2.pdb')) + atomic_numbers = np.array([ a.element.atomic_number for a in traj.topology.atoms ]) + rxyz = traj.xyz[0] * 10.0 + + iso_U = np.random.rand(rxyz.shape[0]) / 10.0 + + ref_A = ref_simulate_shot(rxyz, + atomic_numbers, + 1, + q_grid, + dont_rotate=True, + U=iso_U) + + cpu_A = scatter.simulate_atomic(traj, + 1, + q_grid, + ignore_hydrogens=False, + dont_rotate=True, + U=iso_U) + + assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, + err_msg='scatter: c-cpu/cpu reference mismatch') + assert not np.all( cpu_A == 0.0 ) + assert not np.sum( cpu_A == np.nan ) + + def test_python_interface_anisoU(self): + + nq = 100 # number of detector vectors to do + q_grid = np.loadtxt(ref_file('512_q.xyz'))[:nq] + + traj = mdtraj.load(ref_file('ala2.pdb')) + atomic_numbers = np.array([ a.element.atomic_number for a in traj.topology.atoms ]) + rxyz = traj.xyz[0] * 10.0 + + aniso_U = np.random.rand(rxyz.shape[0], 3, 3) + for i in range(rxyz.shape[0]): + aniso_U[i] += aniso_U[i].T + aniso_U /= 20.0 + + ref_A = ref_simulate_shot(rxyz, + atomic_numbers, + 1, + q_grid, + dont_rotate=True, + U=aniso_U) + + cpu_A = scatter.simulate_atomic(traj, + 1, + q_grid, + ignore_hydrogens=False, + dont_rotate=True, + U=aniso_U) + + assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, + err_msg='scatter: c-cpu/cpu reference mismatch') + assert not np.all( cpu_A == 0.0 ) + assert not np.sum( cpu_A == np.nan ) class TestDiffuseScatter(object): From 80a693f491b25ff719b60e7e2d9c684f63427c09 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Thu, 30 Nov 2017 16:51:19 -0800 Subject: [PATCH 08/15] revised nosetests; now passes GPU tests --- test/scatter/test_scatter.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 6415fdf..e72014b 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -402,7 +402,7 @@ def test_cpu_scatter(self): def test_cpu_scatter_isoU(self): cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) - self.iso_U = self.iso_U[:,None,None] * np.eye(3)[None,None,:] + iso_U = self.iso_U[:,None,None] * np.eye(3)[None,None,:] print 'num_molecules:', self.num_molecules cpu_A_isoU = _cppscatter.cpp_scatter(self.num_molecules, @@ -412,7 +412,7 @@ def test_cpu_scatter_isoU(self): cromermann_parameters, device_id='CPU', random_state=np.random.RandomState(RANDOM_SEED), - U=self.iso_U) + U=iso_U) assert_allclose(cpu_A_isoU, self.ref_A_isoU, rtol=1e-3, atol=1.0, err_msg='scatter: c-cpu/cpu reference mismatch with isotropic B factors') @@ -463,6 +463,7 @@ def test_gpu_scatter_isoU(self): if not GPU: raise SkipTest cromermann_parameters, atom_types = get_cromermann_parameters(self.atomic_numbers) + iso_U = self.iso_U[:,None,None] * np.eye(3)[None,None,:] print 'num_molecules:', self.num_molecules gpu_A_isoU = _cppscatter.cpp_scatter(self.num_molecules, @@ -472,7 +473,7 @@ def test_gpu_scatter_isoU(self): cromermann_parameters, device_id=0, random_state=np.random.RandomState(RANDOM_SEED), - U=self.iso_U) + U=iso_U) assert_allclose(gpu_A_isoU, self.ref_A_isoU, rtol=1e-3, atol=1.0, err_msg='scatter: cuda-gpu/cpu reference mismatch with isotropic B factors') From 3b83ba2c8739d294ab372715bb37564eaadd40eb Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Thu, 30 Nov 2017 17:47:33 -0800 Subject: [PATCH 09/15] updated error messages --- test/scatter/test_scatter.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index e72014b..67e2199 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -564,7 +564,7 @@ def test_python_interface_noU(self): dont_rotate=True) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, - err_msg='scatter: c-cpu/cpu reference mismatch') + err_msg='scatter: python/cpu reference mismatch') assert not np.all( cpu_A == 0.0 ) assert not np.sum( cpu_A == np.nan ) @@ -594,7 +594,7 @@ def test_python_interface_isoU(self): U=iso_U) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, - err_msg='scatter: c-cpu/cpu reference mismatch') + err_msg='scatter: python/cpu reference mismatch') assert not np.all( cpu_A == 0.0 ) assert not np.sum( cpu_A == np.nan ) @@ -627,7 +627,7 @@ def test_python_interface_anisoU(self): U=aniso_U) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, - err_msg='scatter: c-cpu/cpu reference mismatch') + err_msg='scatter: python/cpu reference mismatch') assert not np.all( cpu_A == 0.0 ) assert not np.sum( cpu_A == np.nan ) From 33f0c558a61868650e9ce0e85a671afeafafe44c Mon Sep 17 00:00:00 2001 From: tjlane Date: Sat, 2 Dec 2017 17:58:05 -0800 Subject: [PATCH 10/15] center finding bug --- src/python/parse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/parse.py b/src/python/parse.py index dd9efbd..e3a5ce5 100644 --- a/src/python/parse.py +++ b/src/python/parse.py @@ -894,7 +894,7 @@ def find_center(image2d, mask=None, initial_guess=None, pix_res=0.1, window=15, x_size = image2d.shape[0] y_size = image2d.shape[1] - if mask != None: + if mask is not None: if not mask.shape == image2d.shape: raise ValueError('Mask and image must have same shape! Got %s, %s' ' respectively.' % ( str(mask.shape), str(image2d.shape) )) From 2ae9763f46d6e4eb783a3ad0aedf84466181b6c2 Mon Sep 17 00:00:00 2001 From: tjlane Date: Sat, 2 Dec 2017 18:08:03 -0800 Subject: [PATCH 11/15] added meaningless benchmark test --- src/scatter/cpp_scatter.cpp | 151 +++++++++++++++++++----------------- 1 file changed, 78 insertions(+), 73 deletions(-) diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index f2ad636..e8a38b1 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -284,7 +284,7 @@ void cpu_kernel( int const n_q, int const * const __restrict__ atom_types, float const * const __restrict__ cromermann, - float const * const __restrict__ U, + float const * const __restrict__ U, int const n_rotations, float const * const __restrict__ randN1, @@ -425,7 +425,7 @@ void cpuscatter( int n_q, int * atom_types, float * cromermann, - float * U, + float * U, int n_rotations, float * randN1, @@ -580,77 +580,13 @@ void cpudiffuse( int n_q, } // This is a meaningless test, for speed only... -// #ifndef __CUDACC__ -// int main() { -// -// int nQ_ = 1000; -// int nAtoms_ = 1000; -// int n_atom_types_ = 10; -// int nRot_ = 1000; -// -// float * h_qx_ = new float[nQ_]; -// float * h_qy_ = new float[nQ_]; -// float * h_qz_ = new float[nQ_]; -// -// float * h_rx_ = new float[nAtoms_]; -// float * h_ry_ = new float[nAtoms_]; -// float * h_rz_ = new float[nAtoms_]; -// -// int * atom_types_ = new int[nAtoms_]; -// float * cromermann_ = new float[n_atom_types_ * 9]; -// -// float * h_rand1_ = new float[nRot_]; -// float * h_rand2_ = new float[nRot_]; -// float * h_rand3_ = new float[nRot_]; -// -// float * h_outQ_R = new float[nQ_]; -// float * h_outQ_I = new float[nQ_]; -// -// cpuscatter ( // q vectors -// nQ_, -// h_qx_, -// h_qy_, -// h_qz_, -// -// // atomic positions, ids -// nAtoms_, -// h_rx_, -// h_ry_, -// h_rz_, -// -// // formfactor info -// n_atom_types_, -// atom_types_, -// cromermann_, -// -// // random numbers for rotations -// nRot_, -// h_rand1_, -// h_rand2_, -// h_rand3_, -// -// // output -// h_outQ_R, -// h_outQ_I ); -// -// printf("CPP OUTPUT:\n"); -// printf("%f\n", h_outQ_R[0]); -// printf("%f\n", h_outQ_I[0]); -// -// return 0; -// } -// #endif - #ifndef __CUDACC__ int main() { - int nQ_ = 100; - int nAtoms_ = 1500; - - std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; - std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; - + int nQ_ = 1000; + int nAtoms_ = 1000; int n_atom_types_ = 10; + int nRot_ = 1000; float * h_qx_ = new float[nQ_]; float * h_qy_ = new float[nQ_]; @@ -662,13 +598,17 @@ int main() { int * atom_types_ = new int[nAtoms_]; float * cromermann_ = new float[n_atom_types_ * 9]; + + float * U_ = new float[nAtoms_ * 3]; - float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; + float * h_rand1_ = new float[nRot_]; + float * h_rand2_ = new float[nRot_]; + float * h_rand3_ = new float[nRot_]; float * h_outQ_R = new float[nQ_]; float * h_outQ_I = new float[nQ_]; - cpudiffuse ( // q vectors + cpuscatter ( // q vectors nQ_, h_qx_, h_qy_, @@ -684,9 +624,15 @@ int main() { n_atom_types_, atom_types_, cromermann_, + + // ADPs + U_, - // correlation matrix - V, + // random numbers for rotations + nRot_, + h_rand1_, + h_rand2_, + h_rand3_, // output h_outQ_R, @@ -699,3 +645,62 @@ int main() { return 0; } #endif + +// #ifndef __CUDACC__ +// int main() { +// +// int nQ_ = 100; +// int nAtoms_ = 1500; +// +// std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; +// std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; +// +// int n_atom_types_ = 10; +// +// float * h_qx_ = new float[nQ_]; +// float * h_qy_ = new float[nQ_]; +// float * h_qz_ = new float[nQ_]; +// +// float * h_rx_ = new float[nAtoms_]; +// float * h_ry_ = new float[nAtoms_]; +// float * h_rz_ = new float[nAtoms_]; +// +// int * atom_types_ = new int[nAtoms_]; +// float * cromermann_ = new float[n_atom_types_ * 9]; +// +// float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; +// +// float * h_outQ_R = new float[nQ_]; +// float * h_outQ_I = new float[nQ_]; +// +// cpudiffuse ( // q vectors +// nQ_, +// h_qx_, +// h_qy_, +// h_qz_, +// +// // atomic positions, ids +// nAtoms_, +// h_rx_, +// h_ry_, +// h_rz_, +// +// // formfactor info +// n_atom_types_, +// atom_types_, +// cromermann_, +// +// // correlation matrix +// V, +// +// // output +// h_outQ_R, +// h_outQ_I ); +// +// printf("CPP OUTPUT:\n"); +// printf("%f\n", h_outQ_R[0]); +// printf("%f\n", h_outQ_I[0]); +// +// return 0; +// } +// #endif From 43d1e244288c41f64d4696e4cd6fb075901aa615 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Mon, 4 Dec 2017 21:05:10 -0800 Subject: [PATCH 12/15] modified adp calculation for slight speed improvement on CPU --- src/scatter/cpp_scatter.cpp | 130 ++++++++++++------------------------ 1 file changed, 41 insertions(+), 89 deletions(-) diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index f2ad636..771c90e 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -324,10 +324,11 @@ void cpu_kernel( int const n_q, float mq, qo, fi; // mag of q, formfactor for atom i float q_sum_real, q_sum_imag; // partial sum of real and imaginary amplitude float qr; // dot product of q and r - float qUq; // Debye Waller factor, qT * U_ii * q + // float qUq; // Debye Waller factor, qT * U_ii * q // we will use a small array to store form factors float * formfactors = (float *) malloc(n_atom_types * sizeof(float)); + float * qUq = (float *) malloc(n_atoms * sizeof(float)); // pre-compute rotation quaternions float * q0 = (float *) malloc(n_rotations * sizeof(float)); @@ -372,6 +373,13 @@ void cpu_kernel( int const n_q, } + // precompute Debye-Waller factors for each atom + for( int a = 0; a < n_atoms; a++ ) { + + qUq_product(U, a, qx, qy, qz, qUq[a]); + + } + // for each molecule (2nd nested loop) for( int im = 0; im < n_rotations; im++ ) { @@ -388,11 +396,10 @@ void cpu_kernel( int const n_q, ax, ay, az); qr = ax*qx + ay*qy + az*qz; - qUq_product(U, a, qx, qy, qz, qUq); // FIXME :: swap cos and sin???? e^i*t = cos(t) + i sin(t) - q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq); - q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq); + q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq[a]); + q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq[a]); } // finished one atom (3rd loop) } // finished one molecule (2nd loop) @@ -400,7 +407,7 @@ void cpu_kernel( int const n_q, // add the output to the total intensity array q_out_real[iq] = q_sum_real; q_out_imag[iq] = q_sum_imag; - + } // finished one q vector (1st loop) free(formfactors); @@ -408,6 +415,7 @@ void cpu_kernel( int const n_q, free(q1); free(q2); free(q3); + free(qUq); } @@ -580,122 +588,66 @@ void cpudiffuse( int n_q, } // This is a meaningless test, for speed only... -// #ifndef __CUDACC__ -// int main() { -// -// int nQ_ = 1000; -// int nAtoms_ = 1000; -// int n_atom_types_ = 10; -// int nRot_ = 1000; -// -// float * h_qx_ = new float[nQ_]; -// float * h_qy_ = new float[nQ_]; -// float * h_qz_ = new float[nQ_]; -// -// float * h_rx_ = new float[nAtoms_]; -// float * h_ry_ = new float[nAtoms_]; -// float * h_rz_ = new float[nAtoms_]; -// -// int * atom_types_ = new int[nAtoms_]; -// float * cromermann_ = new float[n_atom_types_ * 9]; -// -// float * h_rand1_ = new float[nRot_]; -// float * h_rand2_ = new float[nRot_]; -// float * h_rand3_ = new float[nRot_]; -// -// float * h_outQ_R = new float[nQ_]; -// float * h_outQ_I = new float[nQ_]; -// -// cpuscatter ( // q vectors -// nQ_, -// h_qx_, -// h_qy_, -// h_qz_, -// -// // atomic positions, ids -// nAtoms_, -// h_rx_, -// h_ry_, -// h_rz_, -// -// // formfactor info -// n_atom_types_, -// atom_types_, -// cromermann_, -// -// // random numbers for rotations -// nRot_, -// h_rand1_, -// h_rand2_, -// h_rand3_, -// -// // output -// h_outQ_R, -// h_outQ_I ); -// -// printf("CPP OUTPUT:\n"); -// printf("%f\n", h_outQ_R[0]); -// printf("%f\n", h_outQ_I[0]); -// -// return 0; -// } -// #endif - #ifndef __CUDACC__ int main() { - - int nQ_ = 100; - int nAtoms_ = 1500; - - std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; - std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; - +// + int nQ_ = 1000; + int nAtoms_ = 1000; int n_atom_types_ = 10; - + int nRot_ = 1000; +// float * h_qx_ = new float[nQ_]; float * h_qy_ = new float[nQ_]; float * h_qz_ = new float[nQ_]; - +// float * h_rx_ = new float[nAtoms_]; float * h_ry_ = new float[nAtoms_]; float * h_rz_ = new float[nAtoms_]; - +// int * atom_types_ = new int[nAtoms_]; float * cromermann_ = new float[n_atom_types_ * 9]; - - float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; - + float * U = new float[nAtoms_ * 3 * 3]; +// + float * h_rand1_ = new float[nRot_]; + float * h_rand2_ = new float[nRot_]; + float * h_rand3_ = new float[nRot_]; +// float * h_outQ_R = new float[nQ_]; float * h_outQ_I = new float[nQ_]; - - cpudiffuse ( // q vectors +// + cpuscatter ( // q vectors nQ_, h_qx_, h_qy_, h_qz_, - +// // atomic positions, ids nAtoms_, h_rx_, h_ry_, h_rz_, - +// // formfactor info n_atom_types_, atom_types_, cromermann_, - - // correlation matrix - V, - + U, +// + // random numbers for rotations + nRot_, + h_rand1_, + h_rand2_, + h_rand3_, +// // output h_outQ_R, h_outQ_I ); - +// printf("CPP OUTPUT:\n"); printf("%f\n", h_outQ_R[0]); printf("%f\n", h_outQ_I[0]); - +// return 0; } #endif + From f8010a29e9909f02326129b0c868ac8fa981f35e Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Mon, 4 Dec 2017 21:19:43 -0800 Subject: [PATCH 13/15] modified adp calculation, slight speed-up on CPU --- src/scatter/cpp_scatter.cpp | 84 +++++++++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 13 deletions(-) diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index 771c90e..df165dd 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -590,64 +590,122 @@ void cpudiffuse( int n_q, // This is a meaningless test, for speed only... #ifndef __CUDACC__ int main() { -// + int nQ_ = 1000; int nAtoms_ = 1000; int n_atom_types_ = 10; int nRot_ = 1000; -// + float * h_qx_ = new float[nQ_]; float * h_qy_ = new float[nQ_]; float * h_qz_ = new float[nQ_]; -// + float * h_rx_ = new float[nAtoms_]; float * h_ry_ = new float[nAtoms_]; float * h_rz_ = new float[nAtoms_]; -// + int * atom_types_ = new int[nAtoms_]; float * cromermann_ = new float[n_atom_types_ * 9]; float * U = new float[nAtoms_ * 3 * 3]; -// + float * h_rand1_ = new float[nRot_]; float * h_rand2_ = new float[nRot_]; float * h_rand3_ = new float[nRot_]; -// + float * h_outQ_R = new float[nQ_]; float * h_outQ_I = new float[nQ_]; -// + cpuscatter ( // q vectors nQ_, h_qx_, h_qy_, h_qz_, -// + // atomic positions, ids nAtoms_, h_rx_, h_ry_, h_rz_, -// + // formfactor info n_atom_types_, atom_types_, cromermann_, U, -// + // random numbers for rotations nRot_, h_rand1_, h_rand2_, h_rand3_, -// + // output h_outQ_R, h_outQ_I ); -// + printf("CPP OUTPUT:\n"); printf("%f\n", h_outQ_R[0]); printf("%f\n", h_outQ_I[0]); -// + return 0; } #endif +// #ifndef __CUDACC__ +// int main() { +// +// int nQ_ = 100; +// int nAtoms_ = 1500; + +// std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; +// std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; + +// int n_atom_types_ = 10; + +// float * h_qx_ = new float[nQ_]; +// float * h_qy_ = new float[nQ_]; +// float * h_qz_ = new float[nQ_]; + +// float * h_rx_ = new float[nAtoms_]; +// float * h_ry_ = new float[nAtoms_]; +// float * h_rz_ = new float[nAtoms_]; + +// int * atom_types_ = new int[nAtoms_]; +// float * cromermann_ = new float[n_atom_types_ * 9]; + +// float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; + +// float * h_outQ_R = new float[nQ_]; +// float * h_outQ_I = new float[nQ_]; + +// cpudiffuse ( // q vectors +// nQ_, +// h_qx_, +// h_qy_, +// h_qz_, + + // atomic positions, ids +// nAtoms_, +// h_rx_, +// h_ry_, +// h_rz_, + + // formfactor info +// n_atom_types_, +// atom_types_, +// cromermann_, + + // correlation matrix +// V, + + // output +// h_outQ_R, +// h_outQ_I ); + +// printf("CPP OUTPUT:\n"); +// printf("%f\n", h_outQ_R[0]); +// printf("%f\n", h_outQ_I[0]); + +// return 0; +// } +// #endif From bbbc5fb0eac4738f12320163a94bc31d2013038b Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 5 Dec 2017 09:41:54 -0800 Subject: [PATCH 14/15] Revert to 3b83ba2 --- src/scatter/cpp_scatter.cpp | 162 +++++++++++++++++------------------- 1 file changed, 76 insertions(+), 86 deletions(-) diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index df165dd..f2ad636 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -324,11 +324,10 @@ void cpu_kernel( int const n_q, float mq, qo, fi; // mag of q, formfactor for atom i float q_sum_real, q_sum_imag; // partial sum of real and imaginary amplitude float qr; // dot product of q and r - // float qUq; // Debye Waller factor, qT * U_ii * q + float qUq; // Debye Waller factor, qT * U_ii * q // we will use a small array to store form factors float * formfactors = (float *) malloc(n_atom_types * sizeof(float)); - float * qUq = (float *) malloc(n_atoms * sizeof(float)); // pre-compute rotation quaternions float * q0 = (float *) malloc(n_rotations * sizeof(float)); @@ -373,13 +372,6 @@ void cpu_kernel( int const n_q, } - // precompute Debye-Waller factors for each atom - for( int a = 0; a < n_atoms; a++ ) { - - qUq_product(U, a, qx, qy, qz, qUq[a]); - - } - // for each molecule (2nd nested loop) for( int im = 0; im < n_rotations; im++ ) { @@ -396,10 +388,11 @@ void cpu_kernel( int const n_q, ax, ay, az); qr = ax*qx + ay*qy + az*qz; + qUq_product(U, a, qx, qy, qz, qUq); // FIXME :: swap cos and sin???? e^i*t = cos(t) + i sin(t) - q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq[a]); - q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq[a]); + q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq); + q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq); } // finished one atom (3rd loop) } // finished one molecule (2nd loop) @@ -407,7 +400,7 @@ void cpu_kernel( int const n_q, // add the output to the total intensity array q_out_real[iq] = q_sum_real; q_out_imag[iq] = q_sum_imag; - + } // finished one q vector (1st loop) free(formfactors); @@ -415,7 +408,6 @@ void cpu_kernel( int const n_q, free(q1); free(q2); free(q3); - free(qUq); } @@ -588,13 +580,77 @@ void cpudiffuse( int n_q, } // This is a meaningless test, for speed only... +// #ifndef __CUDACC__ +// int main() { +// +// int nQ_ = 1000; +// int nAtoms_ = 1000; +// int n_atom_types_ = 10; +// int nRot_ = 1000; +// +// float * h_qx_ = new float[nQ_]; +// float * h_qy_ = new float[nQ_]; +// float * h_qz_ = new float[nQ_]; +// +// float * h_rx_ = new float[nAtoms_]; +// float * h_ry_ = new float[nAtoms_]; +// float * h_rz_ = new float[nAtoms_]; +// +// int * atom_types_ = new int[nAtoms_]; +// float * cromermann_ = new float[n_atom_types_ * 9]; +// +// float * h_rand1_ = new float[nRot_]; +// float * h_rand2_ = new float[nRot_]; +// float * h_rand3_ = new float[nRot_]; +// +// float * h_outQ_R = new float[nQ_]; +// float * h_outQ_I = new float[nQ_]; +// +// cpuscatter ( // q vectors +// nQ_, +// h_qx_, +// h_qy_, +// h_qz_, +// +// // atomic positions, ids +// nAtoms_, +// h_rx_, +// h_ry_, +// h_rz_, +// +// // formfactor info +// n_atom_types_, +// atom_types_, +// cromermann_, +// +// // random numbers for rotations +// nRot_, +// h_rand1_, +// h_rand2_, +// h_rand3_, +// +// // output +// h_outQ_R, +// h_outQ_I ); +// +// printf("CPP OUTPUT:\n"); +// printf("%f\n", h_outQ_R[0]); +// printf("%f\n", h_outQ_I[0]); +// +// return 0; +// } +// #endif + #ifndef __CUDACC__ int main() { - int nQ_ = 1000; - int nAtoms_ = 1000; + int nQ_ = 100; + int nAtoms_ = 1500; + + std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; + std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; + int n_atom_types_ = 10; - int nRot_ = 1000; float * h_qx_ = new float[nQ_]; float * h_qy_ = new float[nQ_]; @@ -606,16 +662,13 @@ int main() { int * atom_types_ = new int[nAtoms_]; float * cromermann_ = new float[n_atom_types_ * 9]; - float * U = new float[nAtoms_ * 3 * 3]; - float * h_rand1_ = new float[nRot_]; - float * h_rand2_ = new float[nRot_]; - float * h_rand3_ = new float[nRot_]; + float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; float * h_outQ_R = new float[nQ_]; float * h_outQ_I = new float[nQ_]; - cpuscatter ( // q vectors + cpudiffuse ( // q vectors nQ_, h_qx_, h_qy_, @@ -631,13 +684,9 @@ int main() { n_atom_types_, atom_types_, cromermann_, - U, - // random numbers for rotations - nRot_, - h_rand1_, - h_rand2_, - h_rand3_, + // correlation matrix + V, // output h_outQ_R, @@ -650,62 +699,3 @@ int main() { return 0; } #endif - -// #ifndef __CUDACC__ -// int main() { -// -// int nQ_ = 100; -// int nAtoms_ = 1500; - -// std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; -// std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; - -// int n_atom_types_ = 10; - -// float * h_qx_ = new float[nQ_]; -// float * h_qy_ = new float[nQ_]; -// float * h_qz_ = new float[nQ_]; - -// float * h_rx_ = new float[nAtoms_]; -// float * h_ry_ = new float[nAtoms_]; -// float * h_rz_ = new float[nAtoms_]; - -// int * atom_types_ = new int[nAtoms_]; -// float * cromermann_ = new float[n_atom_types_ * 9]; - -// float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; - -// float * h_outQ_R = new float[nQ_]; -// float * h_outQ_I = new float[nQ_]; - -// cpudiffuse ( // q vectors -// nQ_, -// h_qx_, -// h_qy_, -// h_qz_, - - // atomic positions, ids -// nAtoms_, -// h_rx_, -// h_ry_, -// h_rz_, - - // formfactor info -// n_atom_types_, -// atom_types_, -// cromermann_, - - // correlation matrix -// V, - - // output -// h_outQ_R, -// h_outQ_I ); - -// printf("CPP OUTPUT:\n"); -// printf("%f\n", h_outQ_R[0]); -// printf("%f\n", h_outQ_I[0]); - -// return 0; -// } -// #endif From ae7ef6d5c25d5f9bd129cc4a86dae7c7b2991925 Mon Sep 17 00:00:00 2001 From: Ariana Peck Date: Tue, 5 Dec 2017 13:25:41 -0800 Subject: [PATCH 15/15] moved adp calculation, slight speed-up on CPU --- src/scatter/cpp_scatter.cpp | 161 ++++++++++++++++++----------------- test/scatter/test_scatter.py | 14 +-- 2 files changed, 92 insertions(+), 83 deletions(-) diff --git a/src/scatter/cpp_scatter.cpp b/src/scatter/cpp_scatter.cpp index f2ad636..5bff417 100644 --- a/src/scatter/cpp_scatter.cpp +++ b/src/scatter/cpp_scatter.cpp @@ -324,10 +324,10 @@ void cpu_kernel( int const n_q, float mq, qo, fi; // mag of q, formfactor for atom i float q_sum_real, q_sum_imag; // partial sum of real and imaginary amplitude float qr; // dot product of q and r - float qUq; // Debye Waller factor, qT * U_ii * q - // we will use a small array to store form factors + // we will use a small array to store form factors and Debye-Waller factors float * formfactors = (float *) malloc(n_atom_types * sizeof(float)); + float * qUq = (float *) malloc(n_atoms * sizeof(float)); // pre-compute rotation quaternions float * q0 = (float *) malloc(n_rotations * sizeof(float)); @@ -372,6 +372,11 @@ void cpu_kernel( int const n_q, } + // precompute Debye-Waller factors for each atom + for( int a = 0; a < n_atoms; a++ ) { + qUq_product(U, a, qx, qy, qz, qUq[a]); + } + // for each molecule (2nd nested loop) for( int im = 0; im < n_rotations; im++ ) { @@ -388,11 +393,10 @@ void cpu_kernel( int const n_q, ax, ay, az); qr = ax*qx + ay*qy + az*qz; - qUq_product(U, a, qx, qy, qz, qUq); // FIXME :: swap cos and sin???? e^i*t = cos(t) + i sin(t) - q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq); - q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq); + q_sum_real += fi * sinf(qr) * exp(- 0.5 * qUq[a]); + q_sum_imag += fi * cosf(qr) * exp(- 0.5 * qUq[a]); } // finished one atom (3rd loop) } // finished one molecule (2nd loop) @@ -408,6 +412,7 @@ void cpu_kernel( int const n_q, free(q1); free(q2); free(q3); + free(qUq); } @@ -580,77 +585,13 @@ void cpudiffuse( int n_q, } // This is a meaningless test, for speed only... -// #ifndef __CUDACC__ -// int main() { -// -// int nQ_ = 1000; -// int nAtoms_ = 1000; -// int n_atom_types_ = 10; -// int nRot_ = 1000; -// -// float * h_qx_ = new float[nQ_]; -// float * h_qy_ = new float[nQ_]; -// float * h_qz_ = new float[nQ_]; -// -// float * h_rx_ = new float[nAtoms_]; -// float * h_ry_ = new float[nAtoms_]; -// float * h_rz_ = new float[nAtoms_]; -// -// int * atom_types_ = new int[nAtoms_]; -// float * cromermann_ = new float[n_atom_types_ * 9]; -// -// float * h_rand1_ = new float[nRot_]; -// float * h_rand2_ = new float[nRot_]; -// float * h_rand3_ = new float[nRot_]; -// -// float * h_outQ_R = new float[nQ_]; -// float * h_outQ_I = new float[nQ_]; -// -// cpuscatter ( // q vectors -// nQ_, -// h_qx_, -// h_qy_, -// h_qz_, -// -// // atomic positions, ids -// nAtoms_, -// h_rx_, -// h_ry_, -// h_rz_, -// -// // formfactor info -// n_atom_types_, -// atom_types_, -// cromermann_, -// -// // random numbers for rotations -// nRot_, -// h_rand1_, -// h_rand2_, -// h_rand3_, -// -// // output -// h_outQ_R, -// h_outQ_I ); -// -// printf("CPP OUTPUT:\n"); -// printf("%f\n", h_outQ_R[0]); -// printf("%f\n", h_outQ_I[0]); -// -// return 0; -// } -// #endif - #ifndef __CUDACC__ int main() { - int nQ_ = 100; - int nAtoms_ = 1500; - - std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; - std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; - + int nQ_ = 1000; + int nAtoms_ = 1000; int n_atom_types_ = 10; + int nRot_ = 1000; float * h_qx_ = new float[nQ_]; float * h_qy_ = new float[nQ_]; @@ -662,13 +603,16 @@ int main() { int * atom_types_ = new int[nAtoms_]; float * cromermann_ = new float[n_atom_types_ * 9]; + float * U = new float[nAtoms_ * 3 * 3]; - float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; + float * h_rand1_ = new float[nRot_]; + float * h_rand2_ = new float[nRot_]; + float * h_rand3_ = new float[nRot_]; float * h_outQ_R = new float[nQ_]; float * h_outQ_I = new float[nQ_]; - cpudiffuse ( // q vectors + cpuscatter ( // q vectors nQ_, h_qx_, h_qy_, @@ -685,8 +629,14 @@ int main() { atom_types_, cromermann_, - // correlation matrix - V, + // ADP info + U, + + // random numbers for rotations + nRot_, + h_rand1_, + h_rand2_, + h_rand3_, // output h_outQ_R, @@ -699,3 +649,62 @@ int main() { return 0; } #endif + +// #ifndef __CUDACC__ +// int main() { +// +// int nQ_ = 100; +// int nAtoms_ = 1500; +// +// std::cout << nQ_ << " q-vectors :: " << nAtoms_ << " atoms\n"; +// std::cout << "remember: linear in q-vectors, quadratic in atoms\n"; +// +// int n_atom_types_ = 10; +// +// float * h_qx_ = new float[nQ_]; +// float * h_qy_ = new float[nQ_]; +// float * h_qz_ = new float[nQ_]; +// +// float * h_rx_ = new float[nAtoms_]; +// float * h_ry_ = new float[nAtoms_]; +// float * h_rz_ = new float[nAtoms_]; +// +// int * atom_types_ = new int[nAtoms_]; +// float * cromermann_ = new float[n_atom_types_ * 9]; +// +// float * V = new float[nAtoms_ * nAtoms_ * 3 * 3]; +// +// float * h_outQ_R = new float[nQ_]; +// float * h_outQ_I = new float[nQ_]; +// +// cpudiffuse ( // q vectors +// nQ_, +// h_qx_, +// h_qy_, +// h_qz_, +// +// // atomic positions, ids +// nAtoms_, +// h_rx_, +// h_ry_, +// h_rz_, +// +// // formfactor info +// n_atom_types_, +// atom_types_, +// cromermann_, +// +// // correlation matrix +// V, +// +// // output +// h_outQ_R, +// h_outQ_I ); +// +// printf("CPP OUTPUT:\n"); +// printf("%f\n", h_outQ_R[0]); +// printf("%f\n", h_outQ_I[0]); +// +// return 0; +// } +// #endif diff --git a/test/scatter/test_scatter.py b/test/scatter/test_scatter.py index 67e2199..31627d9 100644 --- a/test/scatter/test_scatter.py +++ b/test/scatter/test_scatter.py @@ -550,18 +550,18 @@ def test_python_interface_noU(self): traj = mdtraj.load(ref_file('ala2.pdb')) atomic_numbers = np.array([ a.element.atomic_number for a in traj.topology.atoms ]) rxyz = traj.xyz[0] * 10.0 + rs = np.random.RandomState(RANDOM_SEED) ref_A = ref_simulate_shot(rxyz, atomic_numbers, 1, - q_grid, - dont_rotate=True) + q_grid) cpu_A = scatter.simulate_atomic(traj, 1, q_grid, ignore_hydrogens=False, - dont_rotate=True) + random_state=rs,) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, err_msg='scatter: python/cpu reference mismatch') @@ -578,19 +578,19 @@ def test_python_interface_isoU(self): rxyz = traj.xyz[0] * 10.0 iso_U = np.random.rand(rxyz.shape[0]) / 10.0 + rs = np.random.RandomState(RANDOM_SEED) ref_A = ref_simulate_shot(rxyz, atomic_numbers, 1, q_grid, - dont_rotate=True, U=iso_U) cpu_A = scatter.simulate_atomic(traj, 1, q_grid, ignore_hydrogens=False, - dont_rotate=True, + random_state=rs, U=iso_U) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0, @@ -606,6 +606,7 @@ def test_python_interface_anisoU(self): traj = mdtraj.load(ref_file('ala2.pdb')) atomic_numbers = np.array([ a.element.atomic_number for a in traj.topology.atoms ]) rxyz = traj.xyz[0] * 10.0 + rs = np.random.RandomState(RANDOM_SEED) aniso_U = np.random.rand(rxyz.shape[0], 3, 3) for i in range(rxyz.shape[0]): @@ -616,14 +617,13 @@ def test_python_interface_anisoU(self): atomic_numbers, 1, q_grid, - dont_rotate=True, U=aniso_U) cpu_A = scatter.simulate_atomic(traj, 1, q_grid, ignore_hydrogens=False, - dont_rotate=True, + random_state=rs, U=aniso_U) assert_allclose(cpu_A, ref_A, rtol=1e-3, atol=1.0,