- Timestamp:
- Apr 6, 2017 11:04:42 AM (8 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- a470e88
- Parents:
- 00a6ec4
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/dataloader/manipulations.py
rdd11014 rfd5d6eac 5 5 """ 6 6 ##################################################################### 7 # This software was developed by the University of Tennessee as part of the8 # Distributed Data Analysis of Neutron Scattering Experiments (DANSE)9 # project funded by the US National Science Foundation.10 # See the license text in license.txt11 # copyright 2008, University of Tennessee7 # This software was developed by the University of Tennessee as part of the 8 # Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 9 # project funded by the US National Science Foundation. 10 # See the license text in license.txt 11 # copyright 2008, University of Tennessee 12 12 ###################################################################### 13 13 14 #TODO: copy the meta data from the 2D object to the resulting 1D object 14 # If you want to run just a single test from this file: 15 # PYTHONPATH=../src/ python2 -m sasdataloader.test.utest_averaging data_info_tests.test_sectorq_full 16 # TODO: copy the meta data from the 2D object to the resulting 1D object 15 17 import math 16 import numpy 18 import numpy as np 17 19 18 20 #from data_info import plottable_2D … … 82 84 if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: 83 85 raise ValueError, "Can't convert this data: data=None..." 84 new_x = n umpy.tile(data2d.x_bins, (len(data2d.y_bins), 1))85 new_y = n umpy.tile(data2d.y_bins, (len(data2d.x_bins), 1))86 new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 87 new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) 86 88 new_y = new_y.swapaxes(0, 1) 87 89 … … 89 91 qx_data = new_x.flatten() 90 92 qy_data = new_y.flatten() 91 q_data = n umpy.sqrt(qx_data * qx_data + qy_data * qy_data)92 if data2d.err_data is None or n umpy.any(data2d.err_data <= 0):93 new_err_data = n umpy.sqrt(numpy.abs(new_data))93 q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 94 if data2d.err_data is None or np.any(data2d.err_data <= 0): 95 new_err_data = np.sqrt(np.abs(new_data)) 94 96 else: 95 97 new_err_data = data2d.err_data.flatten() 96 mask = n umpy.ones(len(new_data), dtype=bool)97 98 # TODO: make sense of the following two lines...98 mask = np.ones(len(new_data), dtype=bool) 99 100 # TODO: make sense of the following two lines... 99 101 #from sas.sascalc.dataloader.data_info import Data2D 100 102 #output = Data2D() … … 114 116 Compute average I(Q) for a region of interest 115 117 """ 118 116 119 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, 117 120 y_max=0.0, bin_width=0.001): … … 149 152 150 153 # Get data 151 data = data2D.data[n umpy.isfinite(data2D.data)]152 err_data = data2D.err_data[n umpy.isfinite(data2D.data)]153 qx_data = data2D.qx_data[n umpy.isfinite(data2D.data)]154 qy_data = data2D.qy_data[n umpy.isfinite(data2D.data)]154 data = data2D.data[np.isfinite(data2D.data)] 155 err_data = data2D.err_data[np.isfinite(data2D.data)] 156 qx_data = data2D.qx_data[np.isfinite(data2D.data)] 157 qy_data = data2D.qy_data[np.isfinite(data2D.data)] 155 158 156 159 # Build array of Q intervals … … 170 173 raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj) 171 174 172 x = n umpy.zeros(nbins)173 y = n umpy.zeros(nbins)174 err_y = n umpy.zeros(nbins)175 y_counts = n umpy.zeros(nbins)175 x = np.zeros(nbins) 176 y = np.zeros(nbins) 177 err_y = np.zeros(nbins) 178 y_counts = np.zeros(nbins) 176 179 177 180 # Average pixelsize in q space … … 205 208 continue 206 209 207 # TODO: find better definition of x[i_q] based on q_data210 # TODO: find better definition of x[i_q] based on q_data 208 211 # min_value + (i_q + 1) * self.bin_width / 2.0 209 212 x[i_q] += frac * q_value 210 213 y[i_q] += frac * data[npts] 211 214 212 if err_data ==None or err_data[npts] == 0.0:215 if err_data is None or err_data[npts] == 0.0: 213 216 if data[npts] < 0: 214 217 data[npts] = -data[npts] … … 225 228 y = y / y_counts 226 229 x = x / y_counts 227 idx = (n umpy.isfinite(y) & numpy.isfinite(x))230 idx = (np.isfinite(y) & np.isfinite(x)) 228 231 229 232 if not idx.any(): … … 237 240 Compute average I(Qy) for a region of interest 238 241 """ 242 239 243 def __call__(self, data2D): 240 244 """ … … 251 255 Compute average I(Qx) for a region of interest 252 256 """ 257 253 258 def __call__(self, data2D): 254 259 """ … … 264 269 Perform the sum of counts in a 2D region of interest. 265 270 """ 271 266 272 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 267 273 # Minimum Qx value [A-1] … … 304 310 raise RuntimeError, msg 305 311 # Get data 306 data = data2D.data[n umpy.isfinite(data2D.data)]307 err_data = data2D.err_data[n umpy.isfinite(data2D.data)]308 qx_data = data2D.qx_data[n umpy.isfinite(data2D.data)]309 qy_data = data2D.qy_data[n umpy.isfinite(data2D.data)]312 data = data2D.data[np.isfinite(data2D.data)] 313 err_data = data2D.err_data[np.isfinite(data2D.data)] 314 qx_data = data2D.qx_data[np.isfinite(data2D.data)] 315 qy_data = data2D.qy_data[np.isfinite(data2D.data)] 310 316 311 317 y = 0.0 … … 328 334 if self.y_min <= qy and self.y_max > qy: 329 335 frac_y = 1 330 # Find the fraction along each directions336 # Find the fraction along each directions 331 337 frac = frac_x * frac_y 332 338 if frac == 0: 333 339 continue 334 340 y += frac * data[npts] 335 if err_data ==None or err_data[npts] == 0.0:341 if err_data is None or err_data[npts] == 0.0: 336 342 if data[npts] < 0: 337 343 data[npts] = -data[npts] … … 347 353 Perform the average of counts in a 2D region of interest. 348 354 """ 355 349 356 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 350 357 super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, … … 398 405 as a function of Q 399 406 """ 407 400 408 def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 401 409 # Minimum radius included in the average [A-1] … … 414 422 """ 415 423 # Get data W/ finite values 416 data = data2D.data[n umpy.isfinite(data2D.data)]417 q_data = data2D.q_data[n umpy.isfinite(data2D.data)]418 err_data = data2D.err_data[n umpy.isfinite(data2D.data)]419 mask_data = data2D.mask[n umpy.isfinite(data2D.data)]424 data = data2D.data[np.isfinite(data2D.data)] 425 q_data = data2D.q_data[np.isfinite(data2D.data)] 426 err_data = data2D.err_data[np.isfinite(data2D.data)] 427 mask_data = data2D.mask[np.isfinite(data2D.data)] 420 428 421 429 dq_data = None … … 448 456 dq_overlap_y *= dq_overlap_y 449 457 450 dq_overlap = n umpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)458 dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 451 459 # Final protection of dq 452 460 if dq_overlap < 0: 453 461 dq_overlap = y_min 454 dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 455 dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 462 dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 463 dqy_data = data2D.dqy_data[np.isfinite( 464 data2D.data)] - dq_overlap 456 465 # def; dqx_data = dq_r dqy_data = dq_phi 457 466 # Convert dq 2D to 1D here 458 467 dqx = dqx_data * dqx_data 459 468 dqy = dqy_data * dqy_data 460 dq_data = n umpy.add(dqx, dqy)461 dq_data = n umpy.sqrt(dq_data)462 463 #q_data_max = n umpy.max(q_data)469 dq_data = np.add(dqx, dqy) 470 dq_data = np.sqrt(dq_data) 471 472 #q_data_max = np.max(q_data) 464 473 if len(data2D.q_data) == None: 465 474 msg = "Circular averaging: invalid q_data: %g" % data2D.q_data … … 469 478 nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 470 479 471 x = n umpy.zeros(nbins)472 y = n umpy.zeros(nbins)473 err_y = n umpy.zeros(nbins)474 err_x = n umpy.zeros(nbins)475 y_counts = n umpy.zeros(nbins)480 x = np.zeros(nbins) 481 y = np.zeros(nbins) 482 err_y = np.zeros(nbins) 483 err_x = np.zeros(nbins) 484 y_counts = np.zeros(nbins) 476 485 477 486 for npt in range(len(data)): … … 486 495 data_n = data[npt] 487 496 488 # #No need to calculate the frac when all data are within range497 # No need to calculate the frac when all data are within range 489 498 if self.r_min >= self.r_max: 490 499 raise ValueError, "Limit Error: min > max" … … 502 511 # Take dqs from data to get the q_average 503 512 x[i_q] += frac * q_value 504 if err_data ==None or err_data[npt] == 0.0:513 if err_data is None or err_data[npt] == 0.0: 505 514 if data_n < 0: 506 515 data_n = -data_n … … 523 532 err_y[n] = -err_y[n] 524 533 err_y[n] = math.sqrt(err_y[n]) 525 # if err_x != None:534 # if err_x != None: 526 535 # err_x[n] = math.sqrt(err_x[n]) 527 536 528 537 err_y = err_y / y_counts 529 err_y[err_y == 0] = n umpy.average(err_y)538 err_y[err_y == 0] = np.average(err_y) 530 539 y = y / y_counts 531 540 x = x / y_counts 532 idx = (n umpy.isfinite(y)) & (numpy.isfinite(x))541 idx = (np.isfinite(y)) & (np.isfinite(x)) 533 542 534 543 if err_x != None: … … 556 565 in anti-clockwise starting from the x- axis on the left-hand side 557 566 """ 558 #Todo: remove center. 567 # Todo: remove center. 568 559 569 def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): 560 570 # Minimum radius … … 569 579 self.nbins_phi = nbins 570 580 571 572 581 def __call__(self, data2D): 573 582 """ … … 585 594 586 595 # Get data 587 data = data2D.data[n umpy.isfinite(data2D.data)]588 q_data = data2D.q_data[n umpy.isfinite(data2D.data)]589 err_data = data2D.err_data[n umpy.isfinite(data2D.data)]590 qx_data = data2D.qx_data[n umpy.isfinite(data2D.data)]591 qy_data = data2D.qy_data[n umpy.isfinite(data2D.data)]596 data = data2D.data[np.isfinite(data2D.data)] 597 q_data = data2D.q_data[np.isfinite(data2D.data)] 598 err_data = data2D.err_data[np.isfinite(data2D.data)] 599 qx_data = data2D.qx_data[np.isfinite(data2D.data)] 600 qy_data = data2D.qy_data[np.isfinite(data2D.data)] 592 601 593 602 # Set space for 1d outputs 594 phi_bins = n umpy.zeros(self.nbins_phi)595 phi_counts = n umpy.zeros(self.nbins_phi)596 phi_values = n umpy.zeros(self.nbins_phi)597 phi_err = n umpy.zeros(self.nbins_phi)603 phi_bins = np.zeros(self.nbins_phi) 604 phi_counts = np.zeros(self.nbins_phi) 605 phi_values = np.zeros(self.nbins_phi) 606 phi_err = np.zeros(self.nbins_phi) 598 607 599 608 # Shift to apply to calculated phi values in order … … 615 624 continue 616 625 # binning 617 i_phi = int(math.floor((self.nbins_phi) * \626 i_phi = int(math.floor((self.nbins_phi) * 618 627 (phi_value + phi_shift) / (2 * Pi))) 619 628 … … 623 632 phi_bins[i_phi] += frac * data[npt] 624 633 625 if err_data ==None or err_data[npt] == 0.0:634 if err_data is None or err_data[npt] == 0.0: 626 635 if data_n < 0: 627 636 data_n = -data_n … … 636 645 phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 637 646 638 idx = (n umpy.isfinite(phi_bins))647 idx = (np.isfinite(phi_bins)) 639 648 640 649 if not idx.any(): 641 650 msg = "Average Error: No points inside ROI to average..." 642 651 raise ValueError, msg 643 # elif len(phi_bins[idx])!= self.nbins_phi:652 # elif len(phi_bins[idx])!= self.nbins_phi: 644 653 # print "resulted",self.nbins_phi- len(phi_bins[idx]) 645 654 #,"empty bin(s) due to tight binning..." … … 748 757 starting from the x- axis on the left-hand side 749 758 """ 759 750 760 def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20): 751 761 self.r_min = r_min … … 769 779 770 780 # Get the all data & info 771 data = data2D.data[n umpy.isfinite(data2D.data)]772 q_data = data2D.q_data[n umpy.isfinite(data2D.data)]773 err_data = data2D.err_data[n umpy.isfinite(data2D.data)]774 qx_data = data2D.qx_data[n umpy.isfinite(data2D.data)]775 qy_data = data2D.qy_data[n umpy.isfinite(data2D.data)]781 data = data2D.data[np.isfinite(data2D.data)] 782 q_data = data2D.q_data[np.isfinite(data2D.data)] 783 err_data = data2D.err_data[np.isfinite(data2D.data)] 784 qx_data = data2D.qx_data[np.isfinite(data2D.data)] 785 qy_data = data2D.qy_data[np.isfinite(data2D.data)] 776 786 dq_data = None 777 787 … … 803 813 dq_overlap_y *= dq_overlap_y 804 814 805 dq_overlap = n umpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)815 dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 806 816 if dq_overlap < 0: 807 817 dq_overlap = y_min 808 dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 809 dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 818 dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 819 dqy_data = data2D.dqy_data[np.isfinite( 820 data2D.data)] - dq_overlap 810 821 # def; dqx_data = dq_r dqy_data = dq_phi 811 822 # Convert dq 2D to 1D here 812 823 dqx = dqx_data * dqx_data 813 824 dqy = dqy_data * dqy_data 814 dq_data = n umpy.add(dqx, dqy)815 dq_data = n umpy.sqrt(dq_data)816 817 # set space for 1d outputs818 x = n umpy.zeros(self.nbins)819 y = n umpy.zeros(self.nbins)820 y_err = n umpy.zeros(self.nbins)821 x_err = n umpy.zeros(self.nbins)822 y_counts = n umpy.zeros(self.nbins)825 dq_data = np.add(dqx, dqy) 826 dq_data = np.sqrt(dq_data) 827 828 # set space for 1d outputs 829 x = np.zeros(self.nbins) 830 y = np.zeros(self.nbins) 831 y_err = np.zeros(self.nbins) 832 x_err = np.zeros(self.nbins) 833 y_counts = np.zeros(self.nbins) 823 834 824 835 # Get the min and max into the region: 0 <= phi < 2Pi … … 839 850 phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 840 851 841 # #No need to calculate the frac when all data are within range852 # No need to calculate the frac when all data are within range 842 853 if self.r_min <= q_value and q_value <= self.r_max: 843 854 frac = 1 844 855 if frac == 0: 845 856 continue 846 # In case of two ROIs (symmetric major and minor regions)(for 'q2')857 # In case of two ROIs (symmetric major and minor regions)(for 'q2') 847 858 if run.lower() == 'q2': 848 # #For minor sector wing859 # For minor sector wing 849 860 # Calculate the minor wing phis 850 861 phi_min_minor = flip_phi(phi_min - Pi) … … 852 863 # Check if phis of the minor ring is within 0 to 2pi 853 864 if phi_min_minor > phi_max_minor: 854 is_in = (phi_value > phi_min_minor or \855 865 is_in = (phi_value > phi_min_minor or 866 phi_value < phi_max_minor) 856 867 else: 857 is_in = (phi_value > phi_min_minor and \858 859 860 # For all cases(i.e.,for 'q', 'q2', and 'phi')861 # Find pixels within ROI868 is_in = (phi_value > phi_min_minor and 869 phi_value < phi_max_minor) 870 871 # For all cases(i.e.,for 'q', 'q2', and 'phi') 872 # Find pixels within ROI 862 873 if phi_min > phi_max: 863 is_in = is_in or (phi_value > phi_min or \864 874 is_in = is_in or (phi_value > phi_min or 875 phi_value < phi_max) 865 876 else: 866 is_in = is_in or (phi_value >= phi_min and \867 877 is_in = is_in or (phi_value >= phi_min and 878 phi_value < phi_max) 868 879 869 880 if not is_in: … … 885 896 i_bin = self.nbins - 1 886 897 887 # #Get the total y898 # Get the total y 888 899 y[i_bin] += frac * data_n 889 900 x[i_bin] += frac * q_value … … 923 934 #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 924 935 x[i] = x[i] / y_counts[i] 925 y_err[y_err == 0] = n umpy.average(y_err)926 idx = (n umpy.isfinite(y) & numpy.isfinite(y_err))936 y_err[y_err == 0] = np.average(y_err) 937 idx = (np.isfinite(y) & np.isfinite(y_err)) 927 938 if x_err != None: 928 939 d_x = x_err[idx] / y_counts[idx] … … 932 943 msg = "Average Error: No points inside sector of ROI to average..." 933 944 raise ValueError, msg 934 # elif len(y[idx])!= self.nbins:945 # elif len(y[idx])!= self.nbins: 935 946 # print "resulted",self.nbins- len(y[idx]), 936 947 #"empty bin(s) due to tight binning..." … … 946 957 The number of bin in phi also has to be defined. 947 958 """ 959 948 960 def __call__(self, data2D): 949 961 """ … … 965 977 The number of bin in Q also has to be defined. 966 978 """ 979 967 980 def __call__(self, data2D): 968 981 """ … … 987 1000 in anti-clockwise starting from the x- axis on the left-hand side 988 1001 """ 1002 989 1003 def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0): 990 1004 # Minimum radius … … 1012 1026 qx_data = data2D.qx_data 1013 1027 qy_data = data2D.qy_data 1014 q_data = n umpy.sqrt(qx_data * qx_data + qy_data * qy_data)1028 q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 1015 1029 1016 1030 # check whether or not the data point is inside ROI … … 1023 1037 Find a rectangular 2D region of interest. 1024 1038 """ 1039 1025 1040 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 1026 1041 # Minimum Qx value [A-1] … … 1077 1092 and (phi_max-phi_min) should not be larger than pi 1078 1093 """ 1094 1079 1095 def __init__(self, phi_min=0, phi_max=math.pi): 1080 1096 self.phi_min = phi_min … … 1113 1129 1114 1130 # get phi from data 1115 phi_data = n umpy.arctan2(qy_data, qx_data)1131 phi_data = np.arctan2(qy_data, qx_data) 1116 1132 1117 1133 # Get the min and max into the region: -pi <= phi < Pi … … 1120 1136 # check for major sector 1121 1137 if phi_min_major > phi_max_major: 1122 out_major = (phi_min_major <= phi_data) + (phi_max_major > phi_data) 1138 out_major = (phi_min_major <= phi_data) + \ 1139 (phi_max_major > phi_data) 1123 1140 else: 1124 out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data) 1141 out_major = (phi_min_major <= phi_data) & ( 1142 phi_max_major > phi_data) 1125 1143 1126 1144 # minor sector … … 1132 1150 if phi_min_minor > phi_max_minor: 1133 1151 out_minor = (phi_min_minor <= phi_data) + \ 1134 1152 (phi_max_minor >= phi_data) 1135 1153 else: 1136 1154 out_minor = (phi_min_minor <= phi_data) & \ 1137 1155 (phi_max_minor >= phi_data) 1138 1156 out = out_major + out_minor 1139 1157
Note: See TracChangeset
for help on using the changeset viewer.