From d16b6ad71f612c38203546aef0b36a4667df53bd Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 00:08:03 +0100 Subject: [PATCH 01/27] add new ndscan.fitting module --- ndscan/fitting/__init__.py | 2 + ndscan/fitting/fitting.py | 269 +++++++++++++++++++++++++++++++++++++ ndscan/fitting/sinusoid.py | 139 +++++++++++++++++++ 3 files changed, 410 insertions(+) create mode 100644 ndscan/fitting/__init__.py create mode 100644 ndscan/fitting/fitting.py create mode 100644 ndscan/fitting/sinusoid.py diff --git a/ndscan/fitting/__init__.py b/ndscan/fitting/__init__.py new file mode 100644 index 00000000..c8e9af68 --- /dev/null +++ b/ndscan/fitting/__init__.py @@ -0,0 +1,2 @@ +from .fitting import FitBase # noqa F401 +from .sinusoid import Sinusoid # noqa F401 diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py new file mode 100644 index 00000000..b2bbadf3 --- /dev/null +++ b/ndscan/fitting/fitting.py @@ -0,0 +1,269 @@ +import numpy as np +import functools +import collections +from typing import Dict, List, Optional, Tuple, Union +from scipy import optimize + + +# TODO: type annotations! (What's the best way of annotating np arrays?) +# TODO: docstring formatting for sphinx +class FitBase: + def __init__( + self, + x=None, + y=None, + y_err=None, + param_bounds: Optional[Dict[str, Tuple[float, float]]] = None, + fixed_params: Optional[Dict[str, float]] = None, + initial_values: Optional[Dict[str, float]] = None, + scale_factors: Optional[Dict[str, float]] = None, + ): + """ + If `x` and `y` are provided we fit the provided dataset. `x`, `y` and `y_err` + may be updated post construction to fit other datasets (c.f. :meth fit:). + + Fit results are exposed as class properties. For example, a fit parameter + 'foo' is accessed via 'fit.foo', which returns a tuple with the fitted value + of 'foo' and its estimated uncertainty. + + :param x: array of x-axis values. + :param y: array of y-axis values. + :param y_err: array of y-axis uncertainties. + :param param_bounds: dictionary of tuples containing the lower and upper bounds + for each parameter. If not specified, the defaults from + :meth get_default_bounds: are used. + :param fixed_params: dictionary specifying constant values for any non-floated + parameters. If not specified, the defaults from + :meth get_default_fixed_params: are used. + :param initial_values: dictionary specifying initial parameter values to use in + the fit. These override the values found by :meth estimate_parameters: + :param scale_factors: dictionary specifying scale factors for parameters. The + parameter values are normalised by these values during fitting to help the + curve fit (helps when fitting parameters with very large or very small + values). At some point in the future we could consider choosing sensible + default scale factors, but it's not totally trivial to do that in a way + that doesn't blow up for parameters with initial values close to zero. + """ + self._x = x + self._y = y + self._y_err = y_err + + self._validate_param_names(fixed_params) + if fixed_params is not None: + self._fixed_params = fixed_params + else: + self._fixed_params = self.get_default_fixed_params() + self._free_params = [ + param for param in self.get_params() + if param not in self._fixed_params.keys() + ] + + self._validate_param_names(param_bounds) + self._param_bounds = self.get_default_bounds() + self._param_bounds.update(param_bounds or {}) + self._param_bounds = { + param: bound + for param, bound in self._param_bounds.items() if param in self._free_params + } + + self._validate_param_names(initial_values, self._free_params) + self._initial_values = initial_values or {} + + self._validate_param_names(scale_factors, self._free_params) + self._scale_factors = scale_factors or {} + + self._p = {} + self._perr = {} + + def getter(name, instance): + return instance._p[name], instance._p_err[name] + + cls = type(self) + for param in self.get_params() + self.get_derived_params(): + setattr(cls, param, property(functools.partial(getter, param))) + + if self._x is not None and self._y is not None: + self.fit() + + def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: + """ + Fit the dataset and return the fitted parameter values and uncertainties. + """ + valid_pts = np.logical_and(np.isfinite(self._x), np.isfinite(self._y)) + x = self._x[valid_pts] + y = self._y[valid_pts] + y_err = None if self._y_err is None else self._y_err[valid_pts] + + initial_values = self._estimate_parameters() + initial_values.update(self._initial_values) + + lower_bounds = [self._param_bounds[param][0] for param in self._free_params] + upper_bounds = [self._param_bounds[param][1] for param in self._free_params] + + p0 = [initial_values[param] for param in self._free_params] + p0 = [max(estimate, lower) for estimate, lower in zip(p0, lower_bounds)] + p0 = [min(estimate, upper) for estimate, upper in zip(p0, upper_bounds)] + + scale_factors = [ + self._scale_factors.get(param, 1) for param in self._free_params + ] + scale_factors = np.asarray(scale_factors) + + p0 = np.asarray(p0) / scale_factors + lower_bounds = np.asarray(lower_bounds) / scale_factors + upper_bounds = np.asarray(upper_bounds) / scale_factors + + p, p_cov = optimize.curve_fit( + f=functools.partial(self._func_free, scale_factors), + xdata=x, + ydata=y, + p0=p0, + sigma=y_err, + absolute_sigma=True, + bounds=(lower_bounds, upper_bounds), + ) + + p *= scale_factors + p_cov *= scale_factors + + self._p_cov = p_cov + p_err = np.sqrt(np.diag(p_cov)) + + self._p = {param: value for param, value in zip(self._free_params, p)} + self._p_err = {param: value for param, value in zip(self._free_params, p_err)} + self._p.update(self._fixed_params) + self._p_err.update({param: 0 for param in self._fixed_params.keys()}) + + self._calculate_derived_params() + assert set(self.get_params() + self.get_derived_params()) == set(self._p.keys()) + return self._p, self._p_err + + def _func_free(self, scale_factors, x, *args): + """Evaluate the fit function with given values for the free parameters.""" + params = dict(self._fixed_params) + + p = np.asarray(args) + p *= scale_factors + params.update({k: v for k, v in zip(self._free_params, p)}) + + return self.func(x, params) + + def evaluate(self, x_fit: Union[np.array, int] = 100): + """Evaluates the function along x-fit and returns the tuple (x_fit, y_fit) + with the results. + + `x_fit` may either be a scalar or an array. If it is a scalar it gives the + number of equally spaced points between `min(self.x)` `max(self.x)`. Otherwise + it gives the x-axis to use. + """ + if np.isscalar(x_fit): + x_fit = np.linspace(np.min(self._x), np.max(self._x), x_fit) + y_fit = self.func(x_fit, self._p) + return x_fit, y_fit + + @property + def residuals(self): + """Returns the fit residuals.""" + return self._y - self._func_free(self._x, self._p) + + def _calculate_derived_params(self): + """ + Updates fit results with values and uncertainties for derived parameters. + """ + pass + + @staticmethod + def func(x, params): + """The fit function.""" + raise NotImplementedError + + @staticmethod + def get_params(self) -> List[str]: + """Returns list of fit params""" + raise NotImplementedError + + def get_free_params(self) -> List[str]: + """ Returns a list of the free (not held constant) fit parameters. """ + return self._free_params + + @staticmethod + def get_derived_params(self) -> List[str]: + """Returns list of derived parameters""" + return [] + + @staticmethod + def get_default_fixed_params() -> Dict[str, float]: + """Returns a dictionary mapping names of parameters which are not floated by + default to their values. + """ + return {} + + @staticmethod + def get_default_bounds() -> Dict[str, Tuple[float, float]]: + """ + Returns a dictionary mapping parameter names to a tuple of (upper, lower) bounds + """ + raise NotImplementedError + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {} + + def _estimate_parameters(self) -> Dict[str, float]: + """ + Returns a dictionary of estimates for the parameter values for the current + dataset. + """ + raise NotImplementedError + + def _validate_param_names(self, param_names, valid_params=None): + if param_names is None: + return + + valid_params = set( + valid_params if valid_params is not None else self.get_params()) + params = set(param_names) + + duplicates = [ + param for param, count in collections.Counter(params).items() if count > 1 + ] + if duplicates: + raise ValueError(f"Duplicate parameters: {','.join(duplicates)}") + + invalid = params - params.intersection(valid_params) + if invalid: + raise ValueError(f"Invalid parameters: {','.join(invalid)}") + + @property + def x(self): + """ Dataset x axis """ + return self._x + + @x.setter + def x(self, x): + self._x = np.asarray(x) + self._p = {} + self._perr = {} + + @property + def y(self): + """ Dataset y axis """ + return self._y + + @y.setter + def y(self, y): + self._y = np.asarray(y) + self._p = {} + self._perr = {} + + @property + def y_err(self): + """ Dataset y-axis uncertainty """ + return self._y_err + + @y_err.setter + def y_err(self, y_err=None): + self._y_err = np.asarray(y_err) if y_err is not None else None + self._p = {} + self._perr = {} diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py new file mode 100644 index 00000000..ba5de16f --- /dev/null +++ b/ndscan/fitting/sinusoid.py @@ -0,0 +1,139 @@ +import numpy as np +from typing import Dict, List, Tuple +from ndscan import fitting + + +class Sinusoid(fitting.FitBase): + """Sinusoid fit according to: + y = a*sin(2*np.pi*f*x + phi) + offset + + TODO: t_dead, exp decay (but default to fixed at 0) + """ + @staticmethod + def func(x, params): + return (params["a"] * np.sin(2 * np.pi * params["f"] * x + params["phi"]) + + params["offset"]) + + @staticmethod + def get_params() -> List[str]: + """Returns list of fit params""" + return ["a", "f", "phi", "offset", "t_dead"] + + @staticmethod + def get_default_bounds() -> Dict[str, Tuple[float, float]]: + """ + Returns a dictionary mapping parameter names to a tuple of (upper, lower) bounds + """ + return { + "a": (0, 1), + "f": (0, np.inf), + "phi": (0, 2 * np.pi), + "offset": (0, 1), + "t_dead": (0, np.inf) + } + + @staticmethod + def get_default_fixed_params() -> Dict[str, float]: + """Returns a dictionary mapping names of parameters which are not floated by + default to their values. + """ + return {'t_dead': 0} + + @staticmethod + def get_derived_params() -> List[str]: + """Returns a list of derived parameters. + + NB we define t_pi (t_pi_2) as the time taken for a pi-pulse (pi/2-pulse) + including deadtime. As a result, the 2*pi time is not 2*t_pi. + """ + return ["omega", 't_pi', 't_pi_2', 'phi_cosine', "min", "max", "period"] + + def _estimate_parameters(self): + # The OITG code uses a LS periodogram. I found that was giving poor phase + # estimates in some cases so here I hacked something with interpolation + and + # FFT. + # TODO: do something better! + sorted_idx = np.argsort(self.x) + x = self.x[sorted_idx] + y = self.y[sorted_idx] + + dx = np.diff(x).min() + x_axis = np.arange(x.min(), x.max(), dx) + y_interp = np.interp(x_axis, x, y) + N = len(x_axis) + + yf = np.fft.fft(y_interp)[0:N // 2] * 2.0 / N + freq = np.fft.fftfreq(N, dx)[:N // 2] + y0 = yf[0] + yf = yf[1:] + freq = freq[1:] + peak = np.argmax(np.abs(yf)) + param_guesses = {} + param_guesses["a"] = np.abs(yf[peak]) + param_guesses["f"] = freq[peak] + param_guesses["offset"] = np.abs(y0) / 2 + param_guesses["phi"] = np.angle(yf[peak]) + np.pi / 2 + + if param_guesses["phi"] < 0: + param_guesses["phi"] += 2 * np.pi + if param_guesses["phi"] > 2 * np.pi: + param_guesses["phi"] -= 2 * np.pi + + param_guesses["t_dead"] = 0 + return param_guesses + + def _calculate_derived_params(self): + self._p["omega"] = 2 * np.pi * self._p["f"] + self._p['t_pi'] = self._p["t_dead"] + np.pi / self._p["omega"] + self._p["t_pi_2"] = self._p["t_dead"] + np.pi / 2 / self._p["omega"] + self._p["phi_cosine"] = self._p["phi"] + np.pi / 2 + self._p["min"] = self._p["offset"] - np.abs(self._p["a"]) + self._p["max"] = self._p["offset"] + np.abs(self._p["a"]) + self._p["period"] = 2 * np.pi / self._p["omega"] + + # TODO: consider covariances + self._p_err["omega"] = self._p_err["f"] * 2 * np.pi + self._p_err["t_pi"] = np.sqrt(self._p_err["t_dead"]**2 + + (np.pi / self._p["omega"] * + (self._p_err["omega"] / self._p["omega"]))**2) + + self._p_err["t_pi/2"] = np.sqrt(self._p_err["t_dead"]**2 + + (np.pi / 2 / self._p["omega"] * + (self._p_err["omega"] / self._p["omega"]))**2) + + self._p_err["phi_cosine"] = self._p_err["phi"] + self._p_err["min"] = np.sqrt(self._p_err["offset"]**2 + self._p_err["a"]**2) + self._p_err["max"] = np.sqrt(self._p_err["offset"]**2 + self._p_err["a"]**2) + + self._p_err["period"] = \ + 2 * np.pi / self._p["omega"] * (self._p_err["omega"] / self._p["omega"]) + + +# TODO: move this into a unit test! +if __name__ == "__main__": + from matplotlib import pyplot as plt + + x = np.linspace(0, 5, 501) + params = { + "a": np.random.uniform(low=0.25, high=1, size=1), + "f": np.random.uniform(low=0.5, high=3, size=1), + "phi": np.random.uniform(low=0, high=2 * np.pi, size=1), + "offset": np.random.uniform(low=0.25, high=0.75, size=1), + } + y = Sinusoid.func(x, params) + np.random.normal(0, 0.05, x.shape) + fit = Sinusoid(x, y) + + if not np.isclose(fit.a[0], params["a"], rtol=5e-2): + print(f"Amplitude error is {100*(1-fit.a[0]/params['a'])}%") + if not np.isclose(fit.f[0], params["f"], rtol=5e-2): + print(f"Frequency error is {100*(1-fit.f[0]/params['f'])}%") + if not np.isclose(fit.phi[0], params["phi"], rtol=5e-2): + print(f"Phase error is {100*(1-fit.phi[0]/params['phi'])}%") + if not np.isclose(fit.offset[0], params["offset"], rtol=5e-2): + print(f"Offset error is {100*(1-fit.offset[0]/params['offset'])}%") + + plt.plot(x, y, label="data") + x_fit, y_fit = fit.evaluate() + plt.plot(x_fit, y_fit, label="fit") + plt.grid() + plt.show() From 41dccc6173f618593638dda97c54768831f4171d Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 00:10:39 +0100 Subject: [PATCH 02/27] use new fitting framework --- ndscan/experiment/default_analysis.py | 292 +++++++++++++++++--------- ndscan/experiment/utils.py | 10 +- ndscan/plots/annotation_items.py | 15 +- ndscan/plots/model/__init__.py | 5 +- ndscan/plots/model/online_analysis.py | 49 +++-- ndscan/plots/xy_1d.py | 28 +-- ndscan/utils.py | 31 ++- 7 files changed, 276 insertions(+), 154 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 2601a110..e14d9483 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -16,16 +16,20 @@ Both can produce annotations; particular values or plot locations highlighted in the user interface. """ +import numpy as np import logging from typing import Any, Callable, Dict, List, Iterable, Optional, Set, Tuple, Union +import collections +from enum import Enum +import dataclasses -from ..utils import FIT_OBJECTS from .parameters import ParamHandle -from .result_channels import ResultChannel +from .result_channels import ResultChannel, FloatChannel +from .. import utils __all__ = [ "Annotation", "DefaultAnalysis", "CustomAnalysis", "OnlineFit", - "ResultPrefixAnalysisWrapper" + "ResultPrefixAnalysisWrapper", "ExportedResult", "NamedFit" ] logger = logging.getLogger(__name__) @@ -39,6 +43,49 @@ def __init__(self, kind: str, **kwargs): self.spec = {"kind": kind, **kwargs} +class AnalysisType(Enum): + NAMED_FIT = 'named_fit' + + +@dataclasses.dataclass +class Analysis: + kind: AnalysisType = dataclasses.field(init=False) + + def __post_init__(self): + if type(self) == NamedFit: + self.kind = AnalysisType.NAMED_FIT + + +@dataclasses.dataclass +class NamedFit(Analysis): + """ Named fit function used in online analyses. + + Attributes: + fit_class_name: name of the fit class to use + fit_module: module that fit_class resides in (must be in the current python path) + data: maps fit data axis names (``"x"``, ``"y"``) to parameter handles or result + channels that supply the respective data. + param_bounds: dictionary of tuples containing the lower and upper bounds for + each parameter. If not specified, the defaults from the fit class are used. + fixed_params: dictionary specifying constant values for any non-floated + parameters. If not specified, the defaults from the fit class are used. + initial_values: dictionary specifying initial parameter values to use in + the fit. These override the values found by the fit's parameter estimator. + scale_factors: dictionary specifying scale factors for parameters. The + parameter values are normalised by these values during fitting to help the + curve fit (helps when fitting parameters with very large or very small + values). + """ + + fit_class_name: str + fit_module: str + data: dict + param_bounds: Dict[str, Tuple[float, float]] + fixed_params: Dict[str, float] + initial_values: Dict[str, float] + scale_factors: Dict[str, float] + + class AnnotationContext: """Resolves entities in user-specified annotation schemata to stringly-typed dictionary form. @@ -126,8 +173,8 @@ def required_axes(self) -> Set[ParamHandle]: raise NotImplementedError def describe_online_analyses( - self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Dict[str, Any]]]: + self, context: AnnotationContext + ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: """Exceute analysis and serialise information about resulting annotations and online analyses to stringly typed metadata. @@ -211,8 +258,8 @@ def required_axes(self) -> Set[ParamHandle]: return self._required_axis_handles def describe_online_analyses( - self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Dict[str, Any]]]: + self, context: AnnotationContext + ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: "" return [], {} @@ -249,103 +296,130 @@ def execute( return [a.describe(context) for a in annotations] -#: Default points of interest for various fit types (e.g. highlighting the π time for a -#: Rabi flop fit, or the extremum of a parabola. -DEFAULT_FIT_ANNOTATIONS = { - "decaying_sinusoid": { - "pi_time": { - "x": "t_max_transfer" - } - }, - "detuned_square_pulse": { - "centre": { - "x": "offset" - } - }, - "exponential_decay": { - "t_1_e": { - "x": "t_1_e" - } - }, - "gaussian": { - "centre": { - "x": "x0" - } - }, - "lorentzian": { - "extremum": { - "x": "x0" - } - }, - "parabola": { - "extremum": { - "x": "position" - } - }, - "rabi_flop": { - "pi_time": { - "x": "t_pi" - } - }, - "sinusoid": { - "pi_time": { - "x": "t_pi" - } - }, - "v_function": { - "centre": { - "x": "x0" - } - }, -} +@dataclasses.dataclass +class ExportedResult: + """ Analysis result exported from an online fit. + + Attributes: + fit_parameter: fit parameter to export + result_name: name of the new result channel to create. Defaults to + :param fit_parameter: if not specified. + export_err: if True, we export the uncertainty in the fitted parameter value + as an additional result with the name (:param result_name: + `_err`) + result_type: the result channel class for the new result + config: dictionary of kwargs to pass into the constructor of the result + channel + """ + fit_parameter: str + result_name: Optional[str] = None + export_err: bool = True + result_type: ResultChannel = FloatChannel + config: Optional[Dict] = None + + def __post_init__(self): + self.result_name = self.result_name or self.fit_parameter + self.config = self.config or {} class OnlineFit(DefaultAnalysis): """Describes an automatically executed fit for a given combination of scan axes and result channels. - :param fit_type: Fitting procedure name, per :data:`.FIT_OBJECTS`. + :param fit_class: name of the python class within :param fit_module: to use for the + fit. :param data: Maps fit data axis names (``"x"``, ``"y"``) to parameter handles or result channels that supply the respective data. :param annotations: Any points of interest to highlight in the fit results, given in the form of a dictionary mapping (arbitrary) identifiers to dictionaries mapping coordinate names to fit result names. If ``None``, - :data:`DEFAULT_FIT_ANNOTATIONS` will be queried. + the defaults provided by the fit function will be used. :param analysis_identifier: Optional explicit name to use for online analysis. Defaults to ``fit_``, but can be set explicitly to allow more than one fit of a given type at a time. - :param constants: Specifies parameters to be held constant during the fit. This is - a dictionary mapping fit parameter names to the respective constant values, - forwarded to :meth:`oitg.fitting.FitBase.FitBase.fit`. - :param initial_values: Specifies initial values for the fit parameters. This is - a dictionary mapping fit parameter names to the respective values, forwarded to - :meth:`oitg.fitting.FitBase.FitBase.fit`. + :param constants: dictionary specifying constant values for any non-floated + parameters. If not specified, the defaults from the fit class are used. + :param initial_values: dictionary specifying initial parameter values to use in + the fit. These override the values found by the fit's parameter estimator. + :param bounds: dictionary of tuples containing the lower and upper bounds for + each parameter. If not specified, the defaults from the fit class are used. + :param exported_results: Specifies fitted parameter values to export as analysis + results. + :param fit_module: python module containing the fit class. Will default to + `ndscan.fitting` in the future. + :param scale_factors: dictionary specifying scale factors for parameters. The + parameter values are normalised by these values during fitting to help the curve + fit (helps when fitting parameters with very large or very small values). At + some point in the future we could consider choosing sensible default scale + factors, but it's not totally trivial to do that in a way that doesn't blow up + for parameters with initial values close to zero. """ def __init__(self, - fit_type: str, + fit_class: str, data: Dict[str, Union[ParamHandle, ResultChannel]], annotations: Optional[Dict[str, Dict[str, Any]]] = None, analysis_identifier: str = None, constants: Optional[Dict[str, Any]] = None, - initial_values: Optional[Dict[str, Any]] = None): - self.fit_type = fit_type - if fit_type not in FIT_OBJECTS: - logger.warning("Unknown fit type: '%s'", fit_type, exc_info=True) + initial_values: Optional[Dict[str, Any]] = None, + bounds: Optional[Dict[str, Tuple[float, float]]] = None, + exported_results: Optional[List[ExportedResult]] = None, + fit_module: str = 'oitg.fitting', + scale_factors: Optional[Dict[str, float]] = None): + + self.fit_class_name = fit_class + self.fit_module = fit_module self.data = data - if annotations is None: - annotations = DEFAULT_FIT_ANNOTATIONS.get(fit_type, {}) self.annotations = annotations self.analysis_identifier = analysis_identifier - self.constants = {} if constants is None else constants - self.initial_values = {} if initial_values is None else initial_values + self.initial_values = initial_values or {} + self.exported_results = exported_results or [] + self.scale_factors = scale_factors or {} + + self._result_channels = [] + for result in self.exported_results: + channel = result.result_type(path=result.result_name, **result.config) + self._result_channels.append(channel) + + if result.export_err: + err_channel = FloatChannel( + path=result.result_name + '_err', + min=0.0, + display_hints={"error_bar_for": channel.path}) + self._result_channels.append(err_channel) + + duplicate_channels = [ + channel.path + for channel, count in collections.Counter(self._result_channels).items() + if count > 1 + ] + if duplicate_channels: + raise ValueError( + f"Duplicate result channels: {','.join(duplicate_channels)}") + self._result_channels = { + channel.path: channel + for channel in self._result_channels + } + + klass = utils.import_class(self.fit_module, self.fit_class_name) + self.fit_klass = klass + self.bounds = bounds if bounds is not None else klass.get_default_bounds() + + if annotations is not None: + self.annotations = annotations + else: + self.annotations = self.fit_klass.get_default_annotations() + + if constants is not None: + self.constants = constants + else: + self.constants = self.fit_klass.get_default_fixed_params() def required_axes(self) -> Set[ParamHandle]: "" return set(a for a in self.data.values() if isinstance(a, ParamHandle)) def describe_online_analyses( - self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Dict[str, Any]]]: + self, context: AnnotationContext + ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: "" # TODO: Generalise to higher-dimensional fits. channels = [ @@ -358,23 +432,25 @@ def describe_online_analyses( # By default, mangle fit type and channels into a pseudo-unique identifier, # which should work for the vast majority of cases (i.e. unless the user # creates needlessly duplicate analyses). - analysis_identifier = "fit_" + self.fit_type + "_" + "_".join(channels) + analysis_identifier = ("fit_" + f"{self.fit_module}.{self.fit_class_name}" + + "_" + "_".join(channels)) def analysis_ref(key): return AnnotationValueRef("online_result", analysis_name=analysis_identifier, result_key=key) + fit_params = self.fit_klass.get_params() + annotations = [ Annotation("computed_curve", parameters={ - "function_name": self.fit_type, - "associated_channels": channels + "fit_class_name": self.fit_class_name, + "fit_module": self.fit_module, + "associated_channels": channels, }, - data={ - k: analysis_ref(k) - for k in FIT_OBJECTS[self.fit_type].parameter_names - }) + data={param: analysis_ref(param) + for param in fit_params}) ] for a in self.annotations.values(): # TODO: Change API to allow more general annotations. @@ -389,23 +465,25 @@ def analysis_ref(key): }, parameters={"associated_channels": channels})) - return [a.describe(context) for a in annotations], { - analysis_identifier: { - "kind": "named_fit", - "fit_type": self.fit_type, - "data": { - name: context.describe_coordinate(obj) - for name, obj in self.data.items() - }, - "constants": self.constants, - "initial_values": self.initial_values - } + analysis = { + analysis_identifier: + NamedFit(fit_class_name=self.fit_class_name, + fit_module=self.fit_module, + data={ + name: context.describe_coordinate(obj) + for name, obj in self.data.items() + }, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + scale_factors=self.scale_factors) } + return [a.describe(context) for a in annotations], analysis def get_analysis_results(self) -> Dict[str, ResultChannel]: "" # Could return annotation locations in the future. - return {} + return self._result_channels def execute( self, @@ -414,7 +492,27 @@ def execute( context: AnnotationContext, ) -> List[Dict[str, Any]]: "" - # Nothing to do off-line for online fits. + user_axis_data = {} + for handle in self.required_axes(): + user_axis_data[handle] = axis_data[handle._store.identity] + # TODO: Generalise to higher-dimensional fits. + x = axis_data[self.data['x']._store.identity] + y = result_data[self.data['y']] + + fit = self.fit_klass(x=np.asarray(x), + y=np.asarray(y), + y_err=None, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + scale_factors=self.scale_factors) + for result in self.exported_results: + p, p_err = getattr(fit, result.fit_parameter) + channel = self._result_channels[result.result_name] + channel.push(p) + if result.export_err: + err_channel = self._result_channels[result.result_name + '_err'] + err_channel.push(p_err) return [] @@ -437,8 +535,8 @@ def required_axes(self) -> Set[ParamHandle]: return self._wrapped.required_axes() def describe_online_analyses( - self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Dict[str, Any]]]: + self, context: AnnotationContext + ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: return self._wrapped.describe_online_analyses(context) def get_analysis_results(self) -> Dict[str, ResultChannel]: diff --git a/ndscan/experiment/utils.py b/ndscan/experiment/utils.py index 483a0fca..ada94cdd 100644 --- a/ndscan/experiment/utils.py +++ b/ndscan/experiment/utils.py @@ -1,6 +1,8 @@ import json import numpy from typing import Any, Iterable, Optional +import dataclasses +import enum def path_matches_spec(path: Iterable[str], spec: str) -> bool: @@ -20,7 +22,7 @@ def is_kernel(func) -> bool: return meta.core_name is not None and not meta.portable -class NumpyToVanillaEncoder(json.JSONEncoder): +class Encoder(json.JSONEncoder): def default(self, obj): if isinstance(obj, numpy.integer): return int(obj) @@ -28,6 +30,10 @@ def default(self, obj): return float(obj) if isinstance(obj, numpy.ndarray): return obj.tolist() + if dataclasses.is_dataclass(obj): + return dataclasses.asdict(obj) + if isinstance(obj, enum.Enum): + return obj.value return json.JSONEncoder.default(self, obj) @@ -35,7 +41,7 @@ def dump_json(obj: Any) -> str: """Serialise ``obj`` as a JSON string, with NumPy numerical/array types encoded as their vanilla Python counterparts. """ - return json.dumps(obj, cls=NumpyToVanillaEncoder) + return json.dumps(obj, cls=Encoder) def to_metadata_broadcast_type(obj: Any) -> Optional[Any]: diff --git a/ndscan/plots/annotation_items.py b/ndscan/plots/annotation_items.py index cf298a1d..13be06bc 100644 --- a/ndscan/plots/annotation_items.py +++ b/ndscan/plots/annotation_items.py @@ -9,7 +9,7 @@ import pyqtgraph from qasync import QtCore from typing import Dict, Union, Optional, Tuple -from ..utils import FIT_OBJECTS +from .. import fitting from .model import AnnotationDataSource logger = logging.getLogger(__name__) @@ -26,8 +26,7 @@ class ComputedCurveItem(AnnotationItem): """Shows a curve (pyqtgraph.LineItem) that is computed from a given fit function, dynamically adapting to the coordinate region displayed. - :param function_name: The name of the function (see :data:`FIT_OBJECTS`) to - evaluate. + :param fit_obj: Fit object (see ndscan.fitting.FitBase) :param data_sources: A dictionary giving the parameters for the curve. :param view_box: The :class:`pyqtgraph.ViewBox` to add the line item to once there is data. @@ -37,14 +36,10 @@ class ComputedCurveItem(AnnotationItem): :param x_limits: Limits to restrict the drawn horizontal range to even if the viewport extends beyond them. """ - @staticmethod - def is_function_supported(function_name: str) -> bool: - return function_name in FIT_OBJECTS - - def __init__(self, function_name: str, + def __init__(self, fit_obj: fitting.FitBase, data_sources: Dict[str, AnnotationDataSource], view_box, curve_item, x_limits: Tuple[Union[float, None], Union[float, None]]): - self._function = FIT_OBJECTS[function_name].fitting_function + self._fit_obj = fit_obj self._data_sources = data_sources self._view_box = view_box self._curve_item = curve_item @@ -96,7 +91,7 @@ def _redraw(self, *args): # Choose number of points based on width of plot on screen (in pixels). fn_xs = numpy.linspace(*x_lims, int(self._view_box.width())) - fn_ys = self._function(fn_xs, params) + fn_ys = self._fit_obj.func(fn_xs, params) self._curve_item.setData(fn_xs, fn_ys) diff --git a/ndscan/plots/model/__init__.py b/ndscan/plots/model/__init__.py index e9d5e16d..ded78733 100644 --- a/ndscan/plots/model/__init__.py +++ b/ndscan/plots/model/__init__.py @@ -28,6 +28,7 @@ from qasync import QtCore from typing import Any, Callable, Dict, List, Optional from .online_analysis import OnlineNamedFitAnalysis +from ... import experiment logger = logging.getLogger(__name__) @@ -250,7 +251,9 @@ def _set_online_analyses(self, analysis_schemata: Dict[str, Dict[str, for name, schema in analysis_schemata.items(): kind = schema["kind"] if kind == "named_fit": - self._online_analyses[name] = OnlineNamedFitAnalysis(schema, self) + del schema["kind"] + analysis = experiment.NamedFit(**schema) + self._online_analyses[name] = OnlineNamedFitAnalysis(analysis, self) else: logger.warning("Ignoring unsupported online analysis type: '%s'", kind) diff --git a/ndscan/plots/model/online_analysis.py b/ndscan/plots/model/online_analysis.py index d40c4923..beb4b7f9 100644 --- a/ndscan/plots/model/online_analysis.py +++ b/ndscan/plots/model/online_analysis.py @@ -1,9 +1,12 @@ import asyncio +import logging from concurrent.futures import ProcessPoolExecutor from pyqtgraph import SignalProxy from qasync import QtCore -from typing import Any, Dict -from ...utils import FIT_OBJECTS +from ...experiment import NamedFit +from ... import utils + +logger = logging.getLogger(__name__) class OnlineAnalysis(QtCore.QObject): @@ -23,15 +26,16 @@ class OnlineNamedFitAnalysis(OnlineAnalysis): """ _trigger_recompute_fit = QtCore.pyqtSignal() - def __init__(self, schema: Dict[str, Any], parent_model): + def __init__(self, analysis: NamedFit, parent_model): super().__init__() - self._schema = schema self._model = parent_model + self._analysis = analysis - self._fit_type = self._schema["fit_type"] - self._fit_obj = FIT_OBJECTS[self._fit_type] - self._constants = self._schema.get("constants", {}) - self._initial_values = self._schema.get("initial_values", {}) + fit_class = utils.import_class(analysis.fit_module, analysis.fit_class_name) + self._fit_obj = fit_class(param_bounds=self._analysis.param_bounds, + fixed_params=self._analysis.fixed_params, + initial_values=self._analysis.initial_values, + scale_factors=self._analysis.scale_factors) self._last_fit_params = None self._last_fit_errors = None @@ -69,12 +73,12 @@ def _update(self): data = self._model.get_point_data() self._source_data = {} - for param_key, source_key in self._schema["data"].items(): + for param_key, source_key in self._analysis.data.items(): self._source_data[param_key] = data.get(source_key, []) # Truncate the source data to a complete set of points. num_points = min(len(v) for v in self._source_data.values()) - if num_points < len(self._fit_obj.parameter_names): + if num_points < len(self._fit_obj.get_free_params()): # Not enough points yet for the given number of degrees of freedom. return @@ -100,25 +104,28 @@ async def _recompute_fit(self): y_errs = self._source_data.get("y_err", None) loop = asyncio.get_event_loop() - self._last_fit_params, self._last_fit_errors = await loop.run_in_executor( - self._fit_executor, _run_fit, self._fit_type, xs, ys, y_errs, - self._constants, self._initial_values) + + self._fit_obj.x = xs + self._fit_obj.y = ys + self._fit_obj.y_err = y_errs + + fit = await loop.run_in_executor(self._fit_executor, _run_fit, self._fit_obj) + self._last_fit_params, self._last_fit_errors = fit self._recompute_in_progress = False self.updated.emit() -def _run_fit(fit_type, xs, ys, y_errs, constants, initial_values): +def _run_fit(fit_obj): """Fits the given data with the chosen method. This function is intended to be executed on a worker process, hence the primitive API. """ try: - return FIT_OBJECTS[fit_type].fit(x=xs, - y=ys, - y_err=y_errs, - constants=constants, - initialise=initial_values) - except Exception: - return None, None + return fit_obj.fit() + except Exception as e: + logger.warning(f"Exception encountered in fit function: {e}") + pass + + return None, None diff --git a/ndscan/plots/xy_1d.py b/ndscan/plots/xy_1d.py index e688287c..6cf34303 100644 --- a/ndscan/plots/xy_1d.py +++ b/ndscan/plots/xy_1d.py @@ -12,6 +12,7 @@ from .utils import (extract_linked_datasets, extract_scalar_channels, format_param_identity, group_channels_into_axes, setup_axis_item, FIT_COLORS, SERIES_COLORS) +from ..utils import import_class logger = logging.getLogger(__name__) @@ -247,18 +248,21 @@ def make_curve_item(series_idx): continue if a.kind == "computed_curve": - function_name = a.parameters.get("function_name", None) - if ComputedCurveItem.is_function_supported(function_name): - associated_series_idx = max( - channel_ref_to_series_idx(chan) - for chan in a.parameters.get("associated_channels", [None])) - - x_limits = [self.x_param_spec.get(n, None) for n in ("min", "max")] - curve = make_curve_item(associated_series_idx) - vb = self.series[associated_series_idx].view_box - item = ComputedCurveItem(function_name, a.data, vb, curve, x_limits) - self.annotation_items.append(item) - continue + fit_module = a.parameters["fit_module"] + fit_class_name = a.parameters["fit_class_name"] + fit_class = import_class(fit_module, fit_class_name) + fit_obj = fit_class() + + associated_series_idx = max( + channel_ref_to_series_idx(chan) + for chan in a.parameters.get("associated_channels", [None])) + + x_limits = [self.x_param_spec.get(n, None) for n in ("min", "max")] + curve = make_curve_item(associated_series_idx) + vb = self.series[associated_series_idx].view_box + item = ComputedCurveItem(fit_obj, a.data, vb, curve, x_limits) + self.annotation_items.append(item) + continue logger.info("Ignoring annotation of kind '%s' with coordinates %s", a.kind, list(a.coordinates.keys())) diff --git a/ndscan/utils.py b/ndscan/utils.py index 309fc33d..124f30fd 100644 --- a/ndscan/utils.py +++ b/ndscan/utils.py @@ -1,19 +1,10 @@ """Odds and ends common to all of ndscan.""" from enum import Enum, unique -import oitg.fitting +import importlib +import inspect from typing import Any, Callable, Dict, Iterable, Optional, Protocol, TypeVar -#: Registry of well-known fit procedure names. -FIT_OBJECTS = { - n: getattr(oitg.fitting, n) - for n in [ - "cos", "decaying_sinusoid", "detuned_square_pulse", "exponential_decay", - "gaussian", "line", "lorentzian", "rabi_flop", "sinusoid", "v_function" - ] -} -FIT_OBJECTS["parabola"] = oitg.fitting.shifted_parabola - #: Name of the ``artiq.language.HasEnvironment`` argument that is used to confer the #: list of available parameters to the dashboard plugin, and to pass the information #: about scanned and overridden parameters to the :class:`FragmentScanExperiment` @@ -117,3 +108,21 @@ def merge_no_duplicates(target: dict, source: dict, kind: str = "entries") -> No raise ValueError("Duplicate {} of key '{}'".format(kind, k)) target[k] = v return target + + +def import_class(module_name: str, class_name: str): + """ + Imports a named class from a module in the python path or raises an exception. + """ + try: + module = importlib.import_module(module_name) + except ImportError: + raise ValueError(f'Cannot import module "{module_name}"') + + module_classes = [ + name for name, obj in inspect.getmembers(module) if inspect.isclass(obj) + ] + if class_name not in module_classes: + raise ValueError(f'Class "{class_name}" not in module "{module_name}"') + + return getattr(module, class_name) From 815f6c3a6c200a77ed4beda6963c2ca341f8ac6f Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 20:46:24 +0100 Subject: [PATCH 03/27] NamedFit -> FitDescription, remove the Analysis base class --- ndscan/experiment/default_analysis.py | 46 +-- ndscan/plots/model/__init__.py | 9 +- ndscan/plots/model/online_analysis.py | 4 +- poetry.lock | 431 ++++++++++++++++++++++++++ pyproject.toml | 42 +++ 5 files changed, 494 insertions(+), 38 deletions(-) create mode 100644 poetry.lock create mode 100644 pyproject.toml diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index e14d9483..8b765e4b 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -29,7 +29,7 @@ __all__ = [ "Annotation", "DefaultAnalysis", "CustomAnalysis", "OnlineFit", - "ResultPrefixAnalysisWrapper", "ExportedResult", "NamedFit" + "ResultPrefixAnalysisWrapper", "ExportedResult", "FitDescription" ] logger = logging.getLogger(__name__) @@ -43,22 +43,9 @@ def __init__(self, kind: str, **kwargs): self.spec = {"kind": kind, **kwargs} -class AnalysisType(Enum): - NAMED_FIT = 'named_fit' - - -@dataclasses.dataclass -class Analysis: - kind: AnalysisType = dataclasses.field(init=False) - - def __post_init__(self): - if type(self) == NamedFit: - self.kind = AnalysisType.NAMED_FIT - - @dataclasses.dataclass -class NamedFit(Analysis): - """ Named fit function used in online analyses. +class FitDescription(): + """ Description of an online analysis fit. Attributes: fit_class_name: name of the fit class to use @@ -174,7 +161,7 @@ def required_axes(self) -> Set[ParamHandle]: def describe_online_analyses( self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: + ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: """Exceute analysis and serialise information about resulting annotations and online analyses to stringly typed metadata. @@ -259,7 +246,7 @@ def required_axes(self) -> Set[ParamHandle]: def describe_online_analyses( self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: + ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: "" return [], {} @@ -419,7 +406,7 @@ def required_axes(self) -> Set[ParamHandle]: def describe_online_analyses( self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: + ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: "" # TODO: Generalise to higher-dimensional fits. channels = [ @@ -467,16 +454,17 @@ def analysis_ref(key): analysis = { analysis_identifier: - NamedFit(fit_class_name=self.fit_class_name, - fit_module=self.fit_module, - data={ - name: context.describe_coordinate(obj) - for name, obj in self.data.items() + FitDescription( + fit_class_name=self.fit_class_name, + fit_module=self.fit_module, + data={name: context.describe_coordinate(obj) + for name, obj in self.data.items() }, - param_bounds=self.bounds, - fixed_params=self.constants, - initial_values=self.initial_values, - scale_factors=self.scale_factors) + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + scale_factors=self.scale_factors + ) } return [a.describe(context) for a in annotations], analysis @@ -536,7 +524,7 @@ def required_axes(self) -> Set[ParamHandle]: def describe_online_analyses( self, context: AnnotationContext - ) -> Tuple[List[Dict[str, Any]], Dict[str, Analysis]]: + ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: return self._wrapped.describe_online_analyses(context) def get_analysis_results(self) -> Dict[str, ResultChannel]: diff --git a/ndscan/plots/model/__init__.py b/ndscan/plots/model/__init__.py index ded78733..97ee2eaf 100644 --- a/ndscan/plots/model/__init__.py +++ b/ndscan/plots/model/__init__.py @@ -249,13 +249,8 @@ def _set_online_analyses(self, analysis_schemata: Dict[str, Dict[str, self._online_analyses = {} for name, schema in analysis_schemata.items(): - kind = schema["kind"] - if kind == "named_fit": - del schema["kind"] - analysis = experiment.NamedFit(**schema) - self._online_analyses[name] = OnlineNamedFitAnalysis(analysis, self) - else: - logger.warning("Ignoring unsupported online analysis type: '%s'", kind) + analysis = experiment.NamedFit(**schema) + self._online_analyses[name] = OnlineNamedFitAnalysis(analysis, self) # Rebind annotation schemata to new analysis data sources. self._set_annotation_schemata(self._annotation_schemata) diff --git a/ndscan/plots/model/online_analysis.py b/ndscan/plots/model/online_analysis.py index beb4b7f9..6b23d2b6 100644 --- a/ndscan/plots/model/online_analysis.py +++ b/ndscan/plots/model/online_analysis.py @@ -3,7 +3,7 @@ from concurrent.futures import ProcessPoolExecutor from pyqtgraph import SignalProxy from qasync import QtCore -from ...experiment import NamedFit +from ...experiment import FitDescription from ... import utils logger = logging.getLogger(__name__) @@ -26,7 +26,7 @@ class OnlineNamedFitAnalysis(OnlineAnalysis): """ _trigger_recompute_fit = QtCore.pyqtSignal() - def __init__(self, analysis: NamedFit, parent_model): + def __init__(self, analysis: FitDescription, parent_model): super().__init__() self._model = parent_model self._analysis = analysis diff --git a/poetry.lock b/poetry.lock new file mode 100644 index 00000000..34fc6f33 --- /dev/null +++ b/poetry.lock @@ -0,0 +1,431 @@ +[[package]] +name = "artiq" +version = "7.0b0" +description = "Advanced Real-Time Infrastructure for Quantum physics" +category = "main" +optional = false +python-versions = "*" +develop = false + +[package.dependencies] +h5py = "*" +llvmlite = "*" +numpy = "*" +prettytable = "*" +pygit2 = "*" +pyqtgraph = "*" +python-dateutil = "*" +python-Levenshtein = "*" +pythonparser = "*" +qasync = "*" +scipy = "*" + +[package.source] +type = "git" +url = "https://github.com/m-labs/artiq.git" +reference = "master" +resolved_reference = "f31279411ef7f198ee697290efe957544896ebdb" + +[[package]] +name = "cffi" +version = "1.15.1" +description = "Foreign Function Interface for Python calling C code." +category = "main" +optional = false +python-versions = "*" + +[package.dependencies] +pycparser = "*" + +[[package]] +name = "flake8" +version = "4.0.1" +description = "the modular source code checker: pep8 pyflakes and co" +category = "dev" +optional = false +python-versions = ">=3.6" + +[package.dependencies] +mccabe = ">=0.6.0,<0.7.0" +pycodestyle = ">=2.8.0,<2.9.0" +pyflakes = ">=2.4.0,<2.5.0" + +[[package]] +name = "h5py" +version = "3.7.0" +description = "Read and write HDF5 files from Python" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +numpy = ">=1.14.5" + +[[package]] +name = "llvmlite" +version = "0.39.0" +description = "lightweight wrapper around basic LLVM functionality" +category = "main" +optional = false +python-versions = ">=3.7" + +[[package]] +name = "mccabe" +version = "0.6.1" +description = "McCabe checker, plugin for flake8" +category = "dev" +optional = false +python-versions = "*" + +[[package]] +name = "numpy" +version = "1.23.1" +description = "NumPy is the fundamental package for array computing with Python." +category = "main" +optional = false +python-versions = ">=3.8" + +[[package]] +name = "oitg" +version = "0.1" +description = "" +category = "main" +optional = false +python-versions = "*" +develop = false + +[package.dependencies] +statsmodels = "*" + +[package.source] +type = "git" +url = "https://github.com/oxfordiontrapgroup/oitg.git" +reference = "master" +resolved_reference = "937dfaa16a261e60f24c0165ddfc7818756e1683" + +[[package]] +name = "packaging" +version = "21.3" +description = "Core utilities for Python packages" +category = "main" +optional = false +python-versions = ">=3.6" + +[package.dependencies] +pyparsing = ">=2.0.2,<3.0.5 || >3.0.5" + +[[package]] +name = "pandas" +version = "1.4.3" +description = "Powerful data structures for data analysis, time series, and statistics" +category = "main" +optional = false +python-versions = ">=3.8" + +[package.dependencies] +numpy = [ + {version = ">=1.18.5", markers = "platform_machine != \"aarch64\" and platform_machine != \"arm64\" and python_version < \"3.10\""}, + {version = ">=1.19.2", markers = "platform_machine == \"aarch64\" and python_version < \"3.10\""}, + {version = ">=1.20.0", markers = "platform_machine == \"arm64\" and python_version < \"3.10\""}, +] +python-dateutil = ">=2.8.1" +pytz = ">=2020.1" + +[package.extras] +test = ["hypothesis (>=5.5.3)", "pytest (>=6.0)", "pytest-xdist (>=1.31)"] + +[[package]] +name = "patsy" +version = "0.5.2" +description = "A Python package for describing statistical models and for building design matrices." +category = "main" +optional = false +python-versions = "*" + +[package.dependencies] +numpy = ">=1.4" +six = "*" + +[package.extras] +test = ["scipy", "pytest-cov", "pytest"] + +[[package]] +name = "prettytable" +version = "3.3.0" +description = "A simple Python library for easily displaying tabular data in a visually appealing ASCII table format" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +wcwidth = "*" + +[package.extras] +tests = ["pytest", "pytest-cov", "pytest-lazy-fixture"] + +[[package]] +name = "pycodestyle" +version = "2.8.0" +description = "Python style guide checker" +category = "dev" +optional = false +python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*" + +[[package]] +name = "pycparser" +version = "2.21" +description = "C parser in Python" +category = "main" +optional = false +python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" + +[[package]] +name = "pyflakes" +version = "2.4.0" +description = "passive checker of Python programs" +category = "dev" +optional = false +python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" + +[[package]] +name = "pygit2" +version = "1.10.0" +description = "Python bindings for libgit2." +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +cffi = ">=1.9.1" + +[[package]] +name = "pyparsing" +version = "3.0.9" +description = "pyparsing module - Classes and methods to define and execute parsing grammars" +category = "main" +optional = false +python-versions = ">=3.6.8" + +[package.extras] +diagrams = ["railroad-diagrams", "jinja2"] + +[[package]] +name = "pyqt5" +version = "5.15.7" +description = "Python bindings for the Qt cross platform application toolkit" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +PyQt5-Qt5 = ">=5.15.0" +PyQt5-sip = ">=12.11,<13" + +[[package]] +name = "pyqt5-qt5" +version = "5.15.2" +description = "The subset of a Qt installation needed by PyQt5." +category = "main" +optional = false +python-versions = "*" + +[[package]] +name = "pyqt5-sip" +version = "12.11.0" +description = "The sip module support for PyQt5" +category = "main" +optional = false +python-versions = ">=3.7" + +[[package]] +name = "pyqtgraph" +version = "0.12.4" +description = "Scientific Graphics and GUI Library for Python" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +numpy = ">=1.17.0" + +[[package]] +name = "python-dateutil" +version = "2.8.2" +description = "Extensions to the standard Python datetime module" +category = "main" +optional = false +python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7" + +[package.dependencies] +six = ">=1.5" + +[[package]] +name = "python-levenshtein" +version = "0.12.2" +description = "Python extension for computing string edit distances and similarities." +category = "main" +optional = false +python-versions = "*" + +[[package]] +name = "pythonparser" +version = "1.3" +description = "A Python parser intended for use in tooling" +category = "main" +optional = false +python-versions = "*" + +[package.dependencies] +regex = "*" + +[[package]] +name = "pytz" +version = "2022.1" +description = "World timezone definitions, modern and historical" +category = "main" +optional = false +python-versions = "*" + +[[package]] +name = "qasync" +version = "0.23.0" +description = "Implementation of the PEP 3156 Event-Loop with Qt." +category = "main" +optional = false +python-versions = "~=3.6" + +[[package]] +name = "regex" +version = "2022.7.25" +description = "Alternative regular expression module, to replace re." +category = "main" +optional = false +python-versions = ">=3.6" + +[[package]] +name = "scipy" +version = "1.9.0" +description = "SciPy: Scientific Library for Python" +category = "main" +optional = false +python-versions = ">=3.8,<3.12" + +[package.dependencies] +numpy = ">=1.18.5,<1.25.0" + +[[package]] +name = "sipyco" +version = "1.4" +description = "Simple Python communications" +category = "main" +optional = false +python-versions = "*" +develop = false + +[package.source] +type = "git" +url = "https://github.com/m-labs/sipyco.git" +reference = "master" +resolved_reference = "58b0935f7ae47659abee5b5792fa594153328d6f" + +[[package]] +name = "six" +version = "1.16.0" +description = "Python 2 and 3 compatibility utilities" +category = "main" +optional = false +python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*" + +[[package]] +name = "statsmodels" +version = "0.13.2" +description = "Statistical computations and models for Python" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +numpy = ">=1.17" +packaging = ">=21.3" +pandas = ">=0.25" +patsy = ">=0.5.2" +scipy = ">=1.3" + +[package.extras] +build = ["cython (>=0.29.26)"] +develop = ["cython (>=0.29.26)"] +docs = ["sphinx", "nbconvert", "jupyter-client", "ipykernel", "matplotlib", "nbformat", "numpydoc", "pandas-datareader"] + +[[package]] +name = "toml" +version = "0.10.2" +description = "Python Library for Tom's Obvious, Minimal Language" +category = "dev" +optional = false +python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*" + +[[package]] +name = "unittest" +version = "0.0" +description = "" +category = "dev" +optional = false +python-versions = "*" + +[[package]] +name = "wcwidth" +version = "0.2.5" +description = "Measures the displayed width of unicode strings in a terminal" +category = "main" +optional = false +python-versions = "*" + +[[package]] +name = "yapf" +version = "0.32.0" +description = "A formatter for Python code." +category = "dev" +optional = false +python-versions = "*" + +[metadata] +lock-version = "1.1" +python-versions = "~3.8" +content-hash = "c03959ad89fddc6113be37243b6f75e9a887bc51e8064f9659095e7672c587ba" + +[metadata.files] +artiq = [] +cffi = [] +flake8 = [] +h5py = [] +llvmlite = [] +mccabe = [] +numpy = [] +oitg = [] +packaging = [] +pandas = [] +patsy = [] +prettytable = [] +pycodestyle = [] +pycparser = [] +pyflakes = [] +pygit2 = [] +pyparsing = [] +pyqt5 = [] +pyqt5-qt5 = [] +pyqt5-sip = [] +pyqtgraph = [] +python-dateutil = [] +python-levenshtein = [] +pythonparser = [] +pytz = [] +qasync = [] +regex = [] +scipy = [] +sipyco = [] +six = [] +statsmodels = [] +toml = [] +unittest = [] +wcwidth = [] +yapf = [] diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..9f49221f --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,42 @@ +[tool.poetry] +name = "ndscan" +version = "v0.2.0" +description = "N-dimensional scans for ARTIQ" +authors = ["David Nadlinger "] +license = "LGPLv3+" + +[tool.poetry.dependencies] +python = "~3.8" +#numpy = "^1.22.4" +#h5py = "^3.7.0" +#pyqtgraph = "^0.12.4" +#PyQt5 = "^5.15.6" +#PyQt5-Qt5 = "^5.15.2" +#PyQt5-sip = "^12.10.1" +#qasync = "^0.23.0" +#artiq = {git = "https://github.com/m-labs/artiq.git"} +#sipyco = {git = "https://github.com/m-labs/sipyco.git"} +#oitg = {git = "https://github.com/oxfordiontrapgroup/oitg.git"} + +[tool.poetry.dev-dependencies] +unittest = "^0.0" +yapf = "^0.32.0" +flake8 = "^4.0.1" +toml = "^0.10.2" + +[tool.poetry-dynamic-versioning] +enable = true +vcs = "git" +style = "pep440" +format-jinja = "{{ base }}.{{ distance }}+g{{ commit }}" +pattern = "^(?P\\d+(\\.\\d+)*)" + +[build-system] +requires = ["poetry-core>=1.0.0", "poetry-dynamic-versioning"] +build-backend = "poetry.core.masonry.api" + +[tool.poe.tasks] +fmt = "yapf -i -r examples ndscan test" +fmt-test = "yapf -d -r examples ndscan test" +flake = "flake8 examples ndscan test" +test = "python -m unittest discover -v test" From 5973bef4b17a8cd6cb348b251ba6ee52838ada96 Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 20:57:27 +0100 Subject: [PATCH 04/27] fix move to FitDescription, yapf, flake --- ndscan/experiment/default_analysis.py | 33 +++++++++++++-------------- ndscan/experiment/scan_generator.py | 1 - ndscan/plots/model/__init__.py | 2 +- 3 files changed, 17 insertions(+), 19 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 8b765e4b..742a650e 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -20,7 +20,6 @@ import logging from typing import Any, Callable, Dict, List, Iterable, Optional, Set, Tuple, Union import collections -from enum import Enum import dataclasses from .parameters import ParamHandle @@ -49,7 +48,8 @@ class FitDescription(): Attributes: fit_class_name: name of the fit class to use - fit_module: module that fit_class resides in (must be in the current python path) + fit_module: module that fit_class resides in (must be in the current python + path) data: maps fit data axis names (``"x"``, ``"y"``) to parameter handles or result channels that supply the respective data. param_bounds: dictionary of tuples containing the lower and upper bounds for @@ -160,7 +160,7 @@ def required_axes(self) -> Set[ParamHandle]: raise NotImplementedError def describe_online_analyses( - self, context: AnnotationContext + self, context: AnnotationContext ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: """Exceute analysis and serialise information about resulting annotations and online analyses to stringly typed metadata. @@ -245,7 +245,7 @@ def required_axes(self) -> Set[ParamHandle]: return self._required_axis_handles def describe_online_analyses( - self, context: AnnotationContext + self, context: AnnotationContext ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: "" return [], {} @@ -405,7 +405,7 @@ def required_axes(self) -> Set[ParamHandle]: return set(a for a in self.data.values() if isinstance(a, ParamHandle)) def describe_online_analyses( - self, context: AnnotationContext + self, context: AnnotationContext ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: "" # TODO: Generalise to higher-dimensional fits. @@ -454,17 +454,16 @@ def analysis_ref(key): analysis = { analysis_identifier: - FitDescription( - fit_class_name=self.fit_class_name, - fit_module=self.fit_module, - data={name: context.describe_coordinate(obj) - for name, obj in self.data.items() - }, - param_bounds=self.bounds, - fixed_params=self.constants, - initial_values=self.initial_values, - scale_factors=self.scale_factors - ) + FitDescription(fit_class_name=self.fit_class_name, + fit_module=self.fit_module, + data={ + name: context.describe_coordinate(obj) + for name, obj in self.data.items() + }, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + scale_factors=self.scale_factors) } return [a.describe(context) for a in annotations], analysis @@ -523,7 +522,7 @@ def required_axes(self) -> Set[ParamHandle]: return self._wrapped.required_axes() def describe_online_analyses( - self, context: AnnotationContext + self, context: AnnotationContext ) -> Tuple[List[Dict[str, Any]], Dict[str, FitDescription]]: return self._wrapped.describe_online_analyses(context) diff --git a/ndscan/experiment/scan_generator.py b/ndscan/experiment/scan_generator.py index 6682cf5c..952058f7 100644 --- a/ndscan/experiment/scan_generator.py +++ b/ndscan/experiment/scan_generator.py @@ -92,7 +92,6 @@ def __init__(self, def has_level(self, level: int) -> bool: "" - def num_points(limit): return np.floor(abs(self.centre - limit) / self.spacing) diff --git a/ndscan/plots/model/__init__.py b/ndscan/plots/model/__init__.py index 97ee2eaf..14f39e16 100644 --- a/ndscan/plots/model/__init__.py +++ b/ndscan/plots/model/__init__.py @@ -249,7 +249,7 @@ def _set_online_analyses(self, analysis_schemata: Dict[str, Dict[str, self._online_analyses = {} for name, schema in analysis_schemata.items(): - analysis = experiment.NamedFit(**schema) + analysis = experiment.FitDescription(**schema) self._online_analyses[name] = OnlineNamedFitAnalysis(analysis, self) # Rebind annotation schemata to new analysis data sources. From d51fd7f4cffb7079f93114e804a9a70c8d34e364 Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 22:38:51 +0100 Subject: [PATCH 05/27] add oitg fitting compatibility layer --- ndscan/experiment/default_analysis.py | 5 +- ndscan/fitting/fitting.py | 4 +- ndscan/fitting/oitg.py | 152 ++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 4 deletions(-) create mode 100644 ndscan/fitting/oitg.py diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 742a650e..575350c0 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -332,7 +332,8 @@ class OnlineFit(DefaultAnalysis): :param exported_results: Specifies fitted parameter values to export as analysis results. :param fit_module: python module containing the fit class. Will default to - `ndscan.fitting` in the future. + `ndscan.fitting` in the future. To use the oitg fitting functions, `fit_module` + should be set to the default `ndscan.fitting.oitg`. :param scale_factors: dictionary specifying scale factors for parameters. The parameter values are normalised by these values during fitting to help the curve fit (helps when fitting parameters with very large or very small values). At @@ -349,7 +350,7 @@ def __init__(self, initial_values: Optional[Dict[str, Any]] = None, bounds: Optional[Dict[str, Tuple[float, float]]] = None, exported_results: Optional[List[ExportedResult]] = None, - fit_module: str = 'oitg.fitting', + fit_module: str = 'ndscan.fitting.oitg', scale_factors: Optional[Dict[str, float]] = None): self.fit_class_name = fit_class diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index b2bbadf3..02913f5c 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -178,7 +178,7 @@ def func(x, params): raise NotImplementedError @staticmethod - def get_params(self) -> List[str]: + def get_params() -> List[str]: """Returns list of fit params""" raise NotImplementedError @@ -187,7 +187,7 @@ def get_free_params(self) -> List[str]: return self._free_params @staticmethod - def get_derived_params(self) -> List[str]: + def get_derived_params() -> List[str]: """Returns list of derived parameters""" return [] diff --git a/ndscan/fitting/oitg.py b/ndscan/fitting/oitg.py new file mode 100644 index 00000000..d1891912 --- /dev/null +++ b/ndscan/fitting/oitg.py @@ -0,0 +1,152 @@ +import numpy as np +from typing import Dict, List, Tuple +from .fitting import FitBase +import oitg.fitting + + +class OITGFit(FitBase): + _oitg_obj: oitg.fitting.FitBase + """ Compatibility shim for the OITG fitting functions """ + def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: + """ + Fit the dataset and return the fitted parameter values and uncertainties. + """ + if self._scale_factors != {}: + raise ValueError("OITG fitting functions do not accept user-specified" + "scale factors") + return super().fit() + + @classmethod + def func(cls, x, params): + """The fit function.""" + return cls._oitg_obj.fitting_function(x, params) + + @classmethod + def get_params(cls) -> List[str]: + """Returns list of fit params""" + return cls._oitg_obj.parameter_names + + @classmethod + def get_derived_params(cls) -> List[str]: + """Returns list of derived parameters""" + # 🙈🙈🙈🙈🙈🙈 a "not for poets" (TM) implementation + # the oitg code doesn't provide a list of derived parameter values so we feed it + # a series of meaningless parameter values chosen to avoid divide by zero errors + # and then see what it gives us back. + params = cls.get_params() + numbers = np.arange(1, len(params)) + p_dict = {key: number for key, number in zip(params, numbers)} + p_error_dict = dict(p_dict) + cls._oitg_obj.derived_parameter_function(p_dict, p_error_dict) + return list(set(p_dict.keys()) - set(params)) + + @classmethod + def get_default_bounds(cls) -> Dict[str, Tuple[float, float]]: + """ + Returns a dictionary mapping parameter names to a tuple of (upper, lower) bounds + """ + return cls._oitg_obj.parameter_bounds + + def _calculate_derived_params(self): + """ + Updates fit results with values and uncertainties for derived parameters. + """ + self._oitg_obj.derived_parameter_function(self._p, self._p_err) + + def _estimate_parameters(self) -> Dict[str, float]: + """ + Returns a dictionary of estimates for the parameter values for the current + dataset. + """ + initial_values = {} + self._oitg_obj.parameter_initialiser(self._x, self._y, initial_values) + return initial_values + + +class sinusoid(OITGFit): + _oitg_obj = oitg.fitting.sinusoid + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"pi_time": {"x": "t_pi"}} + + +class cos(OITGFit): + _oitg_obj = oitg.fitting.cos + + +class decaying_sinusoid(OITGFit): + _oitg_obj = oitg.fitting.decaying_sinusoid + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"pi_time": {"x": "t_max_transfer"}} + + +class detuned_square_pulse(OITGFit): + _oitg_obj = oitg.fitting.detuned_square_pulse + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"centre": {"x": "offset"}} + + +class exponential_decay(OITGFit): + _oitg_obj = oitg.fitting.exponential_decay + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"t_1_e": {"x": "t_1_e"}} + + +class gaussian(OITGFit): + _oitg_obj = oitg.fitting.gaussian + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"centre": {"x": "x0"}} + + +class line(OITGFit): + _oitg_obj = oitg.fitting.line + + +class lorentzian(OITGFit): + _oitg_obj = oitg.fitting.lorentzian + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"extremum": {"x": "x0"}} + + +class rabi_flop(OITGFit): + _oitg_obj = oitg.fitting.rabi_flop + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"pi_time": {"x": "t_pi"}} + + +class v_function(OITGFit): + _oitg_obj = oitg.fitting.v_function + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"centre": {"x": "x0"}} + + +class parabola(OITGFit): + _oitg_obj = oitg.fitting.shifted_parabola + + @staticmethod + def get_default_annotations(): + """ Returns a dictionary of default annotations for the fit. """ + return {"extremum": {"x": "position"}} From 7b6f536f2c9f73372ce32e47d19232c08f7e9ad5 Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 22:53:06 +0100 Subject: [PATCH 06/27] remove accidentally committed poetry files, change fitting imports to be consistent with the rest of ndscan --- ndscan/fitting/__init__.py | 8 +- ndscan/fitting/fitting.py | 2 + ndscan/fitting/sinusoid.py | 2 + poetry.lock | 431 ------------------------------------- pyproject.toml | 42 ---- 5 files changed, 10 insertions(+), 475 deletions(-) delete mode 100644 poetry.lock delete mode 100644 pyproject.toml diff --git a/ndscan/fitting/__init__.py b/ndscan/fitting/__init__.py index c8e9af68..7a50ec09 100644 --- a/ndscan/fitting/__init__.py +++ b/ndscan/fitting/__init__.py @@ -1,2 +1,6 @@ -from .fitting import FitBase # noqa F401 -from .sinusoid import Sinusoid # noqa F401 +from .fitting import * +from .sinusoid import * + +__all__ = [] +__all__.extend(fitting.__all__) +__all__.extend(sinusoid.__all__) diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 02913f5c..911251fc 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -4,6 +4,8 @@ from typing import Dict, List, Optional, Tuple, Union from scipy import optimize +__all__ = ['FitBase'] + # TODO: type annotations! (What's the best way of annotating np arrays?) # TODO: docstring formatting for sphinx diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index ba5de16f..b6ebb4b0 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -2,6 +2,8 @@ from typing import Dict, List, Tuple from ndscan import fitting +__all__ = ['Sinusoid'] + class Sinusoid(fitting.FitBase): """Sinusoid fit according to: diff --git a/poetry.lock b/poetry.lock deleted file mode 100644 index 34fc6f33..00000000 --- a/poetry.lock +++ /dev/null @@ -1,431 +0,0 @@ -[[package]] -name = "artiq" -version = "7.0b0" -description = "Advanced Real-Time Infrastructure for Quantum physics" -category = "main" -optional = false -python-versions = "*" -develop = false - -[package.dependencies] -h5py = "*" -llvmlite = "*" -numpy = "*" -prettytable = "*" -pygit2 = "*" -pyqtgraph = "*" -python-dateutil = "*" -python-Levenshtein = "*" -pythonparser = "*" -qasync = "*" -scipy = "*" - -[package.source] -type = "git" -url = "https://github.com/m-labs/artiq.git" -reference = "master" -resolved_reference = "f31279411ef7f198ee697290efe957544896ebdb" - -[[package]] -name = "cffi" -version = "1.15.1" -description = "Foreign Function Interface for Python calling C code." -category = "main" -optional = false -python-versions = "*" - -[package.dependencies] -pycparser = "*" - -[[package]] -name = "flake8" -version = "4.0.1" -description = "the modular source code checker: pep8 pyflakes and co" -category = "dev" -optional = false -python-versions = ">=3.6" - -[package.dependencies] -mccabe = ">=0.6.0,<0.7.0" -pycodestyle = ">=2.8.0,<2.9.0" -pyflakes = ">=2.4.0,<2.5.0" - -[[package]] -name = "h5py" -version = "3.7.0" -description = "Read and write HDF5 files from Python" -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -numpy = ">=1.14.5" - -[[package]] -name = "llvmlite" -version = "0.39.0" -description = "lightweight wrapper around basic LLVM functionality" -category = "main" -optional = false -python-versions = ">=3.7" - -[[package]] -name = "mccabe" -version = "0.6.1" -description = "McCabe checker, plugin for flake8" -category = "dev" -optional = false -python-versions = "*" - -[[package]] -name = "numpy" -version = "1.23.1" -description = "NumPy is the fundamental package for array computing with Python." -category = "main" -optional = false -python-versions = ">=3.8" - -[[package]] -name = "oitg" -version = "0.1" -description = "" -category = "main" -optional = false -python-versions = "*" -develop = false - -[package.dependencies] -statsmodels = "*" - -[package.source] -type = "git" -url = "https://github.com/oxfordiontrapgroup/oitg.git" -reference = "master" -resolved_reference = "937dfaa16a261e60f24c0165ddfc7818756e1683" - -[[package]] -name = "packaging" -version = "21.3" -description = "Core utilities for Python packages" -category = "main" -optional = false -python-versions = ">=3.6" - -[package.dependencies] -pyparsing = ">=2.0.2,<3.0.5 || >3.0.5" - -[[package]] -name = "pandas" -version = "1.4.3" -description = "Powerful data structures for data analysis, time series, and statistics" -category = "main" -optional = false -python-versions = ">=3.8" - -[package.dependencies] -numpy = [ - {version = ">=1.18.5", markers = "platform_machine != \"aarch64\" and platform_machine != \"arm64\" and python_version < \"3.10\""}, - {version = ">=1.19.2", markers = "platform_machine == \"aarch64\" and python_version < \"3.10\""}, - {version = ">=1.20.0", markers = "platform_machine == \"arm64\" and python_version < \"3.10\""}, -] -python-dateutil = ">=2.8.1" -pytz = ">=2020.1" - -[package.extras] -test = ["hypothesis (>=5.5.3)", "pytest (>=6.0)", "pytest-xdist (>=1.31)"] - -[[package]] -name = "patsy" -version = "0.5.2" -description = "A Python package for describing statistical models and for building design matrices." -category = "main" -optional = false -python-versions = "*" - -[package.dependencies] -numpy = ">=1.4" -six = "*" - -[package.extras] -test = ["scipy", "pytest-cov", "pytest"] - -[[package]] -name = "prettytable" -version = "3.3.0" -description = "A simple Python library for easily displaying tabular data in a visually appealing ASCII table format" -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -wcwidth = "*" - -[package.extras] -tests = ["pytest", "pytest-cov", "pytest-lazy-fixture"] - -[[package]] -name = "pycodestyle" -version = "2.8.0" -description = "Python style guide checker" -category = "dev" -optional = false -python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*" - -[[package]] -name = "pycparser" -version = "2.21" -description = "C parser in Python" -category = "main" -optional = false -python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" - -[[package]] -name = "pyflakes" -version = "2.4.0" -description = "passive checker of Python programs" -category = "dev" -optional = false -python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" - -[[package]] -name = "pygit2" -version = "1.10.0" -description = "Python bindings for libgit2." -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -cffi = ">=1.9.1" - -[[package]] -name = "pyparsing" -version = "3.0.9" -description = "pyparsing module - Classes and methods to define and execute parsing grammars" -category = "main" -optional = false -python-versions = ">=3.6.8" - -[package.extras] -diagrams = ["railroad-diagrams", "jinja2"] - -[[package]] -name = "pyqt5" -version = "5.15.7" -description = "Python bindings for the Qt cross platform application toolkit" -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -PyQt5-Qt5 = ">=5.15.0" -PyQt5-sip = ">=12.11,<13" - -[[package]] -name = "pyqt5-qt5" -version = "5.15.2" -description = "The subset of a Qt installation needed by PyQt5." -category = "main" -optional = false -python-versions = "*" - -[[package]] -name = "pyqt5-sip" -version = "12.11.0" -description = "The sip module support for PyQt5" -category = "main" -optional = false -python-versions = ">=3.7" - -[[package]] -name = "pyqtgraph" -version = "0.12.4" -description = "Scientific Graphics and GUI Library for Python" -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -numpy = ">=1.17.0" - -[[package]] -name = "python-dateutil" -version = "2.8.2" -description = "Extensions to the standard Python datetime module" -category = "main" -optional = false -python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7" - -[package.dependencies] -six = ">=1.5" - -[[package]] -name = "python-levenshtein" -version = "0.12.2" -description = "Python extension for computing string edit distances and similarities." -category = "main" -optional = false -python-versions = "*" - -[[package]] -name = "pythonparser" -version = "1.3" -description = "A Python parser intended for use in tooling" -category = "main" -optional = false -python-versions = "*" - -[package.dependencies] -regex = "*" - -[[package]] -name = "pytz" -version = "2022.1" -description = "World timezone definitions, modern and historical" -category = "main" -optional = false -python-versions = "*" - -[[package]] -name = "qasync" -version = "0.23.0" -description = "Implementation of the PEP 3156 Event-Loop with Qt." -category = "main" -optional = false -python-versions = "~=3.6" - -[[package]] -name = "regex" -version = "2022.7.25" -description = "Alternative regular expression module, to replace re." -category = "main" -optional = false -python-versions = ">=3.6" - -[[package]] -name = "scipy" -version = "1.9.0" -description = "SciPy: Scientific Library for Python" -category = "main" -optional = false -python-versions = ">=3.8,<3.12" - -[package.dependencies] -numpy = ">=1.18.5,<1.25.0" - -[[package]] -name = "sipyco" -version = "1.4" -description = "Simple Python communications" -category = "main" -optional = false -python-versions = "*" -develop = false - -[package.source] -type = "git" -url = "https://github.com/m-labs/sipyco.git" -reference = "master" -resolved_reference = "58b0935f7ae47659abee5b5792fa594153328d6f" - -[[package]] -name = "six" -version = "1.16.0" -description = "Python 2 and 3 compatibility utilities" -category = "main" -optional = false -python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*" - -[[package]] -name = "statsmodels" -version = "0.13.2" -description = "Statistical computations and models for Python" -category = "main" -optional = false -python-versions = ">=3.7" - -[package.dependencies] -numpy = ">=1.17" -packaging = ">=21.3" -pandas = ">=0.25" -patsy = ">=0.5.2" -scipy = ">=1.3" - -[package.extras] -build = ["cython (>=0.29.26)"] -develop = ["cython (>=0.29.26)"] -docs = ["sphinx", "nbconvert", "jupyter-client", "ipykernel", "matplotlib", "nbformat", "numpydoc", "pandas-datareader"] - -[[package]] -name = "toml" -version = "0.10.2" -description = "Python Library for Tom's Obvious, Minimal Language" -category = "dev" -optional = false -python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*" - -[[package]] -name = "unittest" -version = "0.0" -description = "" -category = "dev" -optional = false -python-versions = "*" - -[[package]] -name = "wcwidth" -version = "0.2.5" -description = "Measures the displayed width of unicode strings in a terminal" -category = "main" -optional = false -python-versions = "*" - -[[package]] -name = "yapf" -version = "0.32.0" -description = "A formatter for Python code." -category = "dev" -optional = false -python-versions = "*" - -[metadata] -lock-version = "1.1" -python-versions = "~3.8" -content-hash = "c03959ad89fddc6113be37243b6f75e9a887bc51e8064f9659095e7672c587ba" - -[metadata.files] -artiq = [] -cffi = [] -flake8 = [] -h5py = [] -llvmlite = [] -mccabe = [] -numpy = [] -oitg = [] -packaging = [] -pandas = [] -patsy = [] -prettytable = [] -pycodestyle = [] -pycparser = [] -pyflakes = [] -pygit2 = [] -pyparsing = [] -pyqt5 = [] -pyqt5-qt5 = [] -pyqt5-sip = [] -pyqtgraph = [] -python-dateutil = [] -python-levenshtein = [] -pythonparser = [] -pytz = [] -qasync = [] -regex = [] -scipy = [] -sipyco = [] -six = [] -statsmodels = [] -toml = [] -unittest = [] -wcwidth = [] -yapf = [] diff --git a/pyproject.toml b/pyproject.toml deleted file mode 100644 index 9f49221f..00000000 --- a/pyproject.toml +++ /dev/null @@ -1,42 +0,0 @@ -[tool.poetry] -name = "ndscan" -version = "v0.2.0" -description = "N-dimensional scans for ARTIQ" -authors = ["David Nadlinger "] -license = "LGPLv3+" - -[tool.poetry.dependencies] -python = "~3.8" -#numpy = "^1.22.4" -#h5py = "^3.7.0" -#pyqtgraph = "^0.12.4" -#PyQt5 = "^5.15.6" -#PyQt5-Qt5 = "^5.15.2" -#PyQt5-sip = "^12.10.1" -#qasync = "^0.23.0" -#artiq = {git = "https://github.com/m-labs/artiq.git"} -#sipyco = {git = "https://github.com/m-labs/sipyco.git"} -#oitg = {git = "https://github.com/oxfordiontrapgroup/oitg.git"} - -[tool.poetry.dev-dependencies] -unittest = "^0.0" -yapf = "^0.32.0" -flake8 = "^4.0.1" -toml = "^0.10.2" - -[tool.poetry-dynamic-versioning] -enable = true -vcs = "git" -style = "pep440" -format-jinja = "{{ base }}.{{ distance }}+g{{ commit }}" -pattern = "^(?P\\d+(\\.\\d+)*)" - -[build-system] -requires = ["poetry-core>=1.0.0", "poetry-dynamic-versioning"] -build-backend = "poetry.core.masonry.api" - -[tool.poe.tasks] -fmt = "yapf -i -r examples ndscan test" -fmt-test = "yapf -d -r examples ndscan test" -flake = "flake8 examples ndscan test" -test = "python -m unittest discover -v test" From 58386362b6eae2b729a338b69198d4a434bc46d8 Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 23:50:13 +0100 Subject: [PATCH 07/27] fitting.oitg: handle the fact that some fit functions do not provide a derived parameter function --- ndscan/fitting/oitg.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ndscan/fitting/oitg.py b/ndscan/fitting/oitg.py index d1891912..91b03c8a 100644 --- a/ndscan/fitting/oitg.py +++ b/ndscan/fitting/oitg.py @@ -33,6 +33,9 @@ def get_derived_params(cls) -> List[str]: # the oitg code doesn't provide a list of derived parameter values so we feed it # a series of meaningless parameter values chosen to avoid divide by zero errors # and then see what it gives us back. + if cls._oitg_obj.derived_parameter_function is None: + return [] + params = cls.get_params() numbers = np.arange(1, len(params)) p_dict = {key: number for key, number in zip(params, numbers)} From 2ce78303ad56c07c72677d5ddf5cb4d12bbcb4e8 Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 23:50:38 +0100 Subject: [PATCH 08/27] test_experiment_entrypoint: update after fitting changes --- test/test_experiment_entrypoint.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index 2837c81f..22a04201 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -114,7 +114,7 @@ def d(key): for i in range(len(timestamps)): prev = 0.0 if i == 0 else timestamps[i - 1] cur = timestamps[i] - self.assertGreater(cur, prev) + self.assertGreaterEqual(cur, prev) # Timestamps are in seconds, so this is a _very_ conservative bound on # running the test loop in-process (which will take much less than a # millisecond). @@ -141,31 +141,33 @@ def test_run_1d_scan(self): self.assertEqual(exp.fragment.num_host_cleanup_calls, 1) self.assertEqual(exp.fragment.num_device_cleanup_calls, 1) + pref = "fit_ndscan.fitting.oitg." curve_annotation = { "kind": "computed_curve", "parameters": { - "function_name": "lorentzian", + "fit_class_name": "lorentzian", + "fit_module": "ndscan.fitting.oitg", "associated_channels": ["channel_result"] }, "coordinates": {}, "data": { "a": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "a" }, "fwhm": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "fwhm" }, "x0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0" }, "y0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "y0" } @@ -178,14 +180,14 @@ def test_run_1d_scan(self): }, "coordinates": { "axis_0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0" } }, "data": { "axis_0_error": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0_error" } @@ -196,19 +198,21 @@ def test_run_1d_scan(self): self.assertEqual( json.loads(self.dataset_db.get("ndscan.rid_0.online_analyses")), { - "fit_lorentzian_channel_result": { - "constants": { + f"{pref}lorentzian_channel_result": { + "fixed_params": { "y0": 1.0 }, "data": { "y": "channel_result", "x": "axis_0" }, - "fit_type": "lorentzian", + "param_bounds": {} + "scale_factors": {} + "fit_class_name": "lorentzian", + "fit_module": "ndscan.fitting.oitg", "initial_values": { "fwhm": 2.0 }, - "kind": "named_fit" } }) From a7d19eeea221f1871b9e9d075d1bea9e5d080f7e Mon Sep 17 00:00:00 2001 From: hartytp Date: Tue, 9 Aug 2022 23:51:28 +0100 Subject: [PATCH 09/27] test_experiment_entrypoint: fix typo --- test/test_experiment_entrypoint.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index 22a04201..130e9a82 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -206,8 +206,8 @@ def test_run_1d_scan(self): "y": "channel_result", "x": "axis_0" }, - "param_bounds": {} - "scale_factors": {} + "param_bounds": {}, + "scale_factors": {}, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { From 6ed72716a9d87a772be920b7943c35563191b5e6 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 00:10:18 +0100 Subject: [PATCH 10/27] ndscan.fitting.oitg: handle the fact that the oitg code does not specify bounds for all parameters, handle (again) the lack of derived parameter functions in some cases --- ndscan/fitting/oitg.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/ndscan/fitting/oitg.py b/ndscan/fitting/oitg.py index 91b03c8a..ae32cc21 100644 --- a/ndscan/fitting/oitg.py +++ b/ndscan/fitting/oitg.py @@ -48,12 +48,19 @@ def get_default_bounds(cls) -> Dict[str, Tuple[float, float]]: """ Returns a dictionary mapping parameter names to a tuple of (upper, lower) bounds """ - return cls._oitg_obj.parameter_bounds + bounds = cls._oitg_obj.parameter_bounds + return { + param: bounds.get(param, (-np.inf, np.inf)) + for param in cls.get_params() + } def _calculate_derived_params(self): """ Updates fit results with values and uncertainties for derived parameters. """ + if self._oitg_obj.derived_parameter_function is None: + return + self._oitg_obj.derived_parameter_function(self._p, self._p_err) def _estimate_parameters(self) -> Dict[str, float]: From 050aaecbea357ae7ee556d4277c803db520aa1bf Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 00:11:39 +0100 Subject: [PATCH 11/27] test: update to new fitting code --- test/test_experiment_entrypoint.py | 9 +++++++-- test/test_experiment_subscan.py | 32 ++++++++++++++++++++---------- 2 files changed, 28 insertions(+), 13 deletions(-) diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index 130e9a82..13b6e60f 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -3,6 +3,7 @@ """ import json +import numpy as np from ndscan.experiment import * from ndscan.experiment.utils import is_kernel from ndscan.utils import PARAMS_ARG_KEY, SCHEMA_REVISION, SCHEMA_REVISION_KEY @@ -195,7 +196,6 @@ def test_run_1d_scan(self): } self.assertEqual(json.loads(self.dataset_db.get("ndscan.rid_0.annotations")), [curve_annotation, location_annotation]) - self.assertEqual( json.loads(self.dataset_db.get("ndscan.rid_0.online_analyses")), { f"{pref}lorentzian_channel_result": { @@ -206,7 +206,12 @@ def test_run_1d_scan(self): "y": "channel_result", "x": "axis_0" }, - "param_bounds": {}, + "param_bounds": { + "a": [-np.inf, np.inf], + "fwhm": [-np.inf, np.inf], + "x0": [-np.inf, np.inf], + "y0": [-np.inf, np.inf], + }, "scale_factors": {}, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index 780d7ac4..244b00bf 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -3,6 +3,7 @@ """ import json +import numpy as np from ndscan.experiment import * from fixtures import (AddOneFragment, ReboundAddOneFragment, AddOneCustomAnalysisFragment, TwoAnalysisFragment, @@ -51,31 +52,33 @@ def test_1d_result_channels(self): self.assertEqual(spec["fragment_fqn"], "fixtures.AddOneFragment") self.assertEqual(spec["seed"], 1234) + pref = "fit_ndscan.fitting.oitg." curve_annotation = { "kind": "computed_curve", "parameters": { - "function_name": "lorentzian", + "fit_class_name": "lorentzian", + "fit_module": "ndscan.fitting.oitg", "associated_channels": ["channel_result"] }, "coordinates": {}, "data": { "a": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "a" }, "fwhm": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "fwhm" }, "x0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0" }, "y0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "y0" } @@ -88,14 +91,14 @@ def test_1d_result_channels(self): }, "coordinates": { "axis_0": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0" } }, "data": { "axis_0_error": { - "analysis_name": "fit_lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_channel_result", "kind": "online_result", "result_key": "x0_error" } @@ -104,19 +107,26 @@ def test_1d_result_channels(self): self.assertEqual(spec["annotations"], [curve_annotation, location_annotation]) self.assertEqual( spec["online_analyses"], { - "fit_lorentzian_channel_result": { - "constants": { + f"{pref}lorentzian_channel_result": { + "fixed_params": { "y0": 1.0 }, + "param_bounds": { + "a": [-np.inf, np.inf], + "fwhm": [-np.inf, np.inf], + "x0": [-np.inf, np.inf], + "y0": [-np.inf, np.inf], + }, "data": { "y": "channel_result", "x": "axis_0" }, - "fit_type": "lorentzian", + "scale_factors": {}, + "fit_class_name": "lorentzian", + "fit_module": "ndscan.fitting.oitg", "initial_values": { "fwhm": 2.0 }, - "kind": "named_fit" } }) self.assertEqual( From 5e1484fffedeb960b9b3540c24ea0b35f6a2f3b1 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 10:26:56 +0100 Subject: [PATCH 12/27] fitting: add a kind parameter to FitDescription --- ndscan/experiment/default_analysis.py | 1 + test/test_experiment_entrypoint.py | 1 + test/test_experiment_subscan.py | 1 + 3 files changed, 3 insertions(+) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 575350c0..b4946016 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -71,6 +71,7 @@ class FitDescription(): fixed_params: Dict[str, float] initial_values: Dict[str, float] scale_factors: Dict[str, float] + kind: str = dataclasses.field(init=False, default="fit_description") class AnnotationContext: diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index 13b6e60f..d54b932c 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -199,6 +199,7 @@ def test_run_1d_scan(self): self.assertEqual( json.loads(self.dataset_db.get("ndscan.rid_0.online_analyses")), { f"{pref}lorentzian_channel_result": { + "kind": "fit_description", "fixed_params": { "y0": 1.0 }, diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index 244b00bf..02e5ac8f 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -108,6 +108,7 @@ def test_1d_result_channels(self): self.assertEqual( spec["online_analyses"], { f"{pref}lorentzian_channel_result": { + "kind": "fit_description", "fixed_params": { "y0": 1.0 }, From ca2daba74010b592834439c0a43fbca213cadb6e Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 11:38:19 +0100 Subject: [PATCH 13/27] FitDescription: add from_dict method --- ndscan/experiment/default_analysis.py | 10 ++++++++++ ndscan/plots/model/__init__.py | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index b4946016..65e5cfb2 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -73,6 +73,16 @@ class FitDescription(): scale_factors: Dict[str, float] kind: str = dataclasses.field(init=False, default="fit_description") + @classmethod + def from_dict(cls, data: dict): + kind = data.get('kind') + if kind != cls.kind: + raise ValueError("Attempt to construct FitDescription from dictionary with" + f"incorrect 'kind': {kind}") + data = dict(data) + del data['kind'] + return cls(**data) + class AnnotationContext: """Resolves entities in user-specified annotation schemata to stringly-typed diff --git a/ndscan/plots/model/__init__.py b/ndscan/plots/model/__init__.py index 14f39e16..245ac4ac 100644 --- a/ndscan/plots/model/__init__.py +++ b/ndscan/plots/model/__init__.py @@ -249,7 +249,7 @@ def _set_online_analyses(self, analysis_schemata: Dict[str, Dict[str, self._online_analyses = {} for name, schema in analysis_schemata.items(): - analysis = experiment.FitDescription(**schema) + analysis = experiment.FitDescription.from_dict(schema) self._online_analyses[name] = OnlineNamedFitAnalysis(analysis, self) # Rebind annotation schemata to new analysis data sources. From 5393a03f0e1ffe59659ef1b2b5a2ca0505882ea3 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 11:39:35 +0100 Subject: [PATCH 14/27] fitting: fix scale factors, sinusoid: add x0 parameter --- ndscan/fitting/fitting.py | 2 +- ndscan/fitting/oitg.py | 14 +++++++++++++- ndscan/fitting/sinusoid.py | 26 ++++++++++++++++++++------ 3 files changed, 34 insertions(+), 8 deletions(-) diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 911251fc..c70da3d3 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -126,7 +126,7 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: ) p *= scale_factors - p_cov *= scale_factors + p_cov *= scale_factors**2 self._p_cov = p_cov p_err = np.sqrt(np.diag(p_cov)) diff --git a/ndscan/fitting/oitg.py b/ndscan/fitting/oitg.py index ae32cc21..e84c1590 100644 --- a/ndscan/fitting/oitg.py +++ b/ndscan/fitting/oitg.py @@ -14,7 +14,19 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: if self._scale_factors != {}: raise ValueError("OITG fitting functions do not accept user-specified" "scale factors") - return super().fit() + + p, p_err = self._oitg_obj.fit(x=self._x, + y=self._y, + y_err=self._y_err, + constants=self._fixed_params, + initialise=self._initial_values, + calculate_residuals=False, + evaluate_function=False, + evaluate_x_limit=[None, None], + evaluate_n=1000) + self._p = p + self._p_err = p_err + return p, p_err @classmethod def func(cls, x, params): diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index b6ebb4b0..5263b037 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -7,19 +7,23 @@ class Sinusoid(fitting.FitBase): """Sinusoid fit according to: - y = a*sin(2*np.pi*f*x + phi) + offset + y = a*sin(2*np.pi*f*(x-x0) + phi) + offset TODO: t_dead, exp decay (but default to fixed at 0) """ @staticmethod def func(x, params): - return (params["a"] * np.sin(2 * np.pi * params["f"] * x + params["phi"]) + - params["offset"]) + w = 2 * np.pi * params["f"] + a = params["a"] + phi = params["phi"] + x0 = params["x0"] + offset = params["offset"] + return a * np.sin(w * (x - x0) + phi) + offset @staticmethod def get_params() -> List[str]: """Returns list of fit params""" - return ["a", "f", "phi", "offset", "t_dead"] + return ["a", "f", "phi", "offset", "t_dead", "x0"] @staticmethod def get_default_bounds() -> Dict[str, Tuple[float, float]]: @@ -31,7 +35,8 @@ def get_default_bounds() -> Dict[str, Tuple[float, float]]: "f": (0, np.inf), "phi": (0, 2 * np.pi), "offset": (0, 1), - "t_dead": (0, np.inf) + "t_dead": (0, np.inf), + "x0": (-np.inf, np.inf) } @staticmethod @@ -39,7 +44,7 @@ def get_default_fixed_params() -> Dict[str, float]: """Returns a dictionary mapping names of parameters which are not floated by default to their values. """ - return {'t_dead': 0} + return {'t_dead': 0, 'x0': 0} @staticmethod def get_derived_params() -> List[str]: @@ -82,6 +87,15 @@ def _estimate_parameters(self): param_guesses["phi"] -= 2 * np.pi param_guesses["t_dead"] = 0 + param_guesses["x0"] = 0 + + # allow the user to float x0 rather than phi + fixed = self._fixed_params + if "x0" not in fixed and "phi" in fixed: + w = 2 * np.pi * param_guesses["f"] + param_guesses["x0"] = (2 * np.pi - param_guesses["phi"]) / w + param_guesses["phi"] + return param_guesses def _calculate_derived_params(self): From 706254856baeda7a6acbe44732364c4c65aadb69 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 12:35:13 +0100 Subject: [PATCH 15/27] Specify an x-axis scale rather than manually specifing scale factors for each parameter. NB potentially contentious design decision: ParamHandles now store a reference to their Parameter. This allows one to get the parameter meta-information from the Handle. --- ndscan/experiment/default_analysis.py | 25 +++++++++---------------- ndscan/experiment/fragment.py | 7 ++++--- ndscan/experiment/parameters.py | 3 ++- ndscan/fitting/fitting.py | 18 +++++++++--------- ndscan/fitting/oitg.py | 4 ---- ndscan/fitting/sinusoid.py | 8 ++++++++ ndscan/plots/model/online_analysis.py | 2 +- test/test_experiment_entrypoint.py | 2 +- test/test_experiment_subscan.py | 2 +- 9 files changed, 35 insertions(+), 36 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 65e5cfb2..674fb370 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -58,10 +58,8 @@ class FitDescription(): parameters. If not specified, the defaults from the fit class are used. initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by the fit's parameter estimator. - scale_factors: dictionary specifying scale factors for parameters. The - parameter values are normalised by these values during fitting to help the - curve fit (helps when fitting parameters with very large or very small - values). + x_scale: x-axis scale factor used to normalise parameter values during fitting + to improve accuracy. """ fit_class_name: str @@ -70,7 +68,7 @@ class FitDescription(): param_bounds: Dict[str, Tuple[float, float]] fixed_params: Dict[str, float] initial_values: Dict[str, float] - scale_factors: Dict[str, float] + x_scale: float kind: str = dataclasses.field(init=False, default="fit_description") @classmethod @@ -345,12 +343,6 @@ class OnlineFit(DefaultAnalysis): :param fit_module: python module containing the fit class. Will default to `ndscan.fitting` in the future. To use the oitg fitting functions, `fit_module` should be set to the default `ndscan.fitting.oitg`. - :param scale_factors: dictionary specifying scale factors for parameters. The - parameter values are normalised by these values during fitting to help the curve - fit (helps when fitting parameters with very large or very small values). At - some point in the future we could consider choosing sensible default scale - factors, but it's not totally trivial to do that in a way that doesn't blow up - for parameters with initial values close to zero. """ def __init__(self, fit_class: str, @@ -361,8 +353,7 @@ def __init__(self, initial_values: Optional[Dict[str, Any]] = None, bounds: Optional[Dict[str, Tuple[float, float]]] = None, exported_results: Optional[List[ExportedResult]] = None, - fit_module: str = 'ndscan.fitting.oitg', - scale_factors: Optional[Dict[str, float]] = None): + fit_module: str = 'ndscan.fitting.oitg'): self.fit_class_name = fit_class self.fit_module = fit_module @@ -371,7 +362,9 @@ def __init__(self, self.analysis_identifier = analysis_identifier self.initial_values = initial_values or {} self.exported_results = exported_results or [] - self.scale_factors = scale_factors or {} + + x_param_handle = self.data['x'] + self.x_scale = x_param_handle.param.scale self._result_channels = [] for result in self.exported_results: @@ -475,7 +468,7 @@ def analysis_ref(key): param_bounds=self.bounds, fixed_params=self.constants, initial_values=self.initial_values, - scale_factors=self.scale_factors) + x_scale=self.x_scale) } return [a.describe(context) for a in annotations], analysis @@ -504,7 +497,7 @@ def execute( param_bounds=self.bounds, fixed_params=self.constants, initial_values=self.initial_values, - scale_factors=self.scale_factors) + x_scale=self.x_scale) for result in self.exported_results: p, p_err = getattr(fit, result.fit_parameter) channel = self._result_channels[result.result_name] diff --git a/ndscan/experiment/fragment.py b/ndscan/experiment/fragment.py index 87a8c9c8..e27284f2 100644 --- a/ndscan/experiment/fragment.py +++ b/ndscan/experiment/fragment.py @@ -317,9 +317,10 @@ def setattr_param(self, name: str, param_class: Type, description: str, *args, assert not hasattr(self, name), "Field '{}' already exists".format(name) fqn = self.fqn + "." + name - self._free_params[name] = param_class(fqn, description, *args, **kwargs) + param = param_class(fqn, description, *args, **kwargs) + self._free_params[name] = param - handle = param_class.HandleType(self, name) + handle = param_class.HandleType(self, name, param) setattr(self, name, handle) return handle @@ -368,7 +369,7 @@ def setattr_param_like(self, for k, v in kwargs.items(): setattr(new_param, k, v) self._free_params[name] = new_param - new_handle = new_param.HandleType(self, name) + new_handle = new_param.HandleType(self, name, new_param) setattr(self, name, new_handle) return new_handle diff --git a/ndscan/experiment/parameters.py b/ndscan/experiment/parameters.py index fcda1eff..0054365b 100644 --- a/ndscan/experiment/parameters.py +++ b/ndscan/experiment/parameters.py @@ -147,7 +147,8 @@ class ParamHandle: :param name: The name of the attribute in the owning fragment bound to this object. """ - def __init__(self, owner: Type["Fragment"], name: str): + def __init__(self, owner: Type["Fragment"], name: str, param: Any = None): + self.param = param self.owner = owner self.name = name assert name.isidentifier(), ("ParamHandle name should be the identifier it is " diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index c70da3d3..f8c07fc7 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -18,7 +18,7 @@ def __init__( param_bounds: Optional[Dict[str, Tuple[float, float]]] = None, fixed_params: Optional[Dict[str, float]] = None, initial_values: Optional[Dict[str, float]] = None, - scale_factors: Optional[Dict[str, float]] = None, + x_scale: float = 1, ): """ If `x` and `y` are provided we fit the provided dataset. `x`, `y` and `y_err` @@ -39,12 +39,8 @@ def __init__( :meth get_default_fixed_params: are used. :param initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by :meth estimate_parameters: - :param scale_factors: dictionary specifying scale factors for parameters. The - parameter values are normalised by these values during fitting to help the - curve fit (helps when fitting parameters with very large or very small - values). At some point in the future we could consider choosing sensible - default scale factors, but it's not totally trivial to do that in a way - that doesn't blow up for parameters with initial values close to zero. + :param x_scale: x-axis scale factor used to normalise parameter values + during fitting to improve accuracy. See :meth get_default_scale_factors:). """ self._x = x self._y = y @@ -71,8 +67,8 @@ def __init__( self._validate_param_names(initial_values, self._free_params) self._initial_values = initial_values or {} - self._validate_param_names(scale_factors, self._free_params) - self._scale_factors = scale_factors or {} + self._x_scale = x_scale + self._scale_factors = self.get_default_scale_factors() self._p = {} self._perr = {} @@ -168,6 +164,10 @@ def residuals(self): """Returns the fit residuals.""" return self._y - self._func_free(self._x, self._p) + def get_default_scale_factors(self) -> Dict[str, float]: + """ Returns a dictionary of default parameter scale factors. """ + return {param: 1 for param in self.get_params()} + def _calculate_derived_params(self): """ Updates fit results with values and uncertainties for derived parameters. diff --git a/ndscan/fitting/oitg.py b/ndscan/fitting/oitg.py index e84c1590..696d6942 100644 --- a/ndscan/fitting/oitg.py +++ b/ndscan/fitting/oitg.py @@ -11,10 +11,6 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: """ Fit the dataset and return the fitted parameter values and uncertainties. """ - if self._scale_factors != {}: - raise ValueError("OITG fitting functions do not accept user-specified" - "scale factors") - p, p_err = self._oitg_obj.fit(x=self._x, y=self._y, y_err=self._y_err, diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 5263b037..61ef5125 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -46,6 +46,14 @@ def get_default_fixed_params() -> Dict[str, float]: """ return {'t_dead': 0, 'x0': 0} + def get_default_scale_factors(self) -> Dict[str, float]: + """ Returns a dictionary of default parameter scale factors. """ + scales = super().get_default_scale_factors() + scales['f'] = 1 / self._x_scale + scales['x0'] = self._x_scale + scales['t_dead'] = self._x_scale + return scales + @staticmethod def get_derived_params() -> List[str]: """Returns a list of derived parameters. diff --git a/ndscan/plots/model/online_analysis.py b/ndscan/plots/model/online_analysis.py index 6b23d2b6..e2691260 100644 --- a/ndscan/plots/model/online_analysis.py +++ b/ndscan/plots/model/online_analysis.py @@ -35,7 +35,7 @@ def __init__(self, analysis: FitDescription, parent_model): self._fit_obj = fit_class(param_bounds=self._analysis.param_bounds, fixed_params=self._analysis.fixed_params, initial_values=self._analysis.initial_values, - scale_factors=self._analysis.scale_factors) + x_scale=self._analysis.x_scale) self._last_fit_params = None self._last_fit_errors = None diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index d54b932c..a2b5e922 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -213,7 +213,7 @@ def test_run_1d_scan(self): "x0": [-np.inf, np.inf], "y0": [-np.inf, np.inf], }, - "scale_factors": {}, + "x_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index 02e5ac8f..a6eac38f 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -122,7 +122,7 @@ def test_1d_result_channels(self): "y": "channel_result", "x": "axis_0" }, - "scale_factors": {}, + "x_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { From 1beaf6b67a399791b302cbd867ad8490381d85af Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 13:15:40 +0100 Subject: [PATCH 16/27] fitting: add y-axis scale factor as well --- ndscan/experiment/default_analysis.py | 13 +++++++++++-- ndscan/fitting/fitting.py | 5 +++++ ndscan/fitting/sinusoid.py | 4 +++- ndscan/plots/model/online_analysis.py | 3 ++- test/test_experiment_entrypoint.py | 1 + test/test_experiment_subscan.py | 1 + 6 files changed, 23 insertions(+), 4 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 674fb370..5d09828f 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -60,6 +60,8 @@ class FitDescription(): the fit. These override the values found by the fit's parameter estimator. x_scale: x-axis scale factor used to normalise parameter values during fitting to improve accuracy. + y_scale: y-axis scale factor used to normalise parameter values during fitting + to improve accuracy. """ fit_class_name: str @@ -69,6 +71,7 @@ class FitDescription(): fixed_params: Dict[str, float] initial_values: Dict[str, float] x_scale: float + y_scale: float kind: str = dataclasses.field(init=False, default="fit_description") @classmethod @@ -363,8 +366,12 @@ def __init__(self, self.initial_values = initial_values or {} self.exported_results = exported_results or [] + # TODO: how to deal with higher dimensional fits? + # TODO: are we happy assuming that x is always a parameter and y a result? x_param_handle = self.data['x'] + y_channel = self.data['y'] self.x_scale = x_param_handle.param.scale + self.y_scale = y_channel.scale self._result_channels = [] for result in self.exported_results: @@ -468,7 +475,8 @@ def analysis_ref(key): param_bounds=self.bounds, fixed_params=self.constants, initial_values=self.initial_values, - x_scale=self.x_scale) + x_scale=self.x_scale, + y_scale=self.y_scale) } return [a.describe(context) for a in annotations], analysis @@ -497,7 +505,8 @@ def execute( param_bounds=self.bounds, fixed_params=self.constants, initial_values=self.initial_values, - x_scale=self.x_scale) + x_scale=self.x_scale, + y_scale=self.y_scale) for result in self.exported_results: p, p_err = getattr(fit, result.fit_parameter) channel = self._result_channels[result.result_name] diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index f8c07fc7..486baf00 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -19,6 +19,7 @@ def __init__( fixed_params: Optional[Dict[str, float]] = None, initial_values: Optional[Dict[str, float]] = None, x_scale: float = 1, + y_scale: float = 1, ): """ If `x` and `y` are provided we fit the provided dataset. `x`, `y` and `y_err` @@ -41,6 +42,8 @@ def __init__( the fit. These override the values found by :meth estimate_parameters: :param x_scale: x-axis scale factor used to normalise parameter values during fitting to improve accuracy. See :meth get_default_scale_factors:). + :param y_scale: y-axis scale factor used to normalise parameter values + during fitting to improve accuracy. See :meth get_default_scale_factors:). """ self._x = x self._y = y @@ -68,7 +71,9 @@ def __init__( self._initial_values = initial_values or {} self._x_scale = x_scale + self._y_scale = y_scale self._scale_factors = self.get_default_scale_factors() + assert set(self._scale_factors) == set(self.get_params()) self._p = {} self._perr = {} diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 61ef5125..21f5428f 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -49,9 +49,11 @@ def get_default_fixed_params() -> Dict[str, float]: def get_default_scale_factors(self) -> Dict[str, float]: """ Returns a dictionary of default parameter scale factors. """ scales = super().get_default_scale_factors() + scales['a'] = self._y_scale scales['f'] = 1 / self._x_scale - scales['x0'] = self._x_scale + scales['offset'] = self._y_scale scales['t_dead'] = self._x_scale + scales['x0'] = self._x_scale return scales @staticmethod diff --git a/ndscan/plots/model/online_analysis.py b/ndscan/plots/model/online_analysis.py index e2691260..d6dfe166 100644 --- a/ndscan/plots/model/online_analysis.py +++ b/ndscan/plots/model/online_analysis.py @@ -35,7 +35,8 @@ def __init__(self, analysis: FitDescription, parent_model): self._fit_obj = fit_class(param_bounds=self._analysis.param_bounds, fixed_params=self._analysis.fixed_params, initial_values=self._analysis.initial_values, - x_scale=self._analysis.x_scale) + x_scale=self._analysis.x_scale, + y_scale=self._analysis.y_scale) self._last_fit_params = None self._last_fit_errors = None diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index a2b5e922..581d3387 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -214,6 +214,7 @@ def test_run_1d_scan(self): "y0": [-np.inf, np.inf], }, "x_scale": 1.0, + "y_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index a6eac38f..c01db217 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -123,6 +123,7 @@ def test_1d_result_channels(self): "x": "axis_0" }, "x_scale": 1.0, + "y_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { From 6fa674a686ddb53f55192ba205b144b7121179d7 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 13:17:32 +0100 Subject: [PATCH 17/27] fitting.sinusoid: docs --- ndscan/fitting/sinusoid.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 21f5428f..68dee92c 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -9,6 +9,10 @@ class Sinusoid(fitting.FitBase): """Sinusoid fit according to: y = a*sin(2*np.pi*f*(x-x0) + phi) + offset + By default, x0 is fixed at 0 and phi is floated. If phi is fixed and x0 floated, we + pick a default initial value for x0 that puts it in the first oscillation of the + sinusoid (equivalent to putting phi in the range [0, 2*pi]). + TODO: t_dead, exp decay (but default to fixed at 0) """ @staticmethod From 6e587a690fa6d2959b4e7ed5c94b5acb1efe41a9 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 13:35:15 +0100 Subject: [PATCH 18/27] Infer fit scale factors from dataset rather than from parameter/result scales. Revert "Specify an x-axis scale rather than manually specifing scale factors for" This reverts commit 706254856baeda7a6acbe44732364c4c65aadb69. --- ndscan/experiment/default_analysis.py | 51 ++++++++++----------------- ndscan/experiment/fragment.py | 7 ++-- ndscan/experiment/parameters.py | 3 +- ndscan/fitting/fitting.py | 32 +++++++++-------- ndscan/plots/model/online_analysis.py | 11 +++--- test/test_experiment_entrypoint.py | 2 -- test/test_experiment_subscan.py | 2 -- 7 files changed, 45 insertions(+), 63 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 5d09828f..9d4d813f 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -58,10 +58,6 @@ class FitDescription(): parameters. If not specified, the defaults from the fit class are used. initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by the fit's parameter estimator. - x_scale: x-axis scale factor used to normalise parameter values during fitting - to improve accuracy. - y_scale: y-axis scale factor used to normalise parameter values during fitting - to improve accuracy. """ fit_class_name: str @@ -70,8 +66,6 @@ class FitDescription(): param_bounds: Dict[str, Tuple[float, float]] fixed_params: Dict[str, float] initial_values: Dict[str, float] - x_scale: float - y_scale: float kind: str = dataclasses.field(init=False, default="fit_description") @classmethod @@ -366,13 +360,6 @@ def __init__(self, self.initial_values = initial_values or {} self.exported_results = exported_results or [] - # TODO: how to deal with higher dimensional fits? - # TODO: are we happy assuming that x is always a parameter and y a result? - x_param_handle = self.data['x'] - y_channel = self.data['y'] - self.x_scale = x_param_handle.param.scale - self.y_scale = y_channel.scale - self._result_channels = [] for result in self.exported_results: channel = result.result_type(path=result.result_name, **result.config) @@ -466,17 +453,17 @@ def analysis_ref(key): analysis = { analysis_identifier: - FitDescription(fit_class_name=self.fit_class_name, - fit_module=self.fit_module, - data={ - name: context.describe_coordinate(obj) - for name, obj in self.data.items() - }, - param_bounds=self.bounds, - fixed_params=self.constants, - initial_values=self.initial_values, - x_scale=self.x_scale, - y_scale=self.y_scale) + FitDescription( + fit_class_name=self.fit_class_name, + fit_module=self.fit_module, + data={ + name: context.describe_coordinate(obj) + for name, obj in self.data.items() + }, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + ) } return [a.describe(context) for a in annotations], analysis @@ -499,14 +486,14 @@ def execute( x = axis_data[self.data['x']._store.identity] y = result_data[self.data['y']] - fit = self.fit_klass(x=np.asarray(x), - y=np.asarray(y), - y_err=None, - param_bounds=self.bounds, - fixed_params=self.constants, - initial_values=self.initial_values, - x_scale=self.x_scale, - y_scale=self.y_scale) + fit = self.fit_klass( + x=np.asarray(x), + y=np.asarray(y), + y_err=None, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + ) for result in self.exported_results: p, p_err = getattr(fit, result.fit_parameter) channel = self._result_channels[result.result_name] diff --git a/ndscan/experiment/fragment.py b/ndscan/experiment/fragment.py index e27284f2..87a8c9c8 100644 --- a/ndscan/experiment/fragment.py +++ b/ndscan/experiment/fragment.py @@ -317,10 +317,9 @@ def setattr_param(self, name: str, param_class: Type, description: str, *args, assert not hasattr(self, name), "Field '{}' already exists".format(name) fqn = self.fqn + "." + name - param = param_class(fqn, description, *args, **kwargs) - self._free_params[name] = param + self._free_params[name] = param_class(fqn, description, *args, **kwargs) - handle = param_class.HandleType(self, name, param) + handle = param_class.HandleType(self, name) setattr(self, name, handle) return handle @@ -369,7 +368,7 @@ def setattr_param_like(self, for k, v in kwargs.items(): setattr(new_param, k, v) self._free_params[name] = new_param - new_handle = new_param.HandleType(self, name, new_param) + new_handle = new_param.HandleType(self, name) setattr(self, name, new_handle) return new_handle diff --git a/ndscan/experiment/parameters.py b/ndscan/experiment/parameters.py index 0054365b..fcda1eff 100644 --- a/ndscan/experiment/parameters.py +++ b/ndscan/experiment/parameters.py @@ -147,8 +147,7 @@ class ParamHandle: :param name: The name of the attribute in the owning fragment bound to this object. """ - def __init__(self, owner: Type["Fragment"], name: str, param: Any = None): - self.param = param + def __init__(self, owner: Type["Fragment"], name: str): self.owner = owner self.name = name assert name.isidentifier(), ("ParamHandle name should be the identifier it is " diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 486baf00..ca93871c 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -18,8 +18,6 @@ def __init__( param_bounds: Optional[Dict[str, Tuple[float, float]]] = None, fixed_params: Optional[Dict[str, float]] = None, initial_values: Optional[Dict[str, float]] = None, - x_scale: float = 1, - y_scale: float = 1, ): """ If `x` and `y` are provided we fit the provided dataset. `x`, `y` and `y_err` @@ -40,10 +38,6 @@ def __init__( :meth get_default_fixed_params: are used. :param initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by :meth estimate_parameters: - :param x_scale: x-axis scale factor used to normalise parameter values - during fitting to improve accuracy. See :meth get_default_scale_factors:). - :param y_scale: y-axis scale factor used to normalise parameter values - during fitting to improve accuracy. See :meth get_default_scale_factors:). """ self._x = x self._y = y @@ -70,11 +64,9 @@ def __init__( self._validate_param_names(initial_values, self._free_params) self._initial_values = initial_values or {} - self._x_scale = x_scale - self._y_scale = y_scale - self._scale_factors = self.get_default_scale_factors() - assert set(self._scale_factors) == set(self.get_params()) - + self._x_scale = None + self._y_scale = None + self._scale_factors = None self._p = {} self._perr = {} @@ -97,6 +89,11 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: y = self._y[valid_pts] y_err = None if self._y_err is None else self._y_err[valid_pts] + self._x_scale = np.max(np.abs(self._x)) + self._y_scale = np.max(np.abs(self._y)) + self._scale_factors = self.get_default_scale_factors() + assert set(self._scale_factors) == set(self.get_params()) + initial_values = self._estimate_parameters() initial_values.update(self._initial_values) @@ -169,10 +166,6 @@ def residuals(self): """Returns the fit residuals.""" return self._y - self._func_free(self._x, self._p) - def get_default_scale_factors(self) -> Dict[str, float]: - """ Returns a dictionary of default parameter scale factors. """ - return {param: 1 for param in self.get_params()} - def _calculate_derived_params(self): """ Updates fit results with values and uncertainties for derived parameters. @@ -252,6 +245,9 @@ def x(self, x): self._x = np.asarray(x) self._p = {} self._perr = {} + self._x_scale = None + self._y_scale = None + self._scale_factors = None @property def y(self): @@ -263,6 +259,9 @@ def y(self, y): self._y = np.asarray(y) self._p = {} self._perr = {} + self._x_scale = None + self._y_scale = None + self._scale_factors = None @property def y_err(self): @@ -274,3 +273,6 @@ def y_err(self, y_err=None): self._y_err = np.asarray(y_err) if y_err is not None else None self._p = {} self._perr = {} + self._x_scale = None + self._y_scale = None + self._scale_factors = None diff --git a/ndscan/plots/model/online_analysis.py b/ndscan/plots/model/online_analysis.py index d6dfe166..6a643d89 100644 --- a/ndscan/plots/model/online_analysis.py +++ b/ndscan/plots/model/online_analysis.py @@ -32,12 +32,11 @@ def __init__(self, analysis: FitDescription, parent_model): self._analysis = analysis fit_class = utils.import_class(analysis.fit_module, analysis.fit_class_name) - self._fit_obj = fit_class(param_bounds=self._analysis.param_bounds, - fixed_params=self._analysis.fixed_params, - initial_values=self._analysis.initial_values, - x_scale=self._analysis.x_scale, - y_scale=self._analysis.y_scale) - + self._fit_obj = fit_class( + param_bounds=self._analysis.param_bounds, + fixed_params=self._analysis.fixed_params, + initial_values=self._analysis.initial_values, + ) self._last_fit_params = None self._last_fit_errors = None diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index 581d3387..ffd8dc4b 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -213,8 +213,6 @@ def test_run_1d_scan(self): "x0": [-np.inf, np.inf], "y0": [-np.inf, np.inf], }, - "x_scale": 1.0, - "y_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index c01db217..389eb6b3 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -122,8 +122,6 @@ def test_1d_result_channels(self): "y": "channel_result", "x": "axis_0" }, - "x_scale": 1.0, - "y_scale": 1.0, "fit_class_name": "lorentzian", "fit_module": "ndscan.fitting.oitg", "initial_values": { From 54786787086951206d8f5c8577c7f1822dbcd420 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 13:42:27 +0100 Subject: [PATCH 19/27] fitting: add get_default_scale_factors back, sinusoid: give scale factors for all params explicitly --- ndscan/fitting/fitting.py | 4 ++++ ndscan/fitting/sinusoid.py | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index ca93871c..35b85617 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -210,6 +210,10 @@ def get_default_annotations(): """ Returns a dictionary of default annotations for the fit. """ return {} + def get_default_scale_factors(self) -> Dict[str, float]: + """ Returns a dictionary of default parameter scale factors. """ + return {param: 1 for param in self.get_params()} + def _estimate_parameters(self) -> Dict[str, float]: """ Returns a dictionary of estimates for the parameter values for the current diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 68dee92c..49520bdf 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -52,9 +52,10 @@ def get_default_fixed_params() -> Dict[str, float]: def get_default_scale_factors(self) -> Dict[str, float]: """ Returns a dictionary of default parameter scale factors. """ - scales = super().get_default_scale_factors() + scales = {} scales['a'] = self._y_scale scales['f'] = 1 / self._x_scale + scales['phi'] = 1 scales['offset'] = self._y_scale scales['t_dead'] = self._x_scale scales['x0'] = self._x_scale From 8fc29382fe9e503dcf802b521177e4168b78ebed Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 17:36:10 +0100 Subject: [PATCH 20/27] fitting: improve sinusoid parameter estimator, don't use abolute sigma without error bars --- ndscan/fitting/fitting.py | 2 +- ndscan/fitting/sinusoid.py | 79 ++++++++++++++++++++++++-------------- 2 files changed, 52 insertions(+), 29 deletions(-) diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 35b85617..0b51fab5 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -119,7 +119,7 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: ydata=y, p0=p0, sigma=y_err, - absolute_sigma=True, + absolute_sigma=y_err is not None, bounds=(lower_bounds, upper_bounds), ) diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 49520bdf..5b2485f9 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -1,4 +1,5 @@ import numpy as np +from scipy import signal, optimize from typing import Dict, List, Tuple from ndscan import fitting @@ -71,44 +72,66 @@ def get_derived_params() -> List[str]: return ["omega", 't_pi', 't_pi_2', 'phi_cosine', "min", "max", "period"] def _estimate_parameters(self): - # The OITG code uses a LS periodogram. I found that was giving poor phase - # estimates in some cases so here I hacked something with interpolation + and - # FFT. - # TODO: do something better! - sorted_idx = np.argsort(self.x) - x = self.x[sorted_idx] - y = self.y[sorted_idx] - - dx = np.diff(x).min() - x_axis = np.arange(x.min(), x.max(), dx) - y_interp = np.interp(x_axis, x, y) - N = len(x_axis) - - yf = np.fft.fft(y_interp)[0:N // 2] * 2.0 / N - freq = np.fft.fftfreq(N, dx)[:N // 2] - y0 = yf[0] - yf = yf[1:] - freq = freq[1:] - peak = np.argmax(np.abs(yf)) param_guesses = {} - param_guesses["a"] = np.abs(yf[peak]) - param_guesses["f"] = freq[peak] - param_guesses["offset"] = np.abs(y0) / 2 - param_guesses["phi"] = np.angle(yf[peak]) + np.pi / 2 - if param_guesses["phi"] < 0: - param_guesses["phi"] += 2 * np.pi - if param_guesses["phi"] > 2 * np.pi: - param_guesses["phi"] -= 2 * np.pi + sorted_inds = np.argsort(self._x) + x = self._x[sorted_inds] + y = self._y[sorted_inds] + y_err = None if self._y_err is None else self._y_err[sorted_inds] + param_guesses["offset"] = np.mean(y) param_guesses["t_dead"] = 0 param_guesses["x0"] = 0 + # Step 1: use a Lomb–Scargle Periodogram to estimate the frequency and amplitude + # of the signal + min_step = np.diff(x).min() + length = x.ptp() + + # Nyquist limit does not apply to irregularly spaced data + # We'll use it as a starting point anyway... + f_max = 0.5 / min_step + # relaxed Fourier limit + f_min = 0.25 / length + + omega_list = 2 * np.pi * np.linspace(f_min, f_max, int(f_max / f_min)) + pgram = signal.lombscargle(x, y, omega_list, precenter=True) + peak = np.argmax(np.abs(pgram)) + + param_guesses["a"] = np.sqrt(pgram[peak] * 4 / len(y)) + param_guesses["f"] = omega_list[peak] / (2 * np.pi) + + # Step 2: crude initial guess of the phase + def fit_fun(x, phi): + params = dict(param_guesses) + params["phi"] = phi + return self.func(x, params) + + def cost_fun(x, phi): + return np.sum(np.power(y - fit_fun(x, phi), 2)) + + phis = np.arange(1, 8) / (2 * np.pi) + costs = np.zeros_like(phis) + for idx, phi in np.ndenumerate(phis): + costs[idx] = cost_fun(x, phi) + param_guesses["phi"] = phis[np.argmin(costs)] + + # Step 3: single-parameter fit to find the phase (more robust than floating + # multiple parameters at once while the phase is off) + p, _ = optimize.curve_fit(f=fit_fun, + xdata=x, + ydata=y, + p0=param_guesses["phi"], + sigma=y_err, + absolute_sigma=True, + bounds=self.get_default_bounds()["phi"]) + param_guesses["phi"] = float(p) + # allow the user to float x0 rather than phi fixed = self._fixed_params if "x0" not in fixed and "phi" in fixed: w = 2 * np.pi * param_guesses["f"] - param_guesses["x0"] = (2 * np.pi - param_guesses["phi"]) / w + param_guesses["x0"] = param_guesses["phi"] / w param_guesses["phi"] return param_guesses From fb47dde2477ac67969db47982a244b30c23ccab3 Mon Sep 17 00:00:00 2001 From: hartytp Date: Wed, 10 Aug 2022 21:25:14 +0100 Subject: [PATCH 21/27] fitting.sinusoid: fix phase estimation and phase wrapping --- ndscan/fitting/sinusoid.py | 47 +++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 5b2485f9..fd77a7b7 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -110,7 +110,7 @@ def fit_fun(x, phi): def cost_fun(x, phi): return np.sum(np.power(y - fit_fun(x, phi), 2)) - phis = np.arange(1, 8) / (2 * np.pi) + phis = np.linspace(0, 2 * np.pi, 16) costs = np.zeros_like(phis) for idx, phi in np.ndenumerate(phis): costs[idx] = cost_fun(x, phi) @@ -124,15 +124,50 @@ def cost_fun(x, phi): p0=param_guesses["phi"], sigma=y_err, absolute_sigma=True, - bounds=self.get_default_bounds()["phi"]) + bounds=(0, 2 * np.pi)) param_guesses["phi"] = float(p) - # allow the user to float x0 rather than phi + # x0 and phi are equivalent parameters, but in some cases it's useful to have + # access to both. If exactly one of them is floated (the main anticipated + # use-case) we can easily pick a sensible initial value for the other. + # If both are floated there is no well-defined solution so we bug out + period = 1 / param_guesses["f"] + w = 2 * np.pi * param_guesses["f"] fixed = self._fixed_params + if "x0" not in fixed and "phi" in fixed: - w = 2 * np.pi * param_guesses["f"] - param_guesses["x0"] = param_guesses["phi"] / w - param_guesses["phi"] + d_phi = param_guesses['phi'] - self._fixed_params['phi'] + param_guesses["x0"] = -d_phi / w + + elif "phi" not in fixed and "x0" in fixed: + d_x0 = -self._fixed_params['x0'] + param_guesses["phi"] = -d_x0 * w + + elif "phi" not in fixed and "x0" not in fixed: + raise ValueError('At most one of "phi" and "x0" can be floated at a time!') + + # handle phase wrapping + if "phi" not in fixed: + phi_bounds = self._param_bounds["phi"] + if param_guesses["phi"] < phi_bounds[0]: + diff = phi_bounds[0] - param_guesses["phi"] + param_guesses["phi"] += ((diff // (2 * np.pi)) + 1) * 2 * np.pi + if param_guesses["phi"] > phi_bounds[1]: + diff = param_guesses["phi"] - phi_bounds[1] + param_guesses["phi"] -= ((diff // (2 * np.pi)) + 1) * 2 * np.pi + else: + param_guesses['phi'] = self._fixed_params['phi'] + + if "x0" not in fixed: + x0_bounds = self._param_bounds["x0"] + if param_guesses["x0"] < x0_bounds[0]: + diff = x0_bounds[0] - param_guesses["x0"] + param_guesses["x0"] += ((diff // period) + 1) * period + if param_guesses["x0"] > x0_bounds[1]: + diff = param_guesses["x0"] - x0_bounds[1] + param_guesses["x0"] -= ((diff // period) + 1) * period + else: + param_guesses['x0'] = self._fixed_params['x0'] return param_guesses From 7261abedae30e867f795eecb42db939050835951 Mon Sep 17 00:00:00 2001 From: hartytp Date: Thu, 11 Aug 2022 09:42:31 +0100 Subject: [PATCH 22/27] fitting.sinusoid: document parameters --- ndscan/fitting/sinusoid.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index fd77a7b7..614ec5d3 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -14,6 +14,22 @@ class Sinusoid(fitting.FitBase): pick a default initial value for x0 that puts it in the first oscillation of the sinusoid (equivalent to putting phi in the range [0, 2*pi]). + Fit parameters: + - a: amplitude of the sinusoid + - f: frequency + - phi: phase offset + - offset: y-axis offset + - t_dead: dead_time + - x0: x-axis offset + + Derived parameters: + - omega: frequency in angular units + - t_pi: pi-time (including dead time) + - t_pi_2: pi/2-time (including dead time) + - phi_cosine: cosine phase (phi + pi/2) + - min/max: min / max values of the sinusoid (including the offset) + - period: period of oscillation + TODO: t_dead, exp decay (but default to fixed at 0) """ @staticmethod From 15702b48853ba474e22a5d9e2ee4faf38485ca7d Mon Sep 17 00:00:00 2001 From: hartytp Date: Thu, 11 Aug 2022 16:12:08 +0100 Subject: [PATCH 23/27] fitting.sinusoid: finish initial implementation --- ndscan/fitting/sinusoid.py | 169 +++++++++++++++++++------------------ 1 file changed, 86 insertions(+), 83 deletions(-) diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index 614ec5d3..eec39154 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -7,44 +7,57 @@ class Sinusoid(fitting.FitBase): - """Sinusoid fit according to: - y = a*sin(2*np.pi*f*(x-x0) + phi) + offset - - By default, x0 is fixed at 0 and phi is floated. If phi is fixed and x0 floated, we - pick a default initial value for x0 that puts it in the first oscillation of the - sinusoid (equivalent to putting phi in the range [0, 2*pi]). - - Fit parameters: - - a: amplitude of the sinusoid - - f: frequency - - phi: phase offset + """Generalised sinusoid fit according to: + x > x_dead: y = Gamma*a*sin(omega*(x - x0 - x_dead) + phi0) + offset + x <= x_dead: y = a*sin(-x0*omega + phi0) + offset + where Gamma = exp(-1/tau * (x - x_dead)) + + Fit parameters (all floated by default unless stated otherwise): + - a: initial (x = x_dead) amplitude of the decaying sinusoid + - omega: angular frequency + - phi0: phase offset - offset: y-axis offset - - t_dead: dead_time - - x0: x-axis offset + - x_dead: "dead time" (fixed to 0 by default) + - x0: x-axis time offset (fixed to 0 by default) + - tau: decay time constant (fixed to np.inf by default) Derived parameters: - - omega: frequency in angular units - - t_pi: pi-time (including dead time) - - t_pi_2: pi/2-time (including dead time) - - phi_cosine: cosine phase (phi + pi/2) - - min/max: min / max values of the sinusoid (including the offset) + - f: frequency + - phi_cosine: cosine phase (phi0 + pi/2) + - min/max: min / max values of the undamped sinusoid (including the offset and + decay). - period: period of oscillation - TODO: t_dead, exp decay (but default to fixed at 0) + All phases are in radians, frequencies are in angular units. + + x0 and phi0 are equivalent parametrisations for the phase offset, but in some cases + it works out convenient to have access to both (e.g. one as a fixed offset, the + other floated). At most one of them may be floated at once. By default, x0 is fixed + at 0 and phi0 is floated. + + TODO: calculate peak values of the damped sinusoid as well as time that the peak + value occurs at. """ @staticmethod def func(x, params): - w = 2 * np.pi * params["f"] a = params["a"] - phi = params["phi"] - x0 = params["x0"] + w = params["omega"] + phi0 = params["phi0"] offset = params["offset"] - return a * np.sin(w * (x - x0) + phi) + offset + x_dead = params["x_dead"] + x0 = params["x0"] + tau = params["tau"] + + Gamma = np.exp(-1 / tau * (x - x_dead)) + y = Gamma * a * np.sin(w * (x - x0 - x_dead) + phi0) + offset + y0 = a * np.sin(-x0 * w + phi0) + offset + y[x <= x_dead] = y0 + return y @staticmethod def get_params() -> List[str]: """Returns list of fit params""" - return ["a", "f", "phi", "offset", "t_dead", "x0"] + return ["a", "omega", "phi0", "offset", "x_dead", "x0", "tau"] @staticmethod def get_default_bounds() -> Dict[str, Tuple[float, float]]: @@ -53,11 +66,12 @@ def get_default_bounds() -> Dict[str, Tuple[float, float]]: """ return { "a": (0, 1), - "f": (0, np.inf), - "phi": (0, 2 * np.pi), + "omega": (0, np.inf), + "phi0": (0, 2 * np.pi), "offset": (0, 1), - "t_dead": (0, np.inf), - "x0": (-np.inf, np.inf) + "x_dead": (0, np.inf), + "x0": (-np.inf, np.inf), + "tau": (0, np.inf) } @staticmethod @@ -65,27 +79,24 @@ def get_default_fixed_params() -> Dict[str, float]: """Returns a dictionary mapping names of parameters which are not floated by default to their values. """ - return {'t_dead': 0, 'x0': 0} + return {'x_dead': 0, 'x0': 0, 'tau': np.inf} def get_default_scale_factors(self) -> Dict[str, float]: """ Returns a dictionary of default parameter scale factors. """ scales = {} scales['a'] = self._y_scale - scales['f'] = 1 / self._x_scale - scales['phi'] = 1 + scales['omega'] = 1 / self._x_scale + scales['phi0'] = 1 scales['offset'] = self._y_scale - scales['t_dead'] = self._x_scale + scales['x_dead'] = self._x_scale scales['x0'] = self._x_scale + scales['tau'] = self._x_scale return scales @staticmethod def get_derived_params() -> List[str]: - """Returns a list of derived parameters. - - NB we define t_pi (t_pi_2) as the time taken for a pi-pulse (pi/2-pulse) - including deadtime. As a result, the 2*pi time is not 2*t_pi. - """ - return ["omega", 't_pi', 't_pi_2', 'phi_cosine', "min", "max", "period"] + """Returns a list of derived parameters. """ + return ["f", "phi_cosine", "min", "max", "period"] def _estimate_parameters(self): param_guesses = {} @@ -96,8 +107,9 @@ def _estimate_parameters(self): y_err = None if self._y_err is None else self._y_err[sorted_inds] param_guesses["offset"] = np.mean(y) - param_guesses["t_dead"] = 0 + param_guesses["x_dead"] = 0 param_guesses["x0"] = 0 + param_guesses["tau"] = max(x) # TODO: more clever estimate here # Step 1: use a Lomb–Scargle Periodogram to estimate the frequency and amplitude # of the signal @@ -115,64 +127,64 @@ def _estimate_parameters(self): peak = np.argmax(np.abs(pgram)) param_guesses["a"] = np.sqrt(pgram[peak] * 4 / len(y)) - param_guesses["f"] = omega_list[peak] / (2 * np.pi) + param_guesses["omega"] = omega_list[peak] # Step 2: crude initial guess of the phase - def fit_fun(x, phi): + def fit_fun(x, phi0): params = dict(param_guesses) - params["phi"] = phi + params["phi0"] = phi0 return self.func(x, params) - def cost_fun(x, phi): - return np.sum(np.power(y - fit_fun(x, phi), 2)) + def cost_fun(x, phi0): + return np.sum(np.power(y - fit_fun(x, phi0), 2)) phis = np.linspace(0, 2 * np.pi, 16) costs = np.zeros_like(phis) - for idx, phi in np.ndenumerate(phis): - costs[idx] = cost_fun(x, phi) - param_guesses["phi"] = phis[np.argmin(costs)] + for idx, phi0 in np.ndenumerate(phis): + costs[idx] = cost_fun(x, phi0) + param_guesses["phi0"] = phis[np.argmin(costs)] # Step 3: single-parameter fit to find the phase (more robust than floating # multiple parameters at once while the phase is off) p, _ = optimize.curve_fit(f=fit_fun, xdata=x, ydata=y, - p0=param_guesses["phi"], + p0=param_guesses["phi0"], sigma=y_err, absolute_sigma=True, bounds=(0, 2 * np.pi)) - param_guesses["phi"] = float(p) + param_guesses["phi0"] = float(p) - # x0 and phi are equivalent parameters, but in some cases it's useful to have + # x0 and phi0 are equivalent parameters, but in some cases it's useful to have # access to both. If exactly one of them is floated (the main anticipated # use-case) we can easily pick a sensible initial value for the other. # If both are floated there is no well-defined solution so we bug out - period = 1 / param_guesses["f"] - w = 2 * np.pi * param_guesses["f"] + period = 2 * np.pi / param_guesses["omega"] + w = param_guesses["omega"] fixed = self._fixed_params - if "x0" not in fixed and "phi" in fixed: - d_phi = param_guesses['phi'] - self._fixed_params['phi'] + if "x0" not in fixed and "phi0" in fixed: + d_phi = param_guesses['phi0'] - self._fixed_params['phi0'] param_guesses["x0"] = -d_phi / w - elif "phi" not in fixed and "x0" in fixed: + elif "phi0" not in fixed and "x0" in fixed: d_x0 = -self._fixed_params['x0'] - param_guesses["phi"] = -d_x0 * w + param_guesses["phi0"] = -d_x0 * w - elif "phi" not in fixed and "x0" not in fixed: - raise ValueError('At most one of "phi" and "x0" can be floated at a time!') + elif "phi0" not in fixed and "x0" not in fixed: + raise ValueError('At most one of "phi0" and "x0" can be floated at a time!') # handle phase wrapping - if "phi" not in fixed: - phi_bounds = self._param_bounds["phi"] - if param_guesses["phi"] < phi_bounds[0]: - diff = phi_bounds[0] - param_guesses["phi"] - param_guesses["phi"] += ((diff // (2 * np.pi)) + 1) * 2 * np.pi - if param_guesses["phi"] > phi_bounds[1]: - diff = param_guesses["phi"] - phi_bounds[1] - param_guesses["phi"] -= ((diff // (2 * np.pi)) + 1) * 2 * np.pi + if "phi0" not in fixed: + phi_bounds = self._param_bounds["phi0"] + if param_guesses["phi0"] < phi_bounds[0]: + diff = phi_bounds[0] - param_guesses["phi0"] + param_guesses["phi0"] += ((diff // (2 * np.pi)) + 1) * 2 * np.pi + if param_guesses["phi0"] > phi_bounds[1]: + diff = param_guesses["phi0"] - phi_bounds[1] + param_guesses["phi0"] -= ((diff // (2 * np.pi)) + 1) * 2 * np.pi else: - param_guesses['phi'] = self._fixed_params['phi'] + param_guesses['phi0'] = self._fixed_params['phi0'] if "x0" not in fixed: x0_bounds = self._param_bounds["x0"] @@ -188,25 +200,16 @@ def cost_fun(x, phi): return param_guesses def _calculate_derived_params(self): - self._p["omega"] = 2 * np.pi * self._p["f"] - self._p['t_pi'] = self._p["t_dead"] + np.pi / self._p["omega"] - self._p["t_pi_2"] = self._p["t_dead"] + np.pi / 2 / self._p["omega"] - self._p["phi_cosine"] = self._p["phi"] + np.pi / 2 + self._p["f"] = self._p["omega"] / (2 * np.pi) + self._p["phi_cosine"] = self._p["phi0"] + np.pi / 2 self._p["min"] = self._p["offset"] - np.abs(self._p["a"]) self._p["max"] = self._p["offset"] + np.abs(self._p["a"]) self._p["period"] = 2 * np.pi / self._p["omega"] # TODO: consider covariances - self._p_err["omega"] = self._p_err["f"] * 2 * np.pi - self._p_err["t_pi"] = np.sqrt(self._p_err["t_dead"]**2 + - (np.pi / self._p["omega"] * - (self._p_err["omega"] / self._p["omega"]))**2) - - self._p_err["t_pi/2"] = np.sqrt(self._p_err["t_dead"]**2 + - (np.pi / 2 / self._p["omega"] * - (self._p_err["omega"] / self._p["omega"]))**2) + self._p_err["f"] = self._p_err["omega"] / (2 * np.pi) - self._p_err["phi_cosine"] = self._p_err["phi"] + self._p_err["phi_cosine"] = self._p_err["phi0"] self._p_err["min"] = np.sqrt(self._p_err["offset"]**2 + self._p_err["a"]**2) self._p_err["max"] = np.sqrt(self._p_err["offset"]**2 + self._p_err["a"]**2) @@ -222,7 +225,7 @@ def _calculate_derived_params(self): params = { "a": np.random.uniform(low=0.25, high=1, size=1), "f": np.random.uniform(low=0.5, high=3, size=1), - "phi": np.random.uniform(low=0, high=2 * np.pi, size=1), + "phi0": np.random.uniform(low=0, high=2 * np.pi, size=1), "offset": np.random.uniform(low=0.25, high=0.75, size=1), } y = Sinusoid.func(x, params) + np.random.normal(0, 0.05, x.shape) @@ -232,8 +235,8 @@ def _calculate_derived_params(self): print(f"Amplitude error is {100*(1-fit.a[0]/params['a'])}%") if not np.isclose(fit.f[0], params["f"], rtol=5e-2): print(f"Frequency error is {100*(1-fit.f[0]/params['f'])}%") - if not np.isclose(fit.phi[0], params["phi"], rtol=5e-2): - print(f"Phase error is {100*(1-fit.phi[0]/params['phi'])}%") + if not np.isclose(fit.phi0[0], params["phi0"], rtol=5e-2): + print(f"Phase error is {100*(1-fit.phi0[0]/params['phi0'])}%") if not np.isclose(fit.offset[0], params["offset"], rtol=5e-2): print(f"Offset error is {100*(1-fit.offset[0]/params['offset'])}%") From 95481347381e49032e9f41e99cf3db57c5daf159 Mon Sep 17 00:00:00 2001 From: hartytp Date: Thu, 11 Aug 2022 22:49:01 +0100 Subject: [PATCH 24/27] fitting: clip initial values to lie within bounds --- ndscan/fitting/fitting.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 0b51fab5..499be923 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -97,6 +97,11 @@ def fit(self) -> Tuple[Dict[str, float], Dict[str, float]]: initial_values = self._estimate_parameters() initial_values.update(self._initial_values) + for param in self.get_free_params(): + param_bounds = self._param_bounds[param] + initial_values[param] = np.clip(initial_values[param], param_bounds[0], + param_bounds[1]) + lower_bounds = [self._param_bounds[param][0] for param in self._free_params] upper_bounds = [self._param_bounds[param][1] for param in self._free_params] From 03c9b367a15d869dafc03d8d90a6bbdc5bad1ac9 Mon Sep 17 00:00:00 2001 From: hartytp Date: Thu, 11 Aug 2022 22:50:43 +0100 Subject: [PATCH 25/27] fitting: add rabi flop fit (untested) --- ndscan/fitting/rabi_flop.py | 179 ++++++++++++++++++++++++++++++++++++ ndscan/fitting/sinusoid.py | 36 +++++--- 2 files changed, 203 insertions(+), 12 deletions(-) create mode 100644 ndscan/fitting/rabi_flop.py diff --git a/ndscan/fitting/rabi_flop.py b/ndscan/fitting/rabi_flop.py new file mode 100644 index 00000000..77a4885d --- /dev/null +++ b/ndscan/fitting/rabi_flop.py @@ -0,0 +1,179 @@ +import numpy as np +from typing import Dict, List, Tuple +from . import FitBase, Sinusoid + + +class RabiFlop(FitBase): + """ Time-/frequency-domain exponentially damped Rabi flop fit according to: + y = P0 * y0 + P1 * y1 + y1 = Gamma * (contrast * (omega * t / 2 * sinc(W*t/2))^2 + P_lower - c) + c + Gamma = exp(-t/tau) + w = sqrt(omega^2 + delta^2) + contrast = P_upper - P_lower + c = 0.5 * (P_upper + P_lower) + P0 = 1 - P1 + y0 = 1 - y1 + + For frequency scans we use: t = t_pulse - t_dead, delta = x - detuning_offset + For time scans we use: {t = x - t_dead for x > t_dead, x = 0 for t <= t_dead}, + delta = detuning_offset + + This class is not intended to be used directly. See RabiFlopTime and + RabiFlopFrequency. + + Fit parameters (all floated by default unless stated otherwise): + - P1: initial upper-state population (fixed to 1 by default) + - P_upper: upper readout level (fixed to 1 by default) + - P_lower: lower readout level (fixed to 0 by default) + - omega: Rabi frequency + - t_pulse: pulse duration (frequency scans only) + - t_dead: dead_time (fixed to 0 by default) + - detuning_offset: detuning offset + - tau: decay time constant (fixed to np.inf by default) + + Derived parameters: + - t_pi: pi-time including dead-time (so t_2pi != 2*t_pi), is not the time for + maximum population transfer for finite tau (TODO: add that as a derived + parameter!) + - t_pi_2: pi/2-time including dead-time (so t_pi != 2*t_pi_2) + - TODO: do we want pulse area error, etc? + + All phases are in radians, detunings are in angular units. + """ + @staticmethod + def func(x, params): + raise NotImplementedError + + @staticmethod + def rabi_func(delta, t, params): + P1 = params["P1"] + P_upper = params["P_upper"] + P_lower = params["P_lower"] + omega = params["omega"] + tau = params["tau"] + + P0 = 1 - P1 + Gamma = np.exp(-t / tau) + W = np.sqrt(np.pow(omega, 2), np.pow(delta, 2)) + contrast = P_upper - P_lower + c = 0.5 * (P_upper + P_lower) + + # NB np.sinc(x) = sin(pi*x)/(pi*x) + y1 = Gamma * (contrast * np.pow(omega * t / 2 * np.sinc(W * t / + (2 * np.pi)), 2) + + P_lower - c) + c + y1[t < 0] = P1 * P_upper + y0 = 1 - y1 + + return P0 * y0 + P1 * y1 + + @staticmethod + def get_params() -> List[str]: + """Returns list of fit params""" + return ["P1", "P_upper", "P_lower", "omega", "t_dead", "tau", "detuning_offset"] + + @classmethod + def get_default_bounds(cls) -> Dict[str, Tuple[float, float]]: + """ + Returns a dictionary mapping parameter names to a tuple of (upper, lower) bounds + """ + bounds = { + "P1": (0, 1), + "P_upper": (0, 1), + "P_lower": (0, 1), + "omega": (0, np.inf), + "t_pulse": (0, np.inf), + "t_dead": (0, np.inf), + "detuning_offset": (-np.inf, np.inf), + "tau": (0, np.inf) + } + return {param: bounds[param] for param in cls.get_params()} + + @classmethod + def get_default_fixed_params(cls) -> Dict[str, float]: + """Returns a dictionary mapping names of parameters which are not floated by + default to their values. + """ + fixed = { + "P1": 0, + "P_upper": 1, + "P_lower": 0, + "t_dead": 0, + "tau": np.inf, + } + return {param: fixed[param] for param in cls.get_params()} + + def get_default_scale_factors(self) -> Dict[str, float]: + """ Returns a dictionary of default parameter scale factors. """ + scales = { + "P1": 1, + "P_upper": 1, + "P_lower": 1, + "omega": 1 / self._x_scale, + "t_pulse": self._x_scale, + "t_dead": self._x_scale, + "detuning_offset": 1 / self._x_scale, + "tau": self._x_scale + } + return {param: scales[param] for param in self.get_params()} + + @staticmethod + def get_derived_params() -> List[str]: + """Returns a list of derived parameters. + + NB we define t_pi (t_pi_2) as the time taken for a pi-pulse (pi/2-pulse) + including deadtime. As a result, the 2*pi time is not 2*t_pi. + """ + return ["t_pi", "t_pi_2"] + + def _estimate_parameters(self): + param_guesses = self._sinusoid_estimator( + x=self._x, + y=self._y, + y_err=self._y_err, + fixed_params=Sinusoid.get_default_fixed_params(), + param_bounds=Sinusoid.get_default_bounds(), + ) + + param_guesses['P1'] = 1 # TODO: estimate from sine phase? + param_guesses['P_upper'] = param_guesses['offset'] + np.abs(param_guesses['a']) + param_guesses['P_lower'] = param_guesses['offset'] + return param_guesses + + +class RabiFlopTime(RabiFlop): + """ Time-domain Rabi flops. See RabiFlop for details. """ + @staticmethod + def func(x, params): + t = x - params["t_dead"] + delta = params["detuning_offset"] + return RabiFlop.rabi_func(delta=delta, t=t, params=params) + + def _estimate_parameters(self): + param_guesses = super()._estimate_parameters() + param_guesses['detuning_offset'] = 0 + return {param: param_guesses[param] for param in self.get_params()} + + +class RabiFlopFrequency(RabiFlop): + """ Frequency-domain Rabi flops. See RabiFlop for details. """ + @staticmethod + def func(x, params): + t = params["t_pulse"] - params["t_dead"] + delta = x - params["detuning_offset"] + return RabiFlop.rabi_func(delta=delta, t=t, params=params) + + @staticmethod + def get_params() -> List[str]: + """Returns list of fit params""" + return super().get_params() + ["t_pulse"] + + def _estimate_parameters(self): + param_guesses = super()._estimate_parameters() + + # sqrt(3) factor derived from assuming pi pulse + param_guesses['t_pulse'] = param_guesses['omega'] / np.sqrt(3) + param_guesses['omega'] = np.pi / param_guesses["t_pulse"] + param_guesses['detuning_offset'] = param_guesses['x_peak'] + + return {param: param_guesses[param] for param in self.get_params()} diff --git a/ndscan/fitting/sinusoid.py b/ndscan/fitting/sinusoid.py index eec39154..24811c59 100644 --- a/ndscan/fitting/sinusoid.py +++ b/ndscan/fitting/sinusoid.py @@ -99,12 +99,23 @@ def get_derived_params() -> List[str]: return ["f", "phi_cosine", "min", "max", "period"] def _estimate_parameters(self): + param_guesses = self._sinusoid_estimator( + x=self._x, + y=self._y, + y_err=self._y_err, + fixed_params=self._fixed_params, + param_bounds=self._param_bounds, + ) + return {param: param_guesses[param] for param in self.get_params()} + + @staticmethod + def _sinusoid_estimator(x, y, y_err, fixed_params, param_bounds): param_guesses = {} - sorted_inds = np.argsort(self._x) - x = self._x[sorted_inds] - y = self._y[sorted_inds] - y_err = None if self._y_err is None else self._y_err[sorted_inds] + sorted_inds = np.argsort(x) + x = x[sorted_inds] + y = y[sorted_inds] + y_err = None if y_err is None else y_err[sorted_inds] param_guesses["offset"] = np.mean(y) param_guesses["x_dead"] = 0 @@ -128,12 +139,13 @@ def _estimate_parameters(self): param_guesses["a"] = np.sqrt(pgram[peak] * 4 / len(y)) param_guesses["omega"] = omega_list[peak] + param_guesses["x_peak"] = x[peak] # used for frequency fits # Step 2: crude initial guess of the phase def fit_fun(x, phi0): params = dict(param_guesses) params["phi0"] = phi0 - return self.func(x, params) + return Sinusoid.func(x, params) def cost_fun(x, phi0): return np.sum(np.power(y - fit_fun(x, phi0), 2)) @@ -161,14 +173,14 @@ def cost_fun(x, phi0): # If both are floated there is no well-defined solution so we bug out period = 2 * np.pi / param_guesses["omega"] w = param_guesses["omega"] - fixed = self._fixed_params + fixed = fixed_params if "x0" not in fixed and "phi0" in fixed: - d_phi = param_guesses['phi0'] - self._fixed_params['phi0'] + d_phi = param_guesses['phi0'] - fixed_params['phi0'] param_guesses["x0"] = -d_phi / w elif "phi0" not in fixed and "x0" in fixed: - d_x0 = -self._fixed_params['x0'] + d_x0 = -fixed_params['x0'] param_guesses["phi0"] = -d_x0 * w elif "phi0" not in fixed and "x0" not in fixed: @@ -176,7 +188,7 @@ def cost_fun(x, phi0): # handle phase wrapping if "phi0" not in fixed: - phi_bounds = self._param_bounds["phi0"] + phi_bounds = param_bounds["phi0"] if param_guesses["phi0"] < phi_bounds[0]: diff = phi_bounds[0] - param_guesses["phi0"] param_guesses["phi0"] += ((diff // (2 * np.pi)) + 1) * 2 * np.pi @@ -184,10 +196,10 @@ def cost_fun(x, phi0): diff = param_guesses["phi0"] - phi_bounds[1] param_guesses["phi0"] -= ((diff // (2 * np.pi)) + 1) * 2 * np.pi else: - param_guesses['phi0'] = self._fixed_params['phi0'] + param_guesses['phi0'] = fixed_params['phi0'] if "x0" not in fixed: - x0_bounds = self._param_bounds["x0"] + x0_bounds = param_bounds["x0"] if param_guesses["x0"] < x0_bounds[0]: diff = x0_bounds[0] - param_guesses["x0"] param_guesses["x0"] += ((diff // period) + 1) * period @@ -195,7 +207,7 @@ def cost_fun(x, phi0): diff = param_guesses["x0"] - x0_bounds[1] param_guesses["x0"] -= ((diff // period) + 1) * period else: - param_guesses['x0'] = self._fixed_params['x0'] + param_guesses['x0'] = fixed_params['x0'] return param_guesses From f3d9222cf842409c11587cba8b98996c90187d8f Mon Sep 17 00:00:00 2001 From: hartytp Date: Fri, 12 Aug 2022 01:46:17 +0100 Subject: [PATCH 26/27] fitting: add chi-squared results, bug fixes --- ndscan/experiment/default_analysis.py | 75 ++++++++++++++++----------- ndscan/fitting/fitting.py | 33 ++++++++++-- 2 files changed, 75 insertions(+), 33 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 9d4d813f..93f1eae4 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -16,14 +16,13 @@ Both can produce annotations; particular values or plot locations highlighted in the user interface. """ -import numpy as np import logging from typing import Any, Callable, Dict, List, Iterable, Optional, Set, Tuple, Union import collections import dataclasses from .parameters import ParamHandle -from .result_channels import ResultChannel, FloatChannel +from .result_channels import ResultChannel, FloatChannel, OpaqueChannel from .. import utils __all__ = [ @@ -58,6 +57,9 @@ class FitDescription(): parameters. If not specified, the defaults from the fit class are used. initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by the fit's parameter estimator. + export_chi2: If True and y-axis error is provided we export the fit Chi-Squared + and p-value as result channels named ``{analysis_identifier}_chi2`` and + ``{analysis_identifier}_p``. """ fit_class_name: str @@ -66,6 +68,7 @@ class FitDescription(): param_bounds: Dict[str, Tuple[float, float]] fixed_params: Dict[str, float] initial_values: Dict[str, float] + export_chi2: bool = True kind: str = dataclasses.field(init=False, default="fit_description") @classmethod @@ -327,8 +330,8 @@ class OnlineFit(DefaultAnalysis): dictionaries mapping coordinate names to fit result names. If ``None``, the defaults provided by the fit function will be used. :param analysis_identifier: Optional explicit name to use for online analysis. - Defaults to ``fit_``, but can be set explicitly to allow more than one - fit of a given type at a time. + Defaults to ``fit_{fit_module}.{fit_class}``, but can be set explicitly to allow + more than one fit of a given type at a time. :param constants: dictionary specifying constant values for any non-floated parameters. If not specified, the defaults from the fit class are used. :param initial_values: dictionary specifying initial parameter values to use in @@ -340,6 +343,9 @@ class OnlineFit(DefaultAnalysis): :param fit_module: python module containing the fit class. Will default to `ndscan.fitting` in the future. To use the oitg fitting functions, `fit_module` should be set to the default `ndscan.fitting.oitg`. + :param export_chi2: If True and y-axis error is provided we export the fit + Chi-Squared and p-value as result channels named + ``{analysis_identifier}_chi2`` and ``{analysis_identifier}_p``. """ def __init__(self, fit_class: str, @@ -350,8 +356,13 @@ def __init__(self, initial_values: Optional[Dict[str, Any]] = None, bounds: Optional[Dict[str, Tuple[float, float]]] = None, exported_results: Optional[List[ExportedResult]] = None, - fit_module: str = 'ndscan.fitting.oitg'): + fit_module: str = 'ndscan.fitting.oitg', + export_chi2: bool = True): + if analysis_identifier is None: + analysis_identifier = f"fit_{fit_module}.{fit_class}" + + # TODO: store a FitDescription object instead of all of these attributes... self.fit_class_name = fit_class self.fit_module = fit_module self.data = data @@ -359,6 +370,7 @@ def __init__(self, self.analysis_identifier = analysis_identifier self.initial_values = initial_values or {} self.exported_results = exported_results or [] + self.export_chi2 = export_chi2 self._result_channels = [] for result in self.exported_results: @@ -372,6 +384,14 @@ def __init__(self, display_hints={"error_bar_for": channel.path}) self._result_channels.append(err_channel) + if self.export_chi2 is not None and 'y_err' in data.keys(): + self._result_channels.append( + OpaqueChannel(path=f"{analysis_identifier}_chi2", + description="Fit Chi-Squared value")) + self._result_channels.append( + OpaqueChannel(path=f"{analysis_identifier}_p", + description="Fit p-value")) + duplicate_channels = [ channel.path for channel, count in collections.Counter(self._result_channels).items() @@ -413,17 +433,9 @@ def describe_online_analyses( if isinstance(v, ResultChannel) ] - analysis_identifier = self.analysis_identifier - if analysis_identifier is None: - # By default, mangle fit type and channels into a pseudo-unique identifier, - # which should work for the vast majority of cases (i.e. unless the user - # creates needlessly duplicate analyses). - analysis_identifier = ("fit_" + f"{self.fit_module}.{self.fit_class_name}" + - "_" + "_".join(channels)) - def analysis_ref(key): return AnnotationValueRef("online_result", - analysis_name=analysis_identifier, + analysis_name=self.analysis_identifier, result_key=key) fit_params = self.fit_klass.get_params() @@ -452,18 +464,17 @@ def analysis_ref(key): parameters={"associated_channels": channels})) analysis = { - analysis_identifier: - FitDescription( - fit_class_name=self.fit_class_name, - fit_module=self.fit_module, - data={ - name: context.describe_coordinate(obj) - for name, obj in self.data.items() - }, - param_bounds=self.bounds, - fixed_params=self.constants, - initial_values=self.initial_values, - ) + self.analysis_identifier: + FitDescription(fit_class_name=self.fit_class_name, + fit_module=self.fit_module, + data={ + name: context.describe_coordinate(obj) + for name, obj in self.data.items() + }, + param_bounds=self.bounds, + fixed_params=self.constants, + initial_values=self.initial_values, + export_chi2=self.export_chi2) } return [a.describe(context) for a in annotations], analysis @@ -485,15 +496,21 @@ def execute( # TODO: Generalise to higher-dimensional fits. x = axis_data[self.data['x']._store.identity] y = result_data[self.data['y']] + y_err = result_data[self.data.get('y_err')] fit = self.fit_klass( - x=np.asarray(x), - y=np.asarray(y), - y_err=None, + x=x, + y=y, + y_err=y_err, param_bounds=self.bounds, fixed_params=self.constants, initial_values=self.initial_values, ) + if self.export_chi2 and y_err is not None: + chi2, fit_p = fit.chi2 + self._result_channels[f"{self.analysis_identifier}_chi2"].push(chi2) + self._result_channels[f"{self.analysis_identifier}_p"].push(fit_p) + for result in self.exported_results: p, p_err = getattr(fit, result.fit_parameter) channel = self._result_channels[result.result_name] diff --git a/ndscan/fitting/fitting.py b/ndscan/fitting/fitting.py index 499be923..5bc34e1f 100644 --- a/ndscan/fitting/fitting.py +++ b/ndscan/fitting/fitting.py @@ -2,7 +2,7 @@ import functools import collections from typing import Dict, List, Optional, Tuple, Union -from scipy import optimize +from scipy import optimize, stats __all__ = ['FitBase'] @@ -39,9 +39,9 @@ def __init__( :param initial_values: dictionary specifying initial parameter values to use in the fit. These override the values found by :meth estimate_parameters: """ - self._x = x - self._y = y - self._y_err = y_err + self._x = None if x is None else np.asarray(x) + self._y = None if x is None else np.asarray(y) + self._y_err = None if y_err is None else np.asarray(y_err) self._validate_param_names(fixed_params) if fixed_params is not None: @@ -171,6 +171,31 @@ def residuals(self): """Returns the fit residuals.""" return self._y - self._func_free(self._x, self._p) + @property + def chi2(self) -> Tuple[float, float]: + """ Chi-Squared test. + + Assumes the values of `y` are independent and normally distributed. + Returns the tuple (Chi-Squared, p-value) + + The the p-value is the probability that the observed deviations from the fit + could be observed by chance assuming the data is normally distributed (values + close to unity indicate high-quality fits). + """ + if self._y_err is None: + raise ValueError("Cannot calculate chi square without knowing y_err") + + n = len(self._x) - len(self._free_params) + if n < 1: + raise ValueError("Cannot calculate chi squared with " + f"{len(self._free_params)} fit parameters and only " + f"{len(self._x)} data points.") + + y_fit = self.evaluate(self._x)[1] + chi_2 = np.sum(np.power((self._y - y_fit) / self._y_err, 2)) + p = stats.chi2.sf(chi_2, n) + return chi_2, p + def _calculate_derived_params(self): """ Updates fit results with values and uncertainties for derived parameters. From f2266d548e97c002167d5a5077eef71c29161877 Mon Sep 17 00:00:00 2001 From: hartytp Date: Fri, 12 Aug 2022 10:08:57 +0100 Subject: [PATCH 27/27] fix tests --- ndscan/experiment/default_analysis.py | 5 +++-- test/test_experiment_entrypoint.py | 15 ++++++++------- test/test_experiment_subscan.py | 15 ++++++++------- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/ndscan/experiment/default_analysis.py b/ndscan/experiment/default_analysis.py index 93f1eae4..5394a2ff 100644 --- a/ndscan/experiment/default_analysis.py +++ b/ndscan/experiment/default_analysis.py @@ -360,7 +360,8 @@ def __init__(self, export_chi2: bool = True): if analysis_identifier is None: - analysis_identifier = f"fit_{fit_module}.{fit_class}" + channels = [v.path for v in data.values() if isinstance(v, ResultChannel)] + analysis_identifier = f"fit_{fit_module}.{fit_class}_" + "_".join(channels) # TODO: store a FitDescription object instead of all of these attributes... self.fit_class_name = fit_class @@ -496,7 +497,7 @@ def execute( # TODO: Generalise to higher-dimensional fits. x = axis_data[self.data['x']._store.identity] y = result_data[self.data['y']] - y_err = result_data[self.data.get('y_err')] + y_err = result_data.get(self.data.get('y_err')) fit = self.fit_klass( x=x, diff --git a/test/test_experiment_entrypoint.py b/test/test_experiment_entrypoint.py index ffd8dc4b..f50131c9 100644 --- a/test/test_experiment_entrypoint.py +++ b/test/test_experiment_entrypoint.py @@ -153,22 +153,22 @@ def test_run_1d_scan(self): "coordinates": {}, "data": { "a": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "a" }, "fwhm": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "fwhm" }, "x0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "x0" }, "y0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "y0" } @@ -181,14 +181,14 @@ def test_run_1d_scan(self): }, "coordinates": { "axis_0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "x0" } }, "data": { "axis_0_error": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_result", "kind": "online_result", "result_key": "x0_error" } @@ -198,11 +198,12 @@ def test_run_1d_scan(self): [curve_annotation, location_annotation]) self.assertEqual( json.loads(self.dataset_db.get("ndscan.rid_0.online_analyses")), { - f"{pref}lorentzian_channel_result": { + f"{pref}lorentzian_result": { "kind": "fit_description", "fixed_params": { "y0": 1.0 }, + 'export_chi2': True, "data": { "y": "channel_result", "x": "axis_0" diff --git a/test/test_experiment_subscan.py b/test/test_experiment_subscan.py index 389eb6b3..eeede8d7 100644 --- a/test/test_experiment_subscan.py +++ b/test/test_experiment_subscan.py @@ -63,22 +63,22 @@ def test_1d_result_channels(self): "coordinates": {}, "data": { "a": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "a" }, "fwhm": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "fwhm" }, "x0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "x0" }, "y0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "y0" } @@ -91,14 +91,14 @@ def test_1d_result_channels(self): }, "coordinates": { "axis_0": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "x0" } }, "data": { "axis_0_error": { - "analysis_name": f"{pref}lorentzian_channel_result", + "analysis_name": f"{pref}lorentzian_child/result", "kind": "online_result", "result_key": "x0_error" } @@ -107,7 +107,7 @@ def test_1d_result_channels(self): self.assertEqual(spec["annotations"], [curve_annotation, location_annotation]) self.assertEqual( spec["online_analyses"], { - f"{pref}lorentzian_channel_result": { + f"{pref}lorentzian_child/result": { "kind": "fit_description", "fixed_params": { "y0": 1.0 @@ -118,6 +118,7 @@ def test_1d_result_channels(self): "x0": [-np.inf, np.inf], "y0": [-np.inf, np.inf], }, + 'export_chi2': True, "data": { "y": "channel_result", "x": "axis_0"