`x` as a matrix

[hartytp]

Currently we the x data for models is nominally an array representing values of a single scan axis (*). In a few places it's been useful to have models which support multiple scan axes (e.g. the Rabi Flop model, which supports varying both time and detuning). Currently this is done by a bit of a hack where we pass a tuple of x vectors.

Since this feels like a broadly useful feature I suggest we make this something that is officially supported. My proposal is to make x into an Array[(num_samples, num_axes), np.float64] everywhere in the package. Then model docstrings should should explicitly state what scan axes they support. Models should also probably have a get_num_axes method (returns 1) by default to allow us to validate data shapes.

NB I do not propose supporting multi-axis fitting (at least for now).

Thoughts @AUTProgram

[hartytp]

(*) I'm struggling a bit with terminology here and don't love the term axis as time and frequency are really model parameters, but that gets confusing with the parameters we're fitting.

I see two options here: the first is to scrap the notion of a special x axis entirely; just make x a normal model parameter like all the others. Then have func expect vectors of model parameters rather than specific values, and tell the fitter which one(s) we want to fit. The second option is to essentially do what we're doing now but have distinct terminology to separate the x parameters from the fitted parameters.

The first of these options is conceptually quite appealing and would make exploring parameter space more ergonomic. However, in practice, I don't think it's actually needed for many use cases and it would make writing models significantly more complex. So, I think we should discount this.

Assuming we stick with the first option we need to settle on terminology. Suggestions welcome, but the best I have come up with so far is axes v parameters.

[AUTProgram]

I think this comes down to the distinction between independent variables and parameters. The function for a Rabi flop, for example, can in the most general case be considered a function of three independent variables $f(t, \delta, \Omega)$ (neglecting readout levels, dead time,...). However, in practice, we often think about it as a function of one independent variable and two parameters, for example $f(t; \delta, \Omega)$, where $\delta$ and $\Omega$ are treated as fixed. So I would actually put it the other way and say that in the most general case $f$ has not got any parameters at all, but only independent variables. As soon as we assume that the data for $f$ has been acquired with one of these "variables" fixed, the "variable" becomes a parameter. So I guess the question is if we want to treat all variables on equal footing or assume that some of them will always be treated as parameters. Will have a think about this.

[hartytp]

Making x a matrix to support multiple x channels (in the same way we treat y) should be relatively straightforward and non-invasive. The downside of that is that the separation between x channels and model parameters is then baked into the model class definition. e.g. one can't take the rabi flop frequency model and float both detuning offset and rabi frequency. However, this might still be a useful capability, e.g. for making a model that fits a 2D gaussian to a camera image to find an ion center position.

One could potentially imagine having some tooling which allows one to automatically generate a subclass of a model with some of its parameters converted into additional x channels. But, even with this, most of the heuristics would break because they are built around the assumption that there is only a single x-axis and all parameters have scalar values.

One could also imagine something more drastic along the lines of your suggestion above @AUTProgram where we replace the notion of an x axis entirely with a set of independent variables. And then allow the division between independent variables and parameters to be set at runtime. That would be more invasive, it would also require support in the heuristics which feels like a lot of work.

[hartytp]

I'm going to close this for now. We can definitely consider a feature along these lines but it needs to be driven by a strong use-case, which we don't have right now.

One of the benefits with the kinds of atomic systems this package is targeted at is that the physics is generally quite separable. This means that data analysis is often best done as a series of hierarchical 1D fits rather than a system where everything couples to everything so n-D fits become important.

[hartytp]

If the main use case here turns out to be doing fits of fits, a nice approach might be to make a MultiFitter wrapper around Fitter that fits along one axis, then takes a fit result and uses that as the y channel for a fix along the next axis.

[hartytp]

Comments:

The general problem of allowing any model parameter to become an independent variable is not something I currently anticipate supporting for two reasons:

1. from the use cases we've seen so far, it's not something that has come up much - most data seems to be well handled by a hierarchical series of 1D scans (fit to one model, take an output of that model and use it as an input to the next model and so on)
2. it is not clear what supporting this would look like. How would we handle parameter estimators? Decide on certain "valid" combinations of parameters to use as independent variables and special-case those?

Whatever we do here, it's going to be an invasive change and the return on effort feels low / the benefit insufficient to justify the additional complexity.

I do, however, see two particular cases which are worth considering:

Hierarchical data

e.g. a 2D scan where we want to fix one model to the first axis and then use fit results as y-channels for a fit to a separate fit model on the second axis. Something like the MultiFitter described above combined with AggregateModel / RepeatedModel would handle this well and without requiring changes to existing code.

Intrinsically multi-dimensional models

The main place this comes up is image analysis. e.g. fitting a 2D gaussian to an image.

The low-effort approach would be to treat the two axes of the image separately (fit to find the x centre, then fit to find the y centre). Then use a MultiFitter as discussed above. That's an okay, if slightly hacky, way of getting a decent answer out which is probably "good enough" for most cases - if not as good as doing the proper multi-d fit.

Given how rarely this comes up, I'm inclined to say that the effort to benefit from implementing this does not currently seem sufficient to justify the effort it would take to do, combined with the amount of additional complexity it would add to the package.

Conclusion

Let's:

• implement the MultiFitter as described above
• add a note about this in the readme
• keep an eye on use-cases and re-open this if there seems to be a genuine need for something else