Computing E^4 from an mtz file
Hi Folks, Trying to write some code to return the average E^4 across (some) bins for an mtz file, I got to def main(): import sys from iotbx import reflection_file_reader for argv in sys.argv[1:]: reflection_file = reflection_file_reader.any_reflection_file( file_name = argv) if reflection_file.file_type() is None: continue miller_arrays = reflection_file.as_miller_arrays() for ma in miller_arrays: E4(ma) def E4(ma): f = ma.as_intensity_array() f = f.array(data=f.data()/f.epsilons().data().as_double()) f.set_observation_type_xray_intensity() f.setup_binner(auto_binning = True) sm = f.second_moment(use_binning = True) sm.show() if __name__ == '__main__': main() but I have no idea how to actually extract the information I want from sm as this only has show() which gives Graemes-MacBook-Pro:PlayPen graeme$ cctbx.python E4.py /Users/graeme/projects/sci-397/DEFAULT/scale/AUTOMATIC_DEFAULT_scaled.mtz unused: - 20.8668 [ 0/7 ] bin 1: 20.8668 - 5.8873 [221/235] 2.3863 bin 2: 5.8873 - 4.6904 [216/221] 1.9823 bin 3: 4.6904 - 4.1026 [221/222] 1.9925 <snip> bin 47: 1.6553 - 1.6435 [153/209] 2.2051 bin 48: 1.6435 - 1.6320 [143/198] 2.2021 bin 49: 1.6320 - 1.6208 [141/204] 2.7085 unused: 1.6208 - [ 0/0 ] unused: - 20.8668 [ 0/9 ] bin 1: 20.8668 - 7.1654 [201/225] 2.2405 bin 2: 7.1654 - 5.7261 [202/213] 1.9109 bin 3: 5.7261 - 5.0137 [216/223] 2.0162 bin 4: 5.0137 - 4.5605 [207/212] 2.0043 <snip> bin 87: 1.6456 - 1.6393 [126/192] 2.5521 bin 88: 1.6393 - 1.6330 [122/210] 2.0195 bin 89: 1.6330 - 1.6269 [131/224] 2.6251 bin 90: 1.6269 - 1.6208 [120/212] 3.0912 unused: 1.6208 - [ 0/0 ] Really I would like the value programatically or at least a way of getting back an array of the binned values so I can average them myself (which is obviously trivial) Any suggestions? Thanks, Graeme ps - most of this code is ripped from iotbx.reflection_statistics :o) -- This e-mail and any attachments may contain confidential, copyright and or privileged material, and are for the use of the intended addressee only. If you are not the intended addressee or an authorised recipient of the addressee please notify us of receipt by returning the e-mail and do not use, copy, retain, distribute or disclose the information in or attached to the e-mail. Any opinions expressed within this e-mail are those of the individual and not necessarily of Diamond Light Source Ltd. Diamond Light Source Ltd. cannot guarantee that this e-mail or any attachments are free from viruses and we cannot accept liability for any damage which you may sustain as a result of software viruses which may be transmitted in or with the message. Diamond Light Source Limited (company no. 4375679). Registered in England and Wales with its registered office at Diamond House, Harwell Science and Innovation Campus, Didcot, Oxfordshire, OX11 0DE, United Kingdom
On Thu, Dec 13, 2012 at 3:00 AM,
Trying to write some code to return the average E^4 across (some) bins for an mtz file, I got to ... but I have no idea how to actually extract the information I want from sm as this only has show() which gives ... Really I would like the value programatically or at least a way of getting back an array of the binned values so I can average them myself (which is obviously trivial)
Won't simply passing use_binning=False to second_moment() do what you want? It returns a float in this mode (which is common to most methods which have a use_binning parameter). If you want to obtain the value for a subset of bins you can use array.select() to pull out the bin(s) and then call second_moment(use_binning=False). Or am I completely misunderstanding the question? -Nat
Hi Nat, That gave me six and a bit: specifically sm.show() print f.second_moment(use_binning = False) 6.46063790107 The E^2 = I/<I> need to be computed in bins. After that (E^2)^2 should be fine. This got me confused - it must be in there somewhere Thanks, Graeme On 13 Dec 2012, at 16:01, Nathaniel Echols wrote:
On Thu, Dec 13, 2012 at 3:00 AM,
wrote: Trying to write some code to return the average E^4 across (some) bins for an mtz file, I got to ... but I have no idea how to actually extract the information I want from sm as this only has show() which gives ... Really I would like the value programatically or at least a way of getting back an array of the binned values so I can average them myself (which is obviously trivial)
Won't simply passing use_binning=False to second_moment() do what you want? It returns a float in this mode (which is common to most methods which have a use_binning parameter). If you want to obtain the value for a subset of bins you can use array.select() to pull out the bin(s) and then call second_moment(use_binning=False).
Or am I completely misunderstanding the question?
-Nat _______________________________________________ cctbxbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/cctbxbb
-- This e-mail and any attachments may contain confidential, copyright and or privileged material, and are for the use of the intended addressee only. If you are not the intended addressee or an authorised recipient of the addressee please notify us of receipt by returning the e-mail and do not use, copy, retain, distribute or disclose the information in or attached to the e-mail. Any opinions expressed within this e-mail are those of the individual and not necessarily of Diamond Light Source Ltd. Diamond Light Source Ltd. cannot guarantee that this e-mail or any attachments are free from viruses and we cannot accept liability for any damage which you may sustain as a result of software viruses which may be transmitted in or with the message. Diamond Light Source Limited (company no. 4375679). Registered in England and Wales with its registered office at Diamond House, Harwell Science and Innovation Campus, Didcot, Oxfordshire, OX11 0DE, United Kingdom
On Thu, Dec 13, 2012 at 8:18 AM,
That gave me six and a bit: specifically
sm.show() print f.second_moment(use_binning = False)
6.46063790107
The E^2 = I/<I> need to be computed in bins. After that (E^2)^2 should be fine. This got me confused - it must be in there somewhere
Ah, I think I see - look at the "data" attribute: sm = f.second_moment(use_binning=True) print sm. data [None, 2.3279309836355044, 2.8099112859261517, 2.1769967857271175, 2.200637763782145, 2.2919308734359407, 2.000356888782398, 1.7555453224108684, 1.4836662688931992, None] (using a dataset from the PDB) I've never fully understood what the first and last elements are supposed to mean (at least one is for "unused" data) - but you can get the elements of interested this way: for i_bin in f.binner().range_used() : # etc. -Nat
Thanks Nat, that hit the spot. For reference the final code I got was def main(): import sys from iotbx import reflection_file_reader for argv in sys.argv[1:]: reflection_file = reflection_file_reader.any_reflection_file( file_name = argv) if reflection_file.file_type() is None: continue miller_arrays = reflection_file.as_miller_arrays() for ma in miller_arrays: print ma.info().label_string(), E4(ma) def E4(ma): f = ma.as_intensity_array() f = f.array(data=f.data()/f.epsilons().data().as_double()) f.set_observation_type_xray_intensity() f.setup_binner(auto_binning = True) sm = f.second_moment(use_binning = True) moments = [sm.data[j] for j in f.binner().range_used()] return sum(moments) / len(moments) if __name__ == '__main__': main() Best wishes, Graeme On 13 Dec 2012, at 16:29, Nathaniel Echols wrote:
On Thu, Dec 13, 2012 at 8:18 AM,
wrote: That gave me six and a bit: specifically
sm.show() print f.second_moment(use_binning = False)
6.46063790107
The E^2 = I/<I> need to be computed in bins. After that (E^2)^2 should be fine. This got me confused - it must be in there somewhere
Ah, I think I see - look at the "data" attribute:
sm = f.second_moment(use_binning=True) print sm. data [None, 2.3279309836355044, 2.8099112859261517, 2.1769967857271175, 2.200637763782145, 2.2919308734359407, 2.000356888782398, 1.7555453224108684, 1.4836662688931992, None]
(using a dataset from the PDB) I've never fully understood what the first and last elements are supposed to mean (at least one is for "unused" data) - but you can get the elements of interested this way:
for i_bin in f.binner().range_used() : # etc.
-Nat _______________________________________________ cctbxbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/cctbxbb
-- This e-mail and any attachments may contain confidential, copyright and or privileged material, and are for the use of the intended addressee only. If you are not the intended addressee or an authorised recipient of the addressee please notify us of receipt by returning the e-mail and do not use, copy, retain, distribute or disclose the information in or attached to the e-mail. Any opinions expressed within this e-mail are those of the individual and not necessarily of Diamond Light Source Ltd. Diamond Light Source Ltd. cannot guarantee that this e-mail or any attachments are free from viruses and we cannot accept liability for any damage which you may sustain as a result of software viruses which may be transmitted in or with the message. Diamond Light Source Limited (company no. 4375679). Registered in England and Wales with its registered office at Diamond House, Harwell Science and Innovation Campus, Didcot, Oxfordshire, OX11 0DE, United Kingdom
participants (2)
-
Graeme.Winter@diamond.ac.uk
-
Nathaniel Echols