ADR suggestion Example of full workflow in lib

[henrikjacobsenfys]

General

I have attempted to write down a complete data analysis workflow in the lib, to figure out what methods are required and what the structure should be. Overall, a complete workflow consists of the following steps:

1. Load resolution data and fit it
2. Load data, inspect data, attach a resolution function to it, and fit it
3. Inspect and fit the resulting fit parameters
4. (Optional) Fit data to a global model (This can also be instead of 2)

Current implementation

Nothing is implemented

Suggested implementation

1. Resolution

import EasySpectroscopy as ESpec
from EasyScience.Objects.variable import Parameter

# Initialize project and add data
Vanadium=ESpec.Project()
Vanadium.add_data('VanadiumData') # Could be a folder with data or a file.. depends on the format.

# Generate a model to fit the data using either built in models (Gaussian, Lorentzian, DHO, a few others?) or user defined models. 
ResModel=ESpec.SampleModel()
ResModel.add_component('Gaussian',
                       amplitude=Parameter('G1amp',value=0.1,unit=None,fixed=False),
                       sigma=Parameter('G1sigma',value=0.1,unit='meV',fixed=False)) # if no center is given, it will be zero. 
 
ResModel.add_component('Lorentzian',
                       amplitude=Parameter('L1amp',value=0.1,unit=None,fixed=False),
                       gamma=Parameter('L1gam',value=0.2,unit='meV',fixed=False))

ResModel.EnergyOffset=Parameter('EOffset',value=0.0,unit='meV',fixed=False) #It defaults to 0 and fixed.

# Add a model of the background
BGModel=ESpec.SampleModel()
BGModel.add_component('Constant',amplitude=Parameter('Const',value=0.1,unit=None,fixed=False))

# Attach models to the project
Vanadium.SampleModel=ResModel #This copies the Res model in some way to match the size of the data
Vanadium.BackgroundModel=BGModel
Vanadium.ResolutionModel=None #This is the default

# Inspect the data and model before fitting
Vanadium.plot_data()
Vanadium.plot_model() 

# Fit data and plot the fit
Vanadium.fit_data(method='sequential')  #Will fit all the data, using the output of one fit as input for the next. Also needs to allow to fit only a select data set using some kind of indexing

Vanadium.plot_fit() #will plot the fit on top of the data

2. Initial fit of data

This will be almost identical to the fit of vanadium, so I will not quote the whole thing here. The important change is the resolution modeL

MyProject=ESpec.Project()
MyProject.add_data('MyData')

SampleModel=ESpec.SampleModel()
SampleModel.add_component('Delta',
                       amplitude=Parameter('D1amp',value=0.1,unit=None,fixed=False)) 
SampleModel.add_component('Lorentzian',
                       amplitude=Parameter('L1amp',value=0.1,unit=None,fixed=False),
                       gamma=Parameter('L1gam',value=0.2,unit='meV',fixed=False))
....

MyProject.ResolutionModel=Vanadium.Model #Should also be possible to do Vanadium.Data or something similar, in case raw data is desired
MyProject.ResolutionModel.fix_all() #We don't want to fit the resolution again
...

Inspect and fit the resulting fit parameters

# Plot fit parameters
MyProject.FitParameter['L1gam'].plot() #will plot the fit parameter as a function of Q. The x axis is the coordinate of the data set. The .plot method can take input to customise the plot

SimpleDiffusionModel=ParameterModel()
SimpleDiffusionModel.add_component('SimpleDiffusion',A=Parameter(value=1,unit=None,fixed=False)) #This is a simple diffusion model, with a single parameter A, giving the width of the Lorentzian as function of Q
                             
MyProject.FitParameter['L1gam'].Model=SimpleDiffusionModel # Attach the parameter model to the fit parameter.

MyProject.FitParameter['L1gam'].fit_data()
MyProject.FitParameter['L1gam'].plot_fit() #will plot the fit on top of the data

Global model

In a global model, there will still be local parameters. I therefore think it makes sense to link parameters via a ParameterModel:

MyProject.link_fit_parameter('L1gam', model=SimpleDiffusionModel)
MyProject.fit_data()

Then all the data will be fitted simultaneously. If we e.g. have data at 10 values of Q, we will fit 10 amplitudes, 10 backgrounds, but only one parameter to describe gamma and how it varies with Q.

Some thoughts

• Each cut of the data should have a model attached. This should happen under the hood when the SampleModel is assigned. However, the data will most likely be a Scipp array, and so some binning of the data is required. If the binning changes, what happens to the models? The assignment probably needs to be explicit using the Analysis class, but I've reached the end of my concentration.
• A data set will in general be intensity as function of Q and E for various temperatures. Need to think carefully about binning.
• Detailed Balancing is very important to get right. Can the analytical convolution even be used in this case?
• Users will probably want to compare fits of different samples

Overall structure and overview of methods

Project

The Project is intended to own a complete experiment on one sample. It contains Data (added by add_data), SampleModel, BackgroundModel, ResolutionModel, FitParameter. It has methods to plot (plot_data) and fit (fit_data) the data. It will use the ResolutionHandler described in #2. It will also have a method to link fit parameters using global models (e.g. to diffusion) models (link_fit_parameter).

A SampleModel added to the Project will attach a SampleModel to each data set in Data.
A BackgroundModel is also a SampleModel, but will not be convoluted with the resolution function
A ResolutionModel will for now be restricted to also being a SampleModel, but can be expanded to be a data set to use raw data for the resolution convolution.

Data

Contains the data with various methods to be defined. Each data set can be given a SampleModel

SampleModel

Also discussed in #2 . Contains a model of the sample, consisting of components. Has methods to add and remove components, as well as to evaluate them on an input x axis.

FitParameter

Fitted parameters are a kind of data set on their own. Perhaps they can even use the same class as Data, but that's unclear for now. It needs methods to select and to plot specific parameters.

[rozyczko]

This looks solid and well structured.
One point I would like to raise is the explicit use of Parameters. We want to hide the Parameter itself away and only create attributes which will internally be expressed as parameters.

E.g.
Instead of

ResModel.add_component('Gaussian',
                       amplitude=Parameter('G1amp',value=0.1,unit=None,fixed=False),
                       sigma=Parameter('G1sigma',value=0.1,unit='meV',fixed=False))

why not use

# add predefined gaussian component, which resides in say, ESpec.models
# and accepts `sigma` and `amplitude` as optional arguments
ResModel.components.add(component_id="g1")
ResModel.components["g1"].type = ESpec.models.gaussian
ResModel.components["g1"].amplitude.value = 0.1
ResModel.components["g1"].sigma.unit = 'meV'
ResModel.components["g1"].sigma.fixed = False

or, create a local component and then assign it with add

component = ESpec.models.gaussian
component.amplitude.value = 0.1
component.sigma.value = 0.1
component.sigma.unit = 'meV'
component.sigma.fixed = False
ResModel.components.add(component)