[WIP] Diffuse

[tjlane]

Adding C++ code for computing the MVN diffuse scatter model.

$ time ./cputest
10000 q-vectors :: 1000 atoms
remember: linear in q-vectors, quadratic in atomsCPP OUTPUT:
0.000000
0.000000

real 0m57.507s
user 0m56.895s
sys 0m0.313s

[tjlane]

@apeck12 here is how this works

You can access the python interface to the new code by doing

from thor import _cpuscatter
_cppscatter.cpp_scatter_diffuse(xyzlist, q_grid, atom_types, cromermann_parameters, V)

The (cython-provided) python interface is here:
https://github.com/tjlane/thor/blob/diffuse/src/scatter/cpp_scatter_wrap.pyx#L384

The C++ code implementing the (almost certainly incorrect) correlation matrix is here:
https://github.com/tjlane/thor/blob/diffuse/src/scatter/cpp_scatter.cpp#L444
https://github.com/tjlane/thor/blob/diffuse/src/scatter/cpp_scatter.cpp#L459

Note that V is a huge array representing the tensor V_ij^ab. In python, as indicated in the docstring in cpp_scatter_wrap.pyx,

        V : ndarray, float
            A 4-d array of shape (n_atoms, n_atoms, 3, 3) representing the
            anisotropic Gaussian correlation between atoms. The value
            V[i,j,a,b] is the correlation between atom i in the a direction with
            j in the b direction, where a & b are one of {x/y/z}.

In C-land, this will be accessed as a flattened array, so V[i,j,a,b] in python should correspond to V[9*(num_atoms * i + j) + 3*a + b]. At least that's what I think -- double check me! And feel free to fix up that part of the c-code.

[tjlane]

@apeck12 regarding testing & what works now:

What works: If you pass a V array of all zeros, right now you do get the same result as a reference implementation. This is tested here:
https://github.com/tjlane/thor/blob/diffuse/test/scatter/test_scatter.py#L368

Adding new test: Add your reference implementation to the top of the test_scatter.py file and then add a new test in the TestDiffuseScatter class that uses the V tensor.

If you add stuff and it doesn't work right off the bat, push the new test an I can help work to make the code correct.

[tjlane]

TODO:

• test correctness of code/implement test
• provide a high level python interface
• benchmark to see if we meet performance requirements

[tjlane]

@apeck12 benchmark for a "large" system and one q-vector

$ time ./cputest 
1 q-vectors :: 15000 atoms
remember: linear in q-vectors, quadratic in atoms
CPP OUTPUT:
0.000000
0.000000

real	0m4.248s
user	0m1.917s
sys	0m1.855s

For 10,000,000 q-vectors, we're talking ~11,000 cpu-hours, which is too many! So we need to keep thinking. Note parallelization over q-vecs is trivial.

[tjlane]

Going to the isotropic approximation seems to make little difference. Going to GPU is possible/easy.

[apeck12]

I've been using the GPU code exclusively for the Thor CypA and ECR simulations. Do you mind if I switch the test case for TestDiffuseScatter from the 512_atom_benchmark.xyz to pentagon.pdb?

[apeck12]

Also, is there a Cython wrapper for the calculation with isotropic V?

[tjlane]

@apeck12

1. 512_atom_benchmark.xyz to pentagon.pdb: go for it
2. There is not (and probably will never be) a specific isotropic interface. I committed a sketch python interface just now (see scatter.py, commit above). Basically the trick will be to simply fill in a sparse V matrix to do the isotropic case. If we can get a big speedup from specific isotropic code, we will do that, but my tests so far indicate that is not the case.

[apeck12]

@tjlane

Python reference implementation is operational in pull request #31 .

For the pentagon tests, _cppscatter.cpp_scatter_diffuse

• produces correct output when V=0, as you had shown for 512_atom_benchmark.xyz.
• fails for the test non-zero V matrix, with error: 'Infs detected in scattering output!'

[tjlane]

@apeck12 cool. I will work on this and fix things up. Thanks!

[tjlane]

@apeck12 latest commit includes functional C++ code splitting diffuse/bragg. Speed improvements to come, but you can start using it if you want. Current best estimate for 1M q-vectors & 1500 atoms is ~30,000 cpu-seconds (8 hrs).

Since we are in a rapidly changing development mode, but have decent tests, I would recommend checking out the latest copy and then always running the tests (at least test_scatter.py) before doing calculations we want to trust.

[apeck12]

@tjlane
Receiving the following error when trying to run scatter.simulate_diffuse:
Error in kernel. CUDA error: invalid argument

This is my call to the Thor code:
scatter.simulate_diffuse(pdb, detector, V, ignore_hydrogens=True, device_id=0)

Any thoughts on what might be up? I think it should be the same version of cuda that I was using on the master branch without issue.

[tjlane]

@apeck12
For a quick fix, just make sure ignore_hydrogens=False does not fix the problem.
If that doesn't work, can you send me the entire example so I can replicate the problem?