Accommodate absorption in a diamond anvil cell

[benjaminhwilliams]

Diamond anvil cells are commonly used to allow measurement structures under applied pressure. Both the incident and diffracted X-ray beam are attenuated by absorption in the anvils between which the sample is sandwiched. The path length of the incident beam through the diamond, and hence the amount of absorption, depends on the orientation of the cell to the incident beam. The path length of the diffracted beam depends on the orientation of the cell and the diffracted angle. Hence a non-uniform absorption factor applies to reflection intensities, which affects the fidelity of data reduction using the existing DIALS scaling tools.

By modelling the anvils as flat plates, of equal thickness and with parallel sides, the overall absorption factor can be calculated. It is a function of the incident and diffracted beam directions, the thickness of each anvil and their normal (itself a function of the mounting orientation and goniometer rotation) and the linear attenuation coefficient of the anvil material (itself a function of wavelength and the material's density).

We can thereby correct the integrated intensities for the anvil absorption and then apply the post-integration components of the data reduction pipeline as normal, just as though the sample had been measured in air.

For diamond anvil cell data, this tool should be called on the output of dials.integrate before any subsequent processing. The default anvil parameters reflect the cell design used on I19-2 at Diamond Light Source.

[benjaminhwilliams]

There are obviously no tests included here. This comes back to earlier discussions about the best way to incorporate regression data for small molecule work. My preference would be to merge this, test it in the school of hard knocks (I19-2) and make adding tests a separate task.

[rjgildea]

For this pull request you shouldn't need a whole dataset added to dials-data - just the integrated experiments/reflections should be enough to exercise this command. To test correct_intensities_for_dac_attenuation you could probably even generate some fake experiments/reflections on the fly for the purposes of exercising the function.

[graeme-winter]

Agree with comments from @rjgildea on tests - they can run on very small data sets, anything high pressure will be fine, just need the integrated reflection files. Looking through the changes now.

[graeme-winter]

Not a precondition on the merge but a proposal for the future - pretty sure everything done in here could easily have been achieved with flex calculations rather than scipy/numpy stuff - which would make it slightly more in line with the rest of the dials code base, though admittedly slightly harder for the "developer in the street"

[graeme-winter]

Modulo comments re: testing and a handful of small suggestions is a useful improvement thank you.

[jbeilstenedmands]

Looks good overall. Wasn't really able to follow the maths just from reading but trust you've got this sorted :) One thing regarding multiple reflection tables, I notice you join them all together at the start, then go through selecting each one to apply the corrections individually.
I think we've shared our differing views on this issue before, but in my opinion its better to split the input reflections into a list of single-dataset tables, then have a function which runs on a single experiment and single-dataset reflection table, as the algorithm is basically doing in this case (or have a function working on an ExperimentList and list of reflection tables the same length). The tables can then be joined together at the end for output. This was probably makes it easier for the algorithm code to be more 'functional' in my opinion.

[benjaminhwilliams]

Not a precondition on the merge but a proposal for the future - pretty sure everything done in here could easily have been achieved with flex calculations rather than scipy/numpy stuff - which would make it slightly more in line with the rest of the dials code base, though admittedly slightly harder for the "developer in the street"

I wasn't sure if there was an efficient way in flex to perform some of the arithmetic here, such as the dot product of each of an array of vectors with each of another array of vectors. This may simply be flex-ignorance on my part but sometimes life's just too short when the alternative is easy, quick and staring me in the face.

[benjaminhwilliams]

Looks good overall. Wasn't really able to follow the maths just from reading but trust you've got this sorted :)

Yes, I started trying to write it out in the PR description and then remembered that maths in GitHub isn't easy. Perhaps it would be better in the wiki or something. I'll finish writing it up so people can refer to it in future.

One thing regarding multiple reflection tables, I notice you join them all together at the start, then go through selecting each one to apply the corrections individually.
I think we've shared our differing views on this issue before, but in my opinion its better to split the input reflections into a list of single-dataset tables, then have a function which runs on a single experiment and single-dataset reflection table, as the algorithm is basically doing in this case (or have a function working on an ExperimentList and list of reflection tables the same length). The tables can then be joined together at the end for output. This was probably makes it easier for the algorithm code to be more 'functional' in my opinion.

I may be coming round to your way of thinking on this. I'll see about refactoring it slightly.

[phyy-nx]

Hi, just to help out a bit :) from dials.array_family import flex a = flex.vec3_double(10, (1,0,0)) b = flex.vec3_double(10, (0,1,0)) a.dot(b) list(a.dot(b)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] Also IPython's ? feature is helpful here. I haven't read the code, but changing between flex and numpy (using as_numpy_array() for example) does do a copy, so I try to avoid it as much as possible. If I start in flex or numpy I try to stay there. -Aaron

…

On Tue, Jan 14, 2020 at 8:28 AM Richard Gildea ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In command_line/rescale_diamond_anvil_cell.py <#1090 (comment)>: > @@ -0,0 +1,302 @@ +# coding: utf-8 "correct" is potentially ambiguous in that it can be an adjective or a verb, whereas at least "rescale" is unambiguous in that sense... — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1090?email_source=notifications&email_token=ADGY5STYFHY3FHDVXEMKM6LQ5XR3LA5CNFSM4KGFW6R2YY3PNVWWK3TUL52HS4DFWFIHK3DMKJSXC5LFON2FEZLWNFSXPKTDN5WW2ZLOORPWSZGOCRWK3UA#discussion_r366440469>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADGY5SVTLTWHNK67T357UQDQ5XR3LANCNFSM4KGFW6RQ> .

[benjaminhwilliams]

Hi, just to help out a bit :)

Ah, thanks @phyy-nx. And if I want to apply each of an array of rotation operators to each of an array of vectors, is there a flexy way to do that?

[phyy-nx]

Yep. I'd do it like this: from dials.array_family import flex from scitbx import matrix a = flex.vec3_double(10, (1,0,0)) r = matrix.col((0,1,0)).axis_and_angle_as_r3_rotation_matrix(90, deg=True) b = (flex.mat3_double(10, r)*a) print(b[0]) (2.220446049250313e-16, 0.0, -1.0) The slight numerical instability there is interesting of course. I haven't tried compiling with fast math off and no optimisations, but I wonder if that would clean it up.

…

On Tue, Jan 14, 2020 at 10:12 AM Ben Williams ***@***.***> wrote: Hi, just to help out a bit :) Ah, thanks Aaron. And if I want to apply each of an array of rotation operators to each of an array of vectors, is there a flexy way to do that? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#1090?email_source=notifications&email_token=ADGY5SVT4KNOY7QLQKK5IDDQ5X56PA5CNFSM4KGFW6R2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEI5SWJA#issuecomment-574303012>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADGY5SWPC2M5ROPQO66P3JTQ5X56PANCNFSM4KGFW6RQ> .

[benjaminhwilliams]

Looks good overall. Wasn't really able to follow the maths just from reading but trust you've got this sorted :)

The changes to the docs in this branch contain details of the arithmetic, in case you want to look that over.

[benjaminhwilliams]

@rjgildea, this now has tests. Does it have your approval?

[rjgildea]

I note that this code makes heavy use of numpy/scipy routines. I don't know whether this is something we want to encourage/discourage, particularly where there are alternative ways to accomplish the same using flex arrays. I appreciate that you probably already know how to accomplish certain things with numpy arrays, whereas you may not know how to do the same with flex arrays.

[Anthchirp]

humbly suggest moving the .dxf and .svg files into the dials.github.io repository (same as other images)

[benjaminhwilliams]

humbly suggest moving the .dxf and .svg files into the dials.github.io repository (same as other images)

Done.

[benjaminhwilliams]

I note that this code makes heavy use of numpy/scipy routines. I don't know whether this is something we want to encourage/discourage, particularly where there are alternative ways to accomplish the same using flex arrays. I appreciate that you probably already know how to accomplish certain things with numpy arrays, whereas you may not know how to do the same with flex arrays.

I leave this to future discussion. I appreciate the desire to keep things uniform, and I know @ndevenish has opinions about this, but I reckon that's another problem for another time.

[benjaminhwilliams]

It looks like the only failures in the Travis builds are Python 3 build tests that seem to be completely unaffected by this PR. Otherwise I'm ready to merge...

[Anthchirp]

Some sort of smtbx build failure, no clear idea why.
Travis on long-running branches can be a problem as it doesn't necessarily have the correct cctbx repository checkout for the time you branched off master.

If you rebased your branch on master this will likely be fixed.