- Timestamp:
- Jul 11, 2016 4:36:55 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.1.1, release-4.1.2, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 7a219e3e
- Parents:
- f2bbabf
- git-author:
- Lewis O'Driscoll <lewis.o'driscoll@…> (07/07/16 11:30:07)
- git-committer:
- Lewis O'Driscoll <lewis.o'driscoll@…> (07/11/16 04:36:55)
- Location:
- src/sas
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/corfunc/corfunc_calculator.py
r77e9ac6 r911dbe4 55 55 56 56 57 def __init__(self, data , lowerq, upperq, background=0, scale=1):57 def __init__(self, data=None, lowerq=None, upperq=None, scale=1): 58 58 """ 59 59 Initialize the class. 60 60 61 61 :param data: Data of the type DataLoader.Data1D 62 :param background: Background value. Will be subtracted from the data 63 before processing 62 :param lowerq: The Q value to use as the boundary for 63 Guinier extrapolation 64 :param upperq: A tuple of the form (lower, upper). 65 Values between lower and upper will be used for Porod extrapolation 64 66 :param scale: Scaling factor for I(q) 65 67 """ 66 print "Before: {}".format(data.y[0]) 67 self._data = self._get_data(data, background, scale) 68 print "After: {}".format(self._data.y[0]) 68 self._data = None 69 self.set_data(data, scale) 69 70 self.lowerq = lowerq 70 71 self.upperq = upperq 72 self.background = 0 73 74 def set_data(self, data, scale=1): 75 """ 76 Prepares the data for analysis 77 78 :return: new_data = data * scale - background 79 """ 80 if data is None: 81 return 82 # Only process data of the class Data1D 83 if not issubclass(data.__class__, Data1D): 84 raise ValueError, "Data must be of the type DataLoader.Data1D" 85 86 # Prepare the data 87 new_data = Data1D(x=data.x, y=data.y) 88 new_data *= scale 89 90 # Ensure the errors are set correctly 91 if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ 92 (min(new_data.dy) == 0 and max(new_data.dy) == 0): 93 new_data.dy = np.ones(len(new_data.x)) 94 95 self._data = new_data 96 97 def compute_background(self, upperq=None): 98 """ 99 Compute the background level from the Porod region of the data 100 """ 101 if self._data is None: return 0 102 elif upperq is None and self.upperq is not None: upperq = self.upperq 103 elif upperq == 0: return 0 104 q = self._data.x 105 mask = np.logical_and(q > upperq[0], q < upperq[1]) 106 _, _, bg = self._fit_porod(q[mask], self._data.y[mask]) 107 108 return bg 71 109 72 110 def compute_extrapolation(self): … … 82 120 return extrapolation 83 121 84 def compute_transform(self, extrapolation ):122 def compute_transform(self, extrapolation, background=None): 85 123 """ 86 124 Transform an extrapolated scattering curve into a correlation function. … … 88 126 qs = extrapolation.x 89 127 iqs = extrapolation.y 128 q = self._data.x 129 if background is None: background = self.background 90 130 91 gamma = dct( iqs*qs**2)131 gamma = dct((iqs-background)*qs**2) 92 132 gamma = gamma / gamma.max() 93 133 xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) … … 97 137 return transform 98 138 99 def _porod(self, q, K, sigma ):139 def _porod(self, q, K, sigma, bg): 100 140 """Equation for the Porod region of the data""" 101 return (K*q**(-4))*np.exp(-q**2*sigma**2)141 return bg + (K*q**(-4))*np.exp(-q**2*sigma**2) 102 142 103 143 def _fit_guinier(self, q, iq): … … 106 146 return lstsq(A, np.log(iq)) 107 147 108 def _get_data(self, data, background, scale): 109 """ 110 Prepares the data for analysis 111 112 :return: new_data = data * scale - background 113 """ 114 # Only process data of the class Data1D 115 if not issubclass(data.__class__, Data1D): 116 raise ValueError, "Data must be of the type DataLoader.Data1D" 117 118 # Prepare the data 119 new_data = (scale * data) 120 new_data.y -= background 121 122 # Check the vector lengths are equal 123 assert len(new_data.x) == len(new_data.y) 124 125 # Ensure the errors are set correctly 126 if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ 127 (min(new_data.dy) == 0 and max(new_data.dy) == 0): 128 new_data.dy = np.ones(len(new_data.x)) 129 130 return new_data 148 def _fit_porod(self, q, iq): 149 """Fit the Porod region of the curve""" 150 fitp = curve_fit(lambda q, k, sig, bg: self._porod(q, k, sig, bg)*q**2, 151 q, iq*q**2)[0] 152 k, sigma, bg = fitp 153 return k, sigma, bg 131 154 132 155 def _fit_data(self, q, iq): … … 137 160 # Returns an array where the 1st and 2nd elements are the values of k 138 161 # and sigma for the best-fit Porod function 139 fitp = curve_fit(lambda q, k, sig: self._porod(q, k, sig)*q**2,140 q[mask], iq[mask]*q[mask]**2)[0]162 k, sigma, bg = self._fit_porod(q[mask], iq[mask]) 163 self.background = bg 141 164 142 165 # Smooths between the best-fit porod function and the data to produce a … … 144 167 data = interp1d(q, iq) 145 168 s1 = self._Interpolator(data, 146 lambda x: self._porod(x, fitp[0], fitp[1]), self.upperq[0], q[-1])169 lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1]) 147 170 148 171 mask = np.logical_and(q < self.lowerq, 0 < q) -
src/sas/sasgui/perspectives/corfunc/corfunc.py
r3b8efec r911dbe4 22 22 IQ_DATA_LABEL = r"$I_{obs}(q)$" 23 23 IQ_EXTRAPOLATED_DATA_LABEL = r"$I_{extrap}(q)$" 24 25 GROUP_ID_TRANSFORM = r"$\Gamma(x)$" 26 TRANSFORM_LABEL = r"$\Gamma(x)$" 24 27 25 28 … … 129 132 new_plot = Data1D(data.x, data.y, dy=data.dy) 130 133 134 group_id = "" 131 135 if label == IQ_DATA_LABEL or label == IQ_EXTRAPOLATED_DATA_LABEL: 132 136 new_plot.xaxis("\\rm{Q}", 'A^{-1}') 133 137 new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}") 134 138 group_id = GROUP_ID_IQ_DATA 139 elif label == TRANSFORM_LABEL: 140 new_plot.xaxis("{x}", 'A') 141 new_plot.yaxis("{\\Gamma}", '') 142 group_id = GROUP_ID_TRANSFORM 135 143 new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM 136 144 new_plot.id = label 137 145 new_plot.name = label 146 new_plot.group_id = group_id 138 147 new_plot.interactive = True 139 new_plot.group_id = GROUP_ID_IQ_DATA 140 new_plot.title = "I(q)" 148 new_plot.title = group_id.replace('$', '').replace('\\', '') 141 149 # Show data on a linear scale 142 150 new_plot.xtransform = 'x' 143 151 new_plot.ytransform = 'y' 144 152 wx.PostEvent(self.parent, 145 NewPlotEvent(plot=new_plot, title="I(q)", reset=reset)) 153 NewPlotEvent(plot=new_plot, title=new_plot.title, 154 reset=reset)) -
src/sas/sasgui/perspectives/corfunc/corfunc_panel.py
rf2bbabf r911dbe4 38 38 self.SetWindowVariant(variant=FONT_VARIANT) 39 39 self._manager = manager 40 self._data = data # The data to be analysed 40 # The data with no correction for background values 41 self._data = data # The data to be analysed (corrected fr background) 41 42 self._extrapolated_data = None # The extrapolated data set 43 self._calculator = CorfuncCalculator() 42 44 self._data_name_box = None # Text box to show name of file 43 45 self._background_input = None … … 45 47 self._qmax1_input = None 46 48 self._qmax2_input = None 49 self._transform_btn = None 50 self._compute_btn = None 47 51 self.qmin = 0 48 52 self.qmax = (0, 0) … … 107 111 :param data: The data that has been loaded 108 112 """ 113 self._transform_btn.Disable() 109 114 if data is None: 110 115 return 111 116 self._data_name_box.SetValue(str(data.title)) 112 117 self._data = data 118 self._calculator.set_data(data) 113 119 if self._manager is not None: 114 120 from sas.sasgui.perspectives.corfunc.corfunc import IQ_DATA_LABEL 115 self._manager.show_data( data, IQ_DATA_LABEL, reset=True)121 self._manager.show_data(self._data, IQ_DATA_LABEL, reset=True) 116 122 if set_qrange: 117 123 lower = data.x[-1]*0.05 … … 120 126 self.set_qmin(lower) 121 127 self.set_qmax((upper1, upper2)) 122 self.set_background( 0.0)128 self.set_background(self._calculator.compute_background(self.qmax)) 123 129 124 130 def get_data(self): … … 134 140 wx.PostEvent(self.parent.parent, StatusEvent(status=msg)) 135 141 return 136 calculator = CorfuncCalculator(self._data, self.qmin, self.qmax, 137 background=self.background) 138 self._extrapolated_data = calculator.compute_extrapolation() 142 self._calculator.set_data(self._data) 143 self._calculator.lowerq = self.qmin 144 self._calculator.upperq = self.qmax 145 self._extrapolated_data = self._calculator.compute_extrapolation() 139 146 # TODO: Find way to set xlim and ylim so full range of data can be 140 147 # plotted … … 148 155 IQ_EXTRAPOLATED_DATA_LABEL 149 156 self._manager.show_data(to_plot, IQ_EXTRAPOLATED_DATA_LABEL) 157 self._transform_btn.Enable() 158 159 def compute_transform(self, event=None): 160 transformed_data = self._calculator.compute_transform( 161 self._extrapolated_data, self.background) 162 from sas.sasgui.perspectives.corfunc.corfunc import TRANSFORM_LABEL 163 import numpy as np 164 plot_x = transformed_data.x[np.where(transformed_data.x <= 200)] 165 plot_y = transformed_data.y[np.where(transformed_data.x <= 200)] 166 self._manager.show_data(Data1D(plot_x, plot_y), TRANSFORM_LABEL) 150 167 151 168 … … 183 200 self._background_input.SetValue(str(bg)) 184 201 202 def _compute_background(self, event=None): 203 self.set_background(self._calculator.compute_background(self.qmax)) 185 204 186 205 def _on_enter_input(self, event=None): … … 200 219 group_id = GROUP_ID_IQ_DATA 201 220 if event is not None: 221 if self._manager is not None and\ 222 event.GetEventObject() == self._background_input: 223 from sas.sasgui.perspectives.corfunc.corfunc\ 224 import IQ_DATA_LABEL 225 self._manager.show_data(self._data, IQ_DATA_LABEL, reset=True) 202 226 wx.PostEvent(self._manager.parent, PlotQrangeEvent( 203 227 ctrl=[self._qmin_input, self._qmax1_input, self._qmax2_input], … … 309 333 wx.RIGHT | wx.EXPAND, 5) 310 334 335 background_button = wx.Button(self, wx.NewId(), "Calculate", 336 size=(75, 25)) 337 background_button.Bind(wx.EVT_BUTTON, self._compute_background) 338 q_sizer.Add(background_button, (1, 2), (1, 1), wx.RIGHT | wx.EXPAND, 5) 339 311 340 qrange_label = wx.StaticText(self, -1, "Q Range:", size=(50,20)) 312 341 q_sizer.Add(qrange_label, (2,0), (1,1), wx.LEFT | wx.EXPAND, 5) … … 395 424 396 425 extrapolate_btn = wx.Button(self, wx.NewId(), "Extrapolate") 397 transform_btn = wx.Button(self, wx.NewId(), "Transform") 398 compute_btn = wx.Button(self, wx.NewId(), "Compute Measuments") 426 self._transform_btn = wx.Button(self, wx.NewId(), "Transform") 427 self._compute_btn = wx.Button(self, wx.NewId(), "Compute Measuments") 428 429 self._transform_btn.Disable() 430 self._compute_btn.Disable() 399 431 400 432 extrapolate_btn.Bind(wx.EVT_BUTTON, self.compute_extrapolation) 433 self._transform_btn.Bind(wx.EVT_BUTTON, self.compute_transform) 401 434 402 435 controls_sizer.Add(extrapolate_btn, wx.CENTER | wx.EXPAND) 403 controls_sizer.Add( transform_btn, wx.CENTER | wx.EXPAND)404 controls_sizer.Add( compute_btn, wx.CENTER | wx.EXPAND)436 controls_sizer.Add(self._transform_btn, wx.CENTER | wx.EXPAND) 437 controls_sizer.Add(self._compute_btn, wx.CENTER | wx.EXPAND) 405 438 406 439 controlbox_sizer.Add(controls_sizer, wx.TOP | wx.EXPAND, 0)
Note: See TracChangeset
for help on using the changeset viewer.