From 6689e94e78a9d3378d24e02d621eedec1c0fcb6b Mon Sep 17 00:00:00 2001 From: "Purnima P. Balakrishnan" Date: Fri, 18 Feb 2022 17:07:00 -0500 Subject: [PATCH 1/3] Fitting 'sum' and difference, plus error in variables Added a probe which contains sum and difference cross-sections, an experiment which calculates residuals based on these rather than pp and mm. Save from gui implemented. Also includes error-in-variables. --- .gitignore | 1 + refl1d/experiment.py | 213 ++++++++++++++++++++++++++++++ refl1d/names.py | 4 +- refl1d/probe.py | 299 ++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 511 insertions(+), 6 deletions(-) diff --git a/.gitignore b/.gitignore index 7ff757439..c706d4a1f 100644 --- a/.gitignore +++ b/.gitignore @@ -23,3 +23,4 @@ /bumps-* /Refl1D.pdf /.eggs +*__pycache__* diff --git a/refl1d/experiment.py b/refl1d/experiment.py index f9b52db68..1557b4e2a 100644 --- a/refl1d/experiment.py +++ b/refl1d/experiment.py @@ -830,3 +830,216 @@ def nice(v, digits=2): place = floor(log10(abs(v))) scale = 10**(place-(digits-1)) return sign*floor(abs(v)/scale+0.5)*scale + +#EIV marginalized_residuals from Paul +def marginalized_residuals(Q, FQ, R, dR, angle_uncertainty=0.002, wavelength=None): + r""" + Returns the residuals from an error-in-variables model marginalized over + the variables. + + *angular_uncertainty* from motor jitter in degrees. + + For error in variables fits with normal uncertainty, start with the + following model: + + .. math:: + + x &=& x_o + \epsilon_1 \text{for} \epsilon_1 ~ N(0, \Delta x^2) \\ + y &=& f(x) + \epsilon_2 \text{for} \epsilon_2 ~ N(0, \Delta y^2) \\ + + Use a linear approximation at the nominal measurement location $x_o$ then + + .. math:: + + f(x) \approx f'(x_?)(x - x_o) + f(x_o) + + and so + + .. math:: + + y &\approx& f'(x_o)(x_o + \epsilon_1 - x_o) + f(x_o) + \epsilon_2 \\ + &=& f(x_o) + f'(x_o)\epsilon_1 + \epsilon_2 + + Therefore, measured $y_o$ is distributed as + + .. math:: + + y_o &~& f(x_o) + f'(x_o)N(0, \Delta x^2) + N(0, \Delta y^2) \\ + &~& N(f(x_o), (f'(x_o)\Delta x)^2 + \Delta y^2) + + That is, assuming that f(x) is approximately linear over $\Delta x$, then + simply add $f'(x_o)\Delta x^2$ to the variance in the data. + + We should be sampling $x$ densely enough that we can approximate + $f'(x)$ using the center point formula, + + .. math:: + + f'(x) \approx \frac{f(x_{k+1}) - f(x_{k-1})}{x_{k+1} - x_{k-1}} + + For reflectometry specifically the motor uncertainty $\delta\theta$ is + in angle and the theory is in $Q$, so we use the following + + .. math:: + + Q &=& \frac{4 \pi}{\lambda} \sin(\theta + \delta\theta) \\ + &=& \frac{4 \pi}{\lambda} (\cos \delta\theta \sin \theta + + \sin \delta\theta \cos \theta) + + Using the small angle approximation and $\cos \theta > 0.96$ + for $\theta < 15^o$, then + + .. math:: + + Q &\approx& \frac{4 \pi}{\lambda} (\sin \theta + \delta\theta\cos\theta) + &\approx& \frac{4 \pi}{\lambda} (\sin \theta + \delta\theta) + &=$ Q + \frac{4 \pi}{\lambda}\delta\theta + + and so + + .. math:: + + $\epsilon_1 = \delta Q = \frac{4 \pi}{\lambda}\delta\theta + + This means that we can compute the residual using + + .. math:: + + \frac{dR}{dQ} &\approx& \frac{R_{k+1} - R_{k-1}}{Q_{k+1} - Q_{k-1}} \\ + \Delta R' &=& \sqrt{\Delta R^2 + + \left(\frac{4\pi\delta\theta}{\lambda} \frac{dR}{dQ}\right)^2 } + + You can arrive at the same expression using marginalization over the + possible incident angles $\theta + \delta\theta$ for each $Q$ point + + .. math:: + + P(R | Q) = \int P(R | Q, \theta) P(\theta)\,\mathrm{d}\theta + + where $P(R | Q, \theta)$ is the usual Gaussian residual for the measurement + and $P(\theta)$ is a Gaussian uncertainty in angle. + """ + + # slope from center point formula + if angle_uncertainty == 0.0: + return (R - FQ)/dR + # Using small angle approximation to Q = 4 pi/L sin(T + d) + # Q = 4 pi / L (cos d sin T + sin d cos T) + # ~ 4 pi / L (sin T + d cos T) since d is small + # ~ 4 pi / L (sin T + d) since cos T > 0.96 for T < 15 degrees + # = Q + 4 pi / L d + # ~ Q + 2.5 d since L in [4, 6] angstroms + DQ = 4*np.pi/wavelength*np.radians(angle_uncertainty) + + # Quick approx to [ log integral P(R,dR;Q') P(Q') dQ'] for motor position + # uncertainty P(Q') and gaussian measurement uncertainty P(R;Q') is to + # increase the size of dR based on the slope at Q and the motor uncertainty. + dRdQ = (R[2:] - R[:-2])/(Q[2:] - Q[:-2]) + # Leave dR untouched for the initial and final point + dRp = dR.copy() + dRp[1:-1] = np.sqrt((dR[1:-1])**2 + (DQ[1:-1]*dRdQ)**2) # add in quadrature + return (R - FQ)/(dRp) + +# TODO: check curvature in marginalized_residuals() +# +# If $|f(x_k) - \hat f(x_k)| \gtrsim |\delta x f'(x_k)|$ where $\hat f(x)$ is +# the line connecting $f(x_k-\delta x)$ to $f(x_k+\delta x)$ then the curvature +# at $x_k$ is too large for the correction factor. Depending on whether it is +# curving toward or away from the measured data point the scaled residuals +# will be too large or too small. +# +# Better: +# Check $(f(x_k) - \hat f(x_k))^2 > C (\delta x f'(x_k))^2 + \delta y^2$. +# That way we ignore the curvature if it is less than some fraction $C$ of +# the combined uncertainty in x and y together. +# +# What should we do when the curvature check fails? Maybe warn the user, +# or maybe further tweak the mean and variance used to compute the residuals. +# +# Perhaps an additional term in the Taylor series around $f(x_k)$ will +# do the trick, adding $f''(x_0) \epsilon_1^2 / 2$ to the approximation of $y$. +# That is, $Z = f''/2 X + f' X + f + Y$. Completing the square, this is +# $Z = f''/2 (X + f'/f'')^2 + f - (f'/f'')^2 + Y$. The X expression +# corresponds to a non-central $\chi'^2_1$ distribution.[1] Using the +# gaussian approximation (i.e., the mean and variance) for this distribution[2], +# we can then find $Z = X' + Y$, and combine them into a single gaussian +# approximation as before. Or use the generalized $\chi^2$ distribution, +# though that may be more difficult to work with numerically.[3] +# +# [1] https://stats.stackexchange.com/questions/127612/polynomial-transform-of-a-gaussian-random-variable#comment244825_127612 +# [2] https://en.wikipedia.org/wiki/Noncentral_chi-squared_distribution +# [3] https://en.wikipedia.org/wiki/Generalized_chi-squared_distribution + + +class SumDiffExperiment(Experiment): + + def residuals(self): + if 'residuals' in self._cache: + return self._cache['residuals'] + + if self.probe.polarized: + have_data = not all(x is None or x.R is None for x in self.probe.xs) + else: + have_data = not (self.probe.R is None) + if not have_data: + resid = np.zeros(0) + self._cache['residuals'] = resid + return resid + + QR = self.reflectivity() + # print('in residuals') + # print(len(QR)) + # for q in QR: + # if q is not None: + # print(len(q[0])) + # else: + # print('None') + # print(len(QR[0][0])) + # print(len(QR[3][0])) + # QRmean = meanreflectivity(self.probe.pp.Q,QR[0][0],None,self.probe.pp.Q,QR[3][0],None) + # print('from data') + # print(len(measRmean[1])) + mm, mp, pm, pp = QR + QRmean = meanreflectivity(pp[0],pp[1],None,mm[0],mm[1],None) + # print(len(QRmean[1])) + QRdf = splitting(pp[0],pp[1],None,mm[0],mm[1],None) + # print(len(QRsa[1])) + # print(len(self.probe.sa.Q)) +# print(self.probe.df.Q) + resid = np.hstack([(xs[1] - QRi[1])/xs[2] + for xs, QRi in zip([(self.probe.sm.Q,self.probe.sm.R,self.probe.sm.dR),(self.probe.df.Q,self.probe.df.R,self.probe.df.dR)], [QRmean,QRdf]) + if xs is not None]) + # print(resid) + self._cache['residuals'] = resid + return resid + + +class SumDiffEIVExperiment(Experiment): + + def residuals(self): + if 'residuals' in self._cache: + return self._cache['residuals'] + + if self.probe.polarized: + have_data = not all([x is None or x.R is None for x in self.probe.xs]) +# print(have_data) + else: + have_data = not (self.probe.R is None) + if not have_data: + resid = np.zeros(0) + self._cache['residuals'] = resid + return resid + + QR = self.reflectivity() + + mm, mp, pm, pp = QR + + QRmean = meanreflectivity(pp[0],pp[1],None,mm[0],mm[1],None) + QRdf = splitting(pp[0],pp[1],None,mm[0],mm[1],None) + + + resid = np.hstack([marginalized_residuals(QRi[0], QRi[1], xs.R, xs.dR, angle_uncertainty=getattr(xs, 'angle_uncertainty', 0.0), wavelength=xs.L) + for xs, QRi in zip([self.probe.sm,self.probe.df], [QRmean,QRdf]) if xs is not None]) + # print(resid) + self._cache['residuals'] = resid + return resid diff --git a/refl1d/names.py b/refl1d/names.py index 133658fdc..2c3dde1c6 100644 --- a/refl1d/names.py +++ b/refl1d/names.py @@ -20,7 +20,7 @@ from bumps.fitproblem import FitProblem from bumps.fitproblem import MultiFitProblem # deprecated -from .experiment import Experiment, plot_sample, MixedExperiment +from .experiment import Experiment, plot_sample, MixedExperiment, SumDiffExperiment, SumDiffEIVExperiment from .flayer import FunctionalProfile, FunctionalMagnetism from .material import SLD, Material, Compound, Mixture from .model import Slab, Stack @@ -30,7 +30,7 @@ from .cheby import FreeformCheby, ChebyVF, cheby_approx, cheby_points from .interface import Erf from .probe import (Probe, ProbeSet, XrayProbe, NeutronProbe, QProbe, - PolarizedNeutronProbe, PolarizedQProbe, load4) + PolarizedNeutronProbe, PolarizedNeutronProbeSumDiff, PolarizedQProbe, load4) from .stajconvert import load_mlayer, save_mlayer from . import ncnrdata as NCNR, snsdata as SNS from .instrument import Monochromatic, Pulsed diff --git a/refl1d/probe.py b/refl1d/probe.py index ada68668a..def9abf85 100644 --- a/refl1d/probe.py +++ b/refl1d/probe.py @@ -578,7 +578,7 @@ def oversample(self, n=20, seed=1): Note: :meth:`oversample` will remove the extra Q calculation points introduced by :meth:`critical_edge`. """ - + rng = numpy.random.RandomState(seed=seed) T = rng.normal(self.T[:, None], self.dT[:, None], size=(len(self.dT), n-1)) L = rng.normal(self.L[:, None], self.dL[:, None], size=(len(self.dL), n-1)) @@ -1607,7 +1607,7 @@ class PolarizedNeutronProbe(object): """ Polarized neutron probe - *xs* (4 x NeutronProbe) is a sequence pp, pm, mp and mm. + *xs* (4 x NeutronProbe) is a sequence mm, mp, pm, and pp. *Aguide* (degrees) is the angle of the applied field relative to the plane of the sample, with angle 270º in the plane of the sample. @@ -1767,7 +1767,7 @@ def _calculate_union(self): self._set_calc(self.T, self.L) else: self._oversample(self.oversampling, self.oversampling_seed) - + self._theta_offsets = theta_offsets @property @@ -1966,7 +1966,7 @@ def spin_asymmetry(Qp, Rp, dRp, Qm, Rm, dRm): .. math:: - \Delta S_A^2 = \frac{4(R_{++}^2\Delta R_{--}^2+R_{--}^2\Delta R_{++})} + \Delta S_A^2 = \frac{4(R_{++}^2\Delta R_{--}^2+R_{--}^2\Delta R_{++}^2)} {(R_{++} + R_{--})^4} """ @@ -1980,7 +1980,298 @@ def spin_asymmetry(Qp, Rp, dRp, Qm, Rm, dRm): else: return Qp, v, None +# based on A.J.Caruana Footprint correction probe / experiment +class PolarizedNeutronProbeSumDiff(PolarizedNeutronProbe): + """ + Polarized neutron probe WITH Sum and Difference + + *xs* (6 x NeutronProbe) is a sequence mm, mp, pm, pp, sum, difference + + *Aguide* (degrees) is the angle of the applied field relative + to the plane of the sample, with angle 270º in the plane of the sample. + + *H* (tesla) is the magnitude of the applied field + """ + view = None # Default to Probe.view when None + show_resolution = None # Default to Probe.show_resolution when None + substrate = surface = None + polarized = True + radiation = 'neutron' + + @property + def sm(self): + return self.xs[4] + + @property + def df(self): + return self.xs[5] + + @property + def xs(self): + return self._xs # Don't let user replace xs + + def parameters(self): + mm, mp, pm, pp, sm, df = [(xsi.parameters() if xsi else None) + for xsi in self.xs] + return { + 'pp': pp, 'pm': pm, 'mp': mp, 'mm': mm, 'sm': sm, 'df': df, + 'Aguide': self.Aguide, 'H': self.H, + } + + + def to_dict(self): + """ Return a dictionary representation of the parameters """ + mm, mp, pm, pp, sm, df = self.xs + return to_dict({ + 'type': type(self).__name__, + 'name': self.name, + 'pp': pp, + 'pm': pm, + 'mp': mp, + 'mm': mm, + 'sm': sm, + 'df': df, + 'a_guide': self.Aguide, + 'h': self.H, + }) + + + + def save(self, filename, theory, substrate=None, surface=None): + # for xsi, xsi_th, suffix in zip(self.xs, theory, ('A', 'B', 'C', 'D', 'SM', 'DF')): + xsi_df = self.df + mm, mp, pm, pp = theory + Q, DF, _ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + xsi_th_df = [Q,DF] + suffix = 'DF' + if xsi_df is not None: + xsi_df.save(filename+suffix, xsi_th_df, + substrate=substrate, surface=surface) + xsi_sm = self.sm + mm, mp, pm, pp = theory + Q, SM, _ = meanreflectivity(pp[0], pp[1], None, mm[0], mm[1], None) + xsi_th_sm = [Q,SM] + suffix = 'SM' + if xsi_sm is not None: + xsi_sm.save(filename+suffix, xsi_th_sm, + substrate=substrate, surface=surface) + +# save.__doc__ = PolarizedNeutronProbe.save.__doc__ + + + def plot_resolution(self, **kwargs): + self._xs_plot('plot_resolution', **kwargs) + + + def plot_linear(self, **kwargs): + self._xs_plot('plot_linear', **kwargs) + + + def plot_log(self, **kwargs): + self._xs_plot('plot_log', **kwargs) + + + def plot_fresnel(self, **kwargs): + self._xs_plot('plot_fresnel', **kwargs) + + + def plot_logfresnel(self, **kwargs): + self._xs_plot('plot_logfresnel', **kwargs) + + + def plot_Q4(self, **kwargs): + self._xs_plot('plot_Q4', **kwargs) + + #plot SA residuals NOT for other xs + def plot_residuals(self, **kwargs): + self._xs_plot('plot_residuals', **kwargs) + + + def plot_SA(self, theory=None, label=None, plot_shift=None, + **kwargs): + import matplotlib.pyplot as plt + # print('plotting Difference') + if self.df is None:# and (self.pp is None or self.mm is None): + # raise TypeError("cannot form difference plot without difference") + raise TypeError("cannot form difference plot without difference") + + plot_shift = plot_shift if plot_shift is not None else Probe.plot_shift + trans = auto_shift(plot_shift) + df = self.df + c = coordinated_colors() + if hasattr(df, 'R') and df.R is not None: #plot measured DF (not from pp, mm) + Q = df.Q + DF = df.R + dDF = df.dR + if dDF is not None: + res = (self.show_resolution if self.show_resolution is not None + else Probe.show_resolution) + dQ = df.dQ if res else None + plt.errorbar(10*Q, ((10*Q)**4)*DF, yerr=((10*Q)**4)*dDF, xerr=10*dQ, fmt='.', capsize=0, + label=df.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + else: + plt.plot(10*Q, ((10*Q)**4)*DF, '.', + label=df.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + if theory is not None: #plot theory curve + mm, mp, pm, pp = theory + Q, DF, _ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + plt.plot(10*Q, ((10*Q)**4)*DF, + label=self.pp.label(prefix=label, gloss='theory'), + transform=trans, + color=c['dark']) + plt.xlabel(r'Q (nm$^{-1}$)') + plt.ylabel(r'Difference ($R^{++} -\, R^{--}$) $\times$ Q$^4$ nm$^{-4}$') + plt.legend(numpoints=1) + + def _xs_plot(self, plotter, theory=None, **kwargs): + import matplotlib.pyplot as plt + # plot DF residuals rather than other xs residuals + # print(plotter) + if plotter == 'plot_residuals': + # plot_shift = plot_shift if plot_shift is not None else Probe.residuals_shift + plot_shift = Probe.residuals_shift + trans = auto_shift(plot_shift) + if theory is not None and self.sm.R is not None: + c = coordinated_colors() + mm, mp, pm, pp = theory + Q, SM, _ = meanreflectivity(pp[0], pp[1], None, mm[0], mm[1], None) + # In case theory curve is evaluated at more/different points... + R = np.interp(self.sm.Q, Q, SM) + residual = (SM - self.sm.R)/self.sm.dR + plt.plot(self.sm.Q, residual, + '.', color=c['light'], + transform=trans, + label=self.sm.label(prefix='', suffix='SM', gloss='resid')) + if theory is not None and self.df.R is not None: + c = coordinated_colors() + mm, mp, pm, pp = theory + Q, DF, _ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + # In case theory curve is evaluated at more/different points... + R = np.interp(self.df.Q, Q, DF) + residual = (DF - self.df.R)/self.df.dR + plt.plot(self.df.Q, residual, + '.', color=c['light'], + transform=trans, + label=self.df.label(prefix='', suffix='DF', gloss='resid')) + plt.axhline(1, color='black', ls='--', lw=1) + plt.axhline(0, color='black', lw=1) + plt.axhline(-1, color='black', ls='--', lw=1) + plt.xlabel('Q (inv A)') + plt.ylabel('(theory-data)/error') + plt.legend(numpoints=1) + #don't plot other cross sections which are placeholders + # else: #don't plot DF on normal plot + # if theory is None: + # theory = (None, None, None, None) + # #only first four cross-sections + # for x_data, x_th, suffix in zip(self.xs[:4], theory, + # ('$^{--}$', '$^{-+}$', '$^{+-}$', '$^{++}$')): + # if x_data is not None: + # fn = getattr(x_data, plotter) + # fn(theory=x_th, suffix=suffix, **kwargs) + + else: # plot NR not other xs + if theory is None: + theory = (None, None, None, None) + #only first four cross-sections + x_data = self.sm + mm, mp, pm, pp = theory + Q,SM,_ = meanreflectivity(pp[0], pp[1], None, mm[0], mm[1], None) + suffix = 'SM' + if x_data is not None: + fn = getattr(x_data, plotter) + fn(theory=[Q,SM], suffix=suffix, **kwargs) + +def meanreflectivity(Qp, Rp, dRp, Qm, Rm, dRm): + r""" + Compute average reflectivity for R++, R--. + + **Parameters:** + + *Qp*, *Rp*, *dRp* : vector + Measured ++ cross section and uncertainty. + *Qm*, *Rm*, *dRm* : vector + Measured -- cross section and uncertainty. + + If *dRp*, *dRm* are None then the returned uncertainty will also be None. + + **Returns:** + + *Q*, *R*, *dR* : vector + Computed reflectivity and uncertainty. + + **Algorithm:** + + Mean Reflectivity, $R$, is: + + .. math:: + + S_A = \frac{R_{++} + R_{--}}{2} + + Uncertainty $\Delta R$ follows from propagation of error: + + .. math:: + + \Delta S_A^2 = \frac{\Delta R_{--}^2+\Delta R_{++}^2)}{4} + + """ + Rm = np.interp(Qp, Qm, Rm) + v = (Rp+Rm)/2 + if dRp is not None: + dRm = np.interp(Qp, Qm, dRm) + dvsq = (dRm**2 + dRp**2) / 4 + dvsq[dvsq < 0] = 0 + return Qp, v, np.sqrt(dvsq) + else: + return Qp, v, None + + +def splitting(Qp, Rp, dRp, Qm, Rm, dRm): + r""" + Compute splitting for R++, R--. + + **Parameters:** + + *Qp*, *Rp*, *dRp* : vector + Measured ++ cross section and uncertainty. + *Qm*, *Rm*, *dRm* : vector + Measured -- cross section and uncertainty. + + If *dRp*, *dRm* are None then the returned uncertainty will also be None. + + **Returns:** + + *Q*, *D*, *dD* : vector + Computed difference and uncertainty. + + **Algorithm:** + + Difference, $D$, is: + .. math:: + + D = R_{++} - R_{--} + + Uncertainty $\Delta S_A$ follows from propagation of error: + + .. math:: + + \Delta D^2 = \Delta R_{--}^2+\Delta R_{++}^2 + + """ + Rm = np.interp(Qp, Qm, Rm) + v = (Rp-Rm) + if dRp is not None: + dRm = np.interp(Qp, Qm, dRm) + dvsq = (dRm)**2 + (dRp)**2 + dvsq[dvsq < 0] = 0 + return Qp, v, sqrt(dvsq) + else: + return Qp, v, None def _interpolate_Q(Q, dQ, n): """ From da72adeb389edaa5709eae3083557b8052008412 Mon Sep 17 00:00:00 2001 From: "Purnima P. Balakrishnan" Date: Tue, 21 Feb 2023 19:12:33 -0500 Subject: [PATCH 2/3] Fix calculation of degrees of freedom Don't count dummy variables used for calculating the correct Q as data points when calculating chi2 --- refl1d/experiment.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/refl1d/experiment.py b/refl1d/experiment.py index 1557b4e2a..f0b97bd4c 100644 --- a/refl1d/experiment.py +++ b/refl1d/experiment.py @@ -25,6 +25,7 @@ from .reflectivity import reflectivity_amplitude as reflamp from .reflectivity import magnetic_amplitude as reflmag from .reflectivity import BASE_GUIDE_ANGLE as DEFAULT_THETA_M +from .probe import PolarizedNeutronProbeSumDiff, meanreflectivity, splitting #print("Using pure python reflectivity calculator") #from .abeles import refl as reflamp from .util import asbytes @@ -1013,6 +1014,15 @@ def residuals(self): self._cache['residuals'] = resid return resid + #probes contain dummy variables for calculating the correct Q - do not want to include + def numpoints(self): + if isinstance(self.probe, PolarizedNeutronProbeSumDiff): + return sum(len(xs.Q) for xs in (self.probe.sm, self.probe.df) if xs is not None) + elif self.probe.polarized: + return sum(len(xs.Q) for xs in self.probe.xs if xs is not None) + else: + return len(self.probe.Q) if self.probe.Q is not None else 0 + class SumDiffEIVExperiment(Experiment): @@ -1042,4 +1052,14 @@ def residuals(self): for xs, QRi in zip([self.probe.sm,self.probe.df], [QRmean,QRdf]) if xs is not None]) # print(resid) self._cache['residuals'] = resid + print("resid", np.sum(resid**2)/2) return resid + + #probes contain dummy variables for calculating the correct Q - do not want to include + def numpoints(self): + if isinstance(self.probe, PolarizedNeutronProbeSumDiff): + return sum(len(xs.Q) for xs in (self.probe.sm, self.probe.df) if xs is not None) + elif self.probe.polarized: + return sum(len(xs.Q) for xs in self.probe.xs if xs is not None) + else: + return len(self.probe.Q) if self.probe.Q is not None else 0 From 606a8c5846bee42d09187c86b1e7c8017cbbe8f8 Mon Sep 17 00:00:00 2001 From: "Purnima P. Balakrishnan" Date: Fri, 24 Feb 2023 10:23:34 -0500 Subject: [PATCH 3/3] Adding additional Difference plot views Added new views to plot difference, Fresnel difference, and difference*Q^4, and correctly calculate spin asymmetry for sum/difference probe --- refl1d/probe.py | 239 +++++++++++++++++++++++++++++++++++---- refl1d/view/data_view.py | 43 ++++++- 2 files changed, 257 insertions(+), 25 deletions(-) diff --git a/refl1d/probe.py b/refl1d/probe.py index def9abf85..2dfdf5a15 100644 --- a/refl1d/probe.py +++ b/refl1d/probe.py @@ -743,8 +743,15 @@ def plot(self, view=None, **kwargs): self.plot_residuals(**kwargs) elif view == 'fft': self.plot_fft(**kwargs) - elif view == 'SA': # SA uses default plot - self.plot(view=None, **kwargs) + # SA and difference (undefined if not polarized) do not plot since they don't exist + elif view == 'SA': + pass + elif view == 'diff': + pass + elif view == 'fresnelDiff': + pass + elif view == 'diffq4': + pass else: raise TypeError("incorrect reflectivity view '%s'"%view) @@ -763,7 +770,8 @@ def plot_linear(self, **kwargs): Plot the data associated with probe. """ import matplotlib.pyplot as plt - self._plot_pair(ylabel='Reflectivity', **kwargs) + kwargs.setdefault('ylabel', 'Reflectivity') + self._plot_pair(**kwargs) plt.yscale('linear') def plot_log(self, **kwargs): @@ -771,7 +779,8 @@ def plot_log(self, **kwargs): Plot the data associated with probe. """ import matplotlib.pyplot as plt - self._plot_pair(ylabel='Reflectivity', **kwargs) + kwargs.setdefault('ylabel','Reflectivity') + self._plot_pair(**kwargs) plt.yscale('log') def plot_logfresnel(self, *args, **kw): @@ -864,7 +873,7 @@ def _plot_pair(self, theory=None, # reflectivity. # Issues with current implementation: # * If the resolution is too tight, there won't be sufficient - # support from _calc_Q to compute theory at Q interpolated. + # support from _cal c_Q to compute theory at Q interpolated. # * dQ for the interpolated points uses linear interpolation # of dQ between neighbours. If measurements with tight and # loose resolution are interleaved the result will look very @@ -1859,6 +1868,12 @@ def plot(self, view=None, **kwargs): self.plot_residuals(**kwargs) elif view == 'SA': self.plot_SA(**kwargs) + elif view == 'diff': + self.plot_Diff(**kwargs) + elif view == 'fresnelDiff': + self.plot_FresnelDiff(**kwargs) + elif view == 'diffq4': + self.plot_DiffQ4(**kwargs) elif view == 'resolution': self.plot_resolution(**kwargs) else: @@ -1925,6 +1940,125 @@ def plot_SA(self, theory=None, label=None, plot_shift=None, plt.ylabel(r'spin asymmetry $(R^{++} -\, R^{--}) / (R^{++} +\, R^{--})$') plt.legend(numpoints=1) + def plot_Diff(self, theory=None, label=None, plot_shift=None, + **kwargs): + import matplotlib.pyplot as plt + if self.pp is None or self.mm is None: + raise TypeError("cannot form difference plot without ++ and --") + + plot_shift = plot_shift if plot_shift is not None else Probe.plot_shift + trans = auto_shift(plot_shift) + pp, mm = self.pp, self.mm + c = coordinated_colors() + if hasattr(pp, 'R') and hasattr(mm, 'R') and pp.R is not None and mm.R is not None: + Q, diff, ddiff = splitting(pp.Q, pp.R, pp.dR, mm.Q, mm.R, mm.dR) + if ddiff is not None: + res = (self.show_resolution if self.show_resolution is not None + else Probe.show_resolution) + dQ = pp.dQ if res else None + plt.errorbar(pp.Q, diff, yerr=ddiff, xerr=dQ, fmt='.', capsize=0, + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + else: + plt.plot(pp.Q, diff, '.', + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + if theory is not None: + mm, mp, pm, pp = theory + Q, diff, ddiff = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + plt.plot(Q, diff, + label=self.pp.label(prefix=label, gloss='theory'), + transform=trans, + color=c['dark']) + plt.xlabel(r'Q (\AA^{-1})') + plt.ylabel(r'Difference $(R^{++} -\, R^{--})$') + plt.legend(numpoints=1) + + def plot_DiffQ4(self, theory=None, label=None, plot_shift=None, + **kwargs): + import matplotlib.pyplot as plt + if self.pp is None or self.mm is None: + raise TypeError("cannot form difference plot without ++ and --") + plot_shift = plot_shift if plot_shift is not None else Probe.plot_shift + trans = auto_shift(plot_shift) + pp, mm = self.pp, self.mm + c = coordinated_colors() + Q4 = 1e-8*self.pp.Q**-4*self.pp.intensity.value + self.pp.background.value + if hasattr(pp, 'R') and hasattr(mm, 'R') and pp.R is not None and mm.R is not None: + Q, diff, ddiff = splitting(pp.Q, pp.R, pp.dR, mm.Q, mm.R, mm.dR) + if ddiff is not None: + res = (self.show_resolution if self.show_resolution is not None + else Probe.show_resolution) + dQ = pp.dQ if res else None + plt.errorbar(Q, diff/Q4, yerr=ddiff/Q4, xerr=dQ, fmt='.', capsize=0, + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + else: + plt.plot(Q, diff/Q4, '.', + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + if theory is not None: + mm, mp, pm, pp = theory + Q, diff, _ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + plt.plot(Q, diff/Q4, + label=self.pp.label(prefix=label, gloss='theory'), + transform=trans, + color=c['dark']) + plt.xlabel(r'Q (\AA^{-1})') + plt.ylabel(r'Difference $(R^{++} -\, R^{--}) (100 Q)^4$') + plt.legend(numpoints=1) + + def plot_FresnelDiff(self, substrate=None, surface=None, theory=None, label=None, plot_shift=None, + **kwargs): + import matplotlib.pyplot as plt + if self.pp is None or self.mm is None: + raise TypeError("cannot form difference plot without ++ and --") + if substrate is None and surface is None: + raise TypeError("Fresnel-normalized reflectivity needs substrate or surface") + F = self.fresnel(substrate=substrate, surface=surface) + #cross-sections could have different intensities - should allow for this + Q,fresnelm = self.mm.apply_beam(self.calc_Q, F(self.calc_Q),interpolation=0) + Q,fresnelp = self.pp.apply_beam(self.calc_Q, F(self.calc_Q),interpolation=0) + if substrate is None: + name = "air:%s" % surface.name + elif surface is None or isinstance(surface, Vacuum): + name = substrate.name + else: + name = "%s:%s" % (substrate.name, surface.name) + plot_shift = plot_shift if plot_shift is not None else Probe.plot_shift + trans = auto_shift(plot_shift) + pp, mm = self.pp, self.mm + c = coordinated_colors() + if hasattr(pp, 'R') and hasattr(mm, 'R') and pp.R is not None and mm.R is not None: + Q, diffF, ddiffF = splitting(pp.Q, pp.R/fresnelp, pp.dR/fresnelp, mm.Q, mm.R/fresnelm, mm.dR/fresnelm) # each cross section could have a different intensity or background theoretically, so divide first before calculating difference + if ddiffF is not None: + res = (self.show_resolution if self.show_resolution is not None + else Probe.show_resolution) + dQ = pp.dQ if res else None + plt.errorbar(Q, diffF, yerr=ddiffF, xerr=dQ, fmt='.', capsize=0, + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + else: + plt.plot(Q, diffF, '.', + label=pp.label(prefix=label, gloss='data'), + transform=trans, + color=c['light']) + if theory is not None: + mm, mp, pm, pp = theory + Q, diffF, _ = splitting(pp[0], pp[1]/fresnelp, None, mm[0], mm[1]/fresnelm, None) + plt.plot(Q, diffF, + label=self.pp.label(prefix=label, gloss='theory'), + transform=trans, + color=c['dark']) + plt.xlabel(r'Q (\AA^{-1})') + plt.ylabel(r'Difference $(R^{++} -\, R^{--})/R($'+name+')') + plt.legend(numpoints=1) + def _xs_plot(self, plotter, theory=None, **kwargs): import matplotlib.pyplot as plt # Plot available cross sections @@ -2064,69 +2198,118 @@ def plot_resolution(self, **kwargs): def plot_linear(self, **kwargs): - self._xs_plot('plot_linear', **kwargs) + self._sm_plot('plot_linear', **kwargs) def plot_log(self, **kwargs): - self._xs_plot('plot_log', **kwargs) + self._sm_plot('plot_log', **kwargs) def plot_fresnel(self, **kwargs): - self._xs_plot('plot_fresnel', **kwargs) + self._sm_plot('plot_fresnel', **kwargs) def plot_logfresnel(self, **kwargs): - self._xs_plot('plot_logfresnel', **kwargs) + self._sm_plot('plot_logfresnel', **kwargs) def plot_Q4(self, **kwargs): - self._xs_plot('plot_Q4', **kwargs) + self._sm_plot('plot_Q4', **kwargs) - #plot SA residuals NOT for other xs def plot_residuals(self, **kwargs): self._xs_plot('plot_residuals', **kwargs) + def plot_Diff(self, **kwargs): + self._df_plot('plot_linear', ylabel=r'Difference $(R^{++} -\, R^{--})$', **kwargs) + + def plot_FresnelDiff(self, **kwargs): + import matplotlib.pyplot as plt + self._df_plot('plot_fresnel', **kwargs) + ylabel = plt.gca().get_ylabel() + plt.ylabel(ylabel.replace('R/',r'Difference $(R^{++} -\, R^{--})$/')) + + def plot_DiffQ4(self, **kwargs): + import matplotlib.pyplot as plt + self._df_plot('plot_Q4', **kwargs) + plt.ylabel(ylabel=r'Difference $(R^{++} -\, R^{--}) (100 Q)^4$') def plot_SA(self, theory=None, label=None, plot_shift=None, **kwargs): import matplotlib.pyplot as plt # print('plotting Difference') - if self.df is None:# and (self.pp is None or self.mm is None): + if self.df is None or self.sm is None:# and (self.pp is None or self.mm is None): # raise TypeError("cannot form difference plot without difference") - raise TypeError("cannot form difference plot without difference") + raise TypeError("cannot form SA plot without sum and difference") plot_shift = plot_shift if plot_shift is not None else Probe.plot_shift trans = auto_shift(plot_shift) df = self.df + sm = self.sm c = coordinated_colors() if hasattr(df, 'R') and df.R is not None: #plot measured DF (not from pp, mm) Q = df.Q - DF = df.R - dDF = df.dR - if dDF is not None: + Rd = df.R + Rs = np.interp(Q, sm.Q, sm.R) + SA = Rd/Rs + dRd = df.dR + if dRd is not None: + dRs = np.interp(Q, sm.Q, sm.dR) + dvsq = (dRd/Rs)**2 + (dRs*Rd)**2/(Rs)**4 + dvsq[dvsq < 0] = 0 + dSA = sqrt(dvsq) + else: + dSA = None + + if dSA is not None: res = (self.show_resolution if self.show_resolution is not None else Probe.show_resolution) dQ = df.dQ if res else None - plt.errorbar(10*Q, ((10*Q)**4)*DF, yerr=((10*Q)**4)*dDF, xerr=10*dQ, fmt='.', capsize=0, + plt.errorbar(Q, SA, yerr=dSA, xerr=dQ, fmt='.', capsize=0, label=df.label(prefix=label, gloss='data'), transform=trans, color=c['light']) else: - plt.plot(10*Q, ((10*Q)**4)*DF, '.', + plt.plot(Q, SA, '.', label=df.label(prefix=label, gloss='data'), transform=trans, color=c['light']) if theory is not None: #plot theory curve mm, mp, pm, pp = theory - Q, DF, _ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) - plt.plot(10*Q, ((10*Q)**4)*DF, + Q, SA, _ = spin_asymmetry(pp[0], pp[1], None, mm[0], mm[1], None) + plt.plot(Q, SA, label=self.pp.label(prefix=label, gloss='theory'), transform=trans, color=c['dark']) - plt.xlabel(r'Q (nm$^{-1}$)') - plt.ylabel(r'Difference ($R^{++} -\, R^{--}$) $\times$ Q$^4$ nm$^{-4}$') + plt.xlabel(r'Q (\AA^{-1})') + plt.ylabel(r'spin asymmetry $(R^{++} -\, R^{--}) / (R^{++} +\, R^{--})$') plt.legend(numpoints=1) + def _sm_plot(self, plotter, theory=None, **kwargs): + import matplotlib.pyplot as plt + if theory is None: + theory = (None, None, None, None) + #only first four cross-sections + x_data = self.sm + mm, mp, pm, pp = theory + Q,SM,_ = meanreflectivity(pp[0], pp[1], None, mm[0], mm[1], None) + suffix = 'SM' + if x_data is not None: + fn = getattr(x_data, plotter) + fn(theory=[Q,SM], suffix=suffix, **kwargs) + + def _df_plot(self, plotter, theory=None, **kwargs): + import matplotlib.pyplot as plt + if theory is None: + theory = (None, None, None, None) + #only first four cross-sections + x_data = self.df + mm, mp, pm, pp = theory + Q,DF,_ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + suffix = 'DF' + if x_data is not None: + fn = getattr(x_data, plotter) + fn(theory=[Q,DF], suffix=suffix, **kwargs) + def _xs_plot(self, plotter, theory=None, **kwargs): import matplotlib.pyplot as plt # plot DF residuals rather than other xs residuals @@ -2173,7 +2356,17 @@ def _xs_plot(self, plotter, theory=None, **kwargs): # if x_data is not None: # fn = getattr(x_data, plotter) # fn(theory=x_th, suffix=suffix, **kwargs) - + elif 'Diff' in plotter: # plot difference not other xs + if theory is None: + theory = (None, None, None, None) + #only first four cross-sections + x_data = self.df + mm, mp, pm, pp = theory + Q,DF,_ = splitting(pp[0], pp[1], None, mm[0], mm[1], None) + suffix = 'DF' + if x_data is not None: + fn = getattr(x_data, plotter) + fn(theory=[Q,DF], suffix=suffix, **kwargs) else: # plot NR not other xs if theory is None: theory = (None, None, None, None) diff --git a/refl1d/view/data_view.py b/refl1d/view/data_view.py index e2416123c..5664740fd 100755 --- a/refl1d/view/data_view.py +++ b/refl1d/view/data_view.py @@ -94,20 +94,47 @@ def menu(self): "Plot R * Q^4") frame.Bind(wx.EVT_MENU, self.OnQ4, _item) _item.Check(Probe.view == 'q4') + + menu.Break() + # menu.AppendSeparator() + _item = menu.AppendRadioItem(wx.ID_ANY, "&SA", "Plot spin asymmetry") frame.Bind(wx.EVT_MENU, self.OnSA, _item) _item.Check(Probe.view == 'SA') - menu.AppendSeparator() + _item = menu.AppendRadioItem(wx.ID_ANY, + "Difference", + "Plot difference") + frame.Bind(wx.EVT_MENU, self.OnDiff, _item) + _item.Check(Probe.view == 'diff') + + _item = menu.AppendRadioItem(wx.ID_ANY, + "Fresnel Difference", + "Plot difference/R_F") + frame.Bind(wx.EVT_MENU, self.OnFresnelDiff, _item) + _item.Check(Probe.view == 'fresnelDiff') - _item = menu.Append(wx.ID_ANY, + _item = menu.AppendRadioItem(wx.ID_ANY, + "Difference Q4", + "Plot difference * Q^4") + frame.Bind(wx.EVT_MENU, self.OnDiffQ4, _item) + _item.Check(Probe.view == 'diffq4') + + # menu.AppendSeparator() + + _item = menu.AppendRadioItem(wx.ID_ANY, "&Residuals", "Plot residuals (R_theory - R)/dR") frame.Bind(wx.EVT_MENU, self.OnResiduals, _item) + _item.Check(Probe.view == 'residuals') + + # menu.InsertSeparator(size_t 6) + menu.Enable(id=_item.GetId(), enable=True) + return menu # ==== Views ==== @@ -135,6 +162,18 @@ def OnSA(self, event): self.view = "SA" self.redraw(reset=True) + def OnDiff(self, event): + self.view = "diff" + self.redraw(reset=True) + + def OnFresnelDiff(self, event): + self.view = "fresnelDiff" + self.redraw(reset=True) + + def OnDiffQ4(self, event): + self.view = "diffq4" + self.redraw(reset=True) + def OnResiduals(self, event): self.view = "residuals" self.redraw(reset=True)