cctbx.miller.complete_array() bug?
Hi, I think I found a bug or missing feature in cctbx.miller.complete_array: As I suppose the complete_array function should check if a set of reflections is complete and add missing indices with the specified values. The bug is now when a set contains symmetrically equivalent reflections e.g. in P4 (1, 0, 2) and (0, 1, 2). The sample code below generates this output for the atached data file: --- f_obs read in: Miller array info: None Observation type: xray.amplitude Type of data: double, size=8 Type of sigmas: None Number of Miller indices: 8 Anomalous flag: False Unit cell: (4, 4, 4.44, 90, 90, 90) Space group: P 4 (No. 75) Systematic absences: 0 Centric reflections: 3 Resolution range: 4.44 1.94109 Completeness in resolution range: 1 Completeness with d_max=infinity: 1 (0, 0, 1) 4.10852771683 (0, 1, 0) 4.34856298103 (0, 1, 1) 10.0 (1, 1, 0) 0.911043357914 (1, 1, 1) 9.12523972288 (0, 0, 2) 5.97578446733 (0, 2, 0) 7.10492786733 (1, 0, 2) 2.17485631709 --- f_obs completed: Miller array info: None Observation type: xray.amplitude Type of data: double, size=9 Type of sigmas: None Number of Miller indices: 9 Anomalous flag: False Unit cell: (4, 4, 4.44, 90, 90, 90) Space group: P 4 (No. 75) Systematic absences: 0 Centric reflections: 3 Resolution range: 4.44 1.94109 (0, 0, 1) 4.10852771683 (0, 1, 0) 4.34856298103 (0, 1, 1) 10.0 (1, 1, 0) 0.911043357914 (1, 1, 1) 9.12523972288 (0, 0, 2) 5.97578446733 (0, 2, 0) 7.10492786733 (1, 0, 2) 2.17485631709 (0, 1, 2) 0.0 ---- sample code ------- obs_file = "Model02.hkl" xtal_symm = crystal.symmetry( unit_cell=(4.000, 4.0, 4.44, 90.0, 90.0, 90.0), space_group_symbol="P 4") # load (integral) intensities from diffraction data from iotbx import reflection_file_reader rfl = reflection_file_reader.any_reflection_file("intensities="+obs_file, ensure_read_access=True) I_obs = rfl.as_miller_arrays(crystal_symmetry=xtal_symm, force_symmetry=False, merge_equivalents=True, base_array_info=None)[0] I_obs = I_obs.discard_sigmas() I_obs.merge_equivalents() f_obs = I_obs.as_amplitude_array() print("\n--- f_obs read in:") f_obs.show_comprehensive_summary() f_obs.show_array() d_max, d_min = f_obs.resolution_range() # add missing reflections with 0.0 intensity (= not observed) f_obs = f_obs.complete_array(d_min=d_min, d_max=d_max, new_data_value=0.0, new_sigmas_value=1.0).sort() print("\n--- f_obs completed:") f_obs.show_comprehensive_summary() f_obs.show_array() ---- end sample code ------- As you can see the index (0, 1, 2) is added with zero amplitude where instead the index (1, 0, 2) should be renamed to (0, 1, 2) as those two are symmetrically equivalent for any tetragonal spacegroup (like P4). So to me it looks like the complete_array method does not honor symmetrically equivalent indices at all, where it should imho. What do you think on this issue? With regards, Dipl. Phys. Jan M. Simons Institute of Crystallography RWTH Aachen University
Hi Jan,
As you can see the index (0, 1, 2) is added with zero amplitude where instead the index (1, 0, 2) should be renamed to (0, 1, 2) as those two are symmetrically equivalent for any tetragonal spacegroup (like P4). So to me it looks like the complete_array method does not honor symmetrically equivalent indices at all, where it should imho. What do you think on this issue?
Considering the current implementation of method complete_array, it seems to me that the problem lies in method complete_set of class miller.set, that replaces (1, 0, 2) with (0, 1, 2) in your example for no particularly good reason I can see. The implementation of complete_array clearly assumes that complete_set should not do that. But then, the way complete_set is implemented, such surprises are bound to happen. A more involved algorithm is necessary there. I put that on the todo list. Best wishes, Luc Bourhis
Hi,
In general, I use a 'map_to_asu' whenever possible, this avoids some
problems you illustrate, but isn't elegant.
P
On Thu, Jul 12, 2012 at 8:37 AM, Luc Bourhis
Hi Jan,
As you can see the index (0, 1, 2) is added with zero amplitude where instead the index (1, 0, 2) should be renamed to (0, 1, 2) as those two are symmetrically equivalent for any tetragonal spacegroup (like P4). So to me it looks like the complete_array method does not honor symmetrically equivalent indices at all, where it should imho. What do you think on this issue?
Considering the current implementation of method complete_array, it seems to me that the problem lies in method complete_set of class miller.set, that replaces (1, 0, 2) with (0, 1, 2) in your example for no particularly good reason I can see. The implementation of complete_array clearly assumes that complete_set should not do that. But then, the way complete_set is implemented, such surprises are bound to happen. A more involved algorithm is necessary there.
I put that on the todo list.
Best wishes,
Luc Bourhis
_______________________________________________ cctbxbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/cctbxbb
Hi Peter and Nat,
On Thu, Jul 12, 2012 at 8:45 AM, Petrus Zwart
wrote: In general, I use a 'map_to_asu' whenever possible, this avoids some problems you illustrate, but isn't elegant.
+1
This is pretty essential for much of our code, assuming you're working with merged/symmetry-unique data.
It would make sense to automagically map to asu in complete_set and complete_array in order to avoid surprises, wouldn't it? Would you see any adverse effect? Best wishes, Luc Bourhis
I floated this once to Ralf and his objection was performance. I guess
it also depends on the situation if you want this to happen.
P
On 12 July 2012 23:06, Luc Bourhis
Hi Peter and Nat,
On Thu, Jul 12, 2012 at 8:45 AM, Petrus Zwart
wrote: In general, I use a 'map_to_asu' whenever possible, this avoids some problems you illustrate, but isn't elegant.
+1
This is pretty essential for much of our code, assuming you're working with merged/symmetry-unique data.
It would make sense to automagically map to asu in complete_set and complete_array in order to avoid surprises, wouldn't it? Would you see any adverse effect?
Best wishes,
Luc Bourhis
_______________________________________________ cctbxbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/cctbxbb
-- ----------------------------------------------------------------- P.H. Zwart Research Scientist Berkeley Center for Structural Biology Lawrence Berkeley National Laboratories 1 Cyclotron Road, Berkeley, CA-94703, USA Cell: 510 289 9246 BCSB: http://bcsb.als.lbl.gov PHENIX: http://www.phenix-online.org SASTBX: http://sastbx.als.lbl.gov -----------------------------------------------------------------
Yes, Ralf's concern was performance. If you hide away map_to_asu in a function that you or someone else runs repetitively then one can run into performance surprises. So it might be better to call map_to_asu consciously, when it's really needed. Pavel On 7/12/12 11:11 PM, Peter Zwart wrote:
I floated this once to Ralf and his objection was performance. I guess it also depends on the situation if you want this to happen.
P
On 12 July 2012 23:06, Luc Bourhis
wrote: Hi Peter and Nat,
On Thu, Jul 12, 2012 at 8:45 AM, Petrus Zwart
wrote: In general, I use a 'map_to_asu' whenever possible, this avoids some problems you illustrate, but isn't elegant. +1
This is pretty essential for much of our code, assuming you're working with merged/symmetry-unique data. It would make sense to automagically map to asu in complete_set and complete_array in order to avoid surprises, wouldn't it? Would you see any adverse effect?
Best wishes,
Luc Bourhis
_______________________________________________ cctbxbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/cctbxbb
Maybe it would be best just to add a comment to the complete_array and
complete_set functions to make clear that the implementation assumes that
map_to_asu has been called previously?
Richard
On 12 July 2012 23:23, Pavel Afonine
Yes, Ralf's concern was performance. If you hide away map_to_asu in a function that you or someone else runs repetitively then one can run into performance surprises. So it might be better to call map_to_asu consciously, when it's really needed.
Pavel
On 7/12/12 11:11 PM, Peter Zwart wrote:
I floated this once to Ralf and his objection was performance. I guess it also depends on the situation if you want this to happen.
P
On 12 July 2012 23:06, Luc Bourhis
wrote: Hi Peter and Nat,
On Thu, Jul 12, 2012 at 8:45 AM, Petrus Zwart
wrote: In general, I use a 'map_to_asu' whenever possible, this avoids some problems you illustrate, but isn't elegant.
+1
This is pretty essential for much of our code, assuming you're working with merged/symmetry-unique data.
It would make sense to automagically map to asu in complete_set and complete_array in order to avoid surprises, wouldn't it? Would you see any adverse effect?
Best wishes,
Luc Bourhis
______________________________**_________________ cctbxbb mailing list [email protected] http://phenix-online.org/**mailman/listinfo/cctbxbbhttp://phenix-online.org/mailman/listinfo/cctbxbb
______________________________**_________________ cctbxbb mailing list [email protected] http://phenix-online.org/**mailman/listinfo/cctbxbbhttp://phenix-online.org/mailman/listinfo/cctbxbb
participants (7)
-
Jan Marten Simons
-
Luc Bourhis
-
Nathaniel Echols
-
Pavel Afonine
-
Peter Zwart
-
Petrus Zwart
-
Richard Gildea