ADR suggestion SampleModel

[henrikjacobsenfys]

General

The SampleModel is briefly described in #2 , but needs to be fleshed out more. I first repeat what was written in #2

This holds the model, which consists of components, which can be Gaussian, Lorentzian, DHO, etc., or user defined. (Care must be taken to also allow 2d models often used in QENS, such as JumpDiffusion. This will probably be handled by a separate class, which can then generate a SampleModel for each probed Q).
It has methods to create and remove components.

I will discuss the convolution and some issues relating to it in another ADR.

Current Implementation

Only a basic skeleton has been implemented.

Proposed implementation

SampleModel

The SampleModel will contain components as discussed below. I think it will inherit from BaseObj. It will at least have the following methods:
add_component: adds a ModelComponent. It optionally takes a new name of the component, otherwise the name of the component will be used. It should check to avoid duplicate names.
remove_component: removes a component using its name
evaluate: evaluates the full model at given x, including detailed balancing. Should give a warning if the temperature varies between components
set_temperature: sets temperature of all components (for detailed balancing). temperature is a Parameter, which is shared between components.
fit_temperature True or False, defaults to False as it`s rare (but not unheard of) to fit the temperature.

It is almost always necessary to allow for a small offset in the x axis. This offset should be identical for all components. It should have a setter, getter, and a method to fix/free it>
offset
fix_offset: True or False. Defaults to False

ModelComponent

Each component will inherit from ModelComponent, which inherits from BaseObj. It will at least have the following methods:
evaluate: Evaluates the component at given x
set_temperature: sets temperature of this component (useful if people want to play around with the components without making a SampleModel). temperature is a Parameter, which is shared between components. It will by default be fixed.
use_detailed_balancing True or False. Default True if temperature is non-zero.
fix_all: fixes all Parameters of this component, except temperature
free_all: frees all Parametersof this component, excepttemperature`

args:
name a name to refer to the component. There should be checks to disallow duplicate names in a SampleModel

GaussianComponent

Here is one example of a ModelComponent. It will take the following arguments and have the following methods:
args:
name
amplitude: amplitude of the Gaussian, defaults to None
width: sigma of the Gaussian
center: center of the Gaussian
area: area of the Gaussian, defaults to None. Exactly one of amplitude or area must be given. Internally, we will work with amplitude, but area is very useful both internally and externally since its conserved during convolution unit: unit of the width, centerandarea`, as a string.

The args are given as doubles, and Parameters will be created from them.

methods:
evaluate evaluates the model
area (both setter and getter):
unit (setter and getter): change/get the unit
fix_amplitude
fix_width
fix_center

LorentzianComponent

Will be nearly identical to the GaussianComponent

DeltaFunctionComponent

Will be nearly identical to the GaussianComponent, except its evaluate will be amplitude if x is 0, and 0 otherwise

[damskii9992]

Since the temperature attribute of a SampleModel is a Parameter, you do not need specific methods to set or make it fittable. Just use its own fixed attribute.

[damskii9992]

I would expect a free_all method to also free the temperature parameter if the temperature parameter has been set on the component. If it is set on the overarching SampleModel class, then it makes sense that it is unchanged, as it then doesn't belong to the component. And same for fix_all.

[henrikjacobsenfys]

This is a difficult choice, and either option would make some users unhappy. I'm inclined to only fit temperature if people explicitly ask for it, since it shouls be a fixed parameter for 95% of experiments

[damskii9992]

Regarding the name of each component, I would suggest simply using the unique_names as those are already checked for duplicates and can be customized by the user (but not re-used). Since you will inherit from object_base, the unique_name parameter will be there already.

[henrikjacobsenfys]

Thanks, I'll do this

[damskii9992]

Regarding the GaussianComponent, the unit is already an attribute on each parameter. I would simply add a check on the initialization that the unit of center matches the unit of width.
You also cannot have unit setters as this is something we have explicitly not allowed due to models potentially breaking. This is our design choice. Users can use the convert_unit method to convert between ex. m and cm, but if they want another dimension they will have to make a new model object.

[henrikjacobsenfys]

It's my understanding that we don't want to encourage users to use Parameters directly, hence a "global" unit for the component.

[damskii9992]

It's true we don't want users to create their own Parameters, this should be done for them when the classes are created, but they can still interact with the created Parameters.
But yes, since they don't create their own Parameters, the unit should be arguments to the constructors, my bad. But don't you also need a unit for the amplitude?
And you should still not be allowed to change the unit once the model has been created, only convert it between units of the same dimension.

[henrikjacobsenfys]

The unit for Amplitude is a.u. (arbitrary units), so it's not needed.

I'll think some more about setting units. I agree with you, but I can also see users be annoyed that they can't specify the units later. Perhaps I'll just default to meV if no unit is given. I suppose It'll evolve with the Measurement class, since the main point is to make sure the experiment and model units match/are converted properly

[henrikjacobsenfys]

I suppose some instruments might allow easy normalisation to absolute units, but I'll ignore that for now

[henrikjacobsenfys]

There will also be a DoubleGaussian and DoubleLorentzian, with pairs of peaks equidistant from zero, with identical width and detailed balancing

[rozyczko]

The unit for Amplitude is a.u. (arbitrary units), so it's not needed.

I'd propose arb. rather than a.u. which can get confused with commonly used atomic unit

[henrikjacobsenfys]

Yes, it'll probably be labelled arb. unit or similar. But in terms of unit handling, its unit is None or Dimensionless.

a.u. is more or less the standard, and, amusingly, sometimes also used for absolute units.

[damskii9992]

This bothers and disturbs me to no end . . .
Even if the community has a standard, can we please choose to ignore it in favour of not being idiots? I.e. a.u. is only atomic units throughout our codebase. Or better yet, don't ever use a.u. for anything, always write out atomic units, arbitrary units or absolute units?
Btw, what is absolute units?

[henrikjacobsenfys]

Don't worry, we'll do arb. units or something similar that is not ambiguous (I know that arbitrary units and "not ambiguous" a bit of an oxymoron, but it is the accepted standard that everyone knows and understands). Arbitrary units can be a bit long, hence the arb.
Arbitrary units are arbitrary, because they're obtained by taking the counts and dividing by the monitor count. The monitor detects a small, generally unknown, fraction of the neutrons at some point before the sample, so it's a measure of the lux of the incoming beam. But the actual numbers depend on incoming energy, chopper settings, slit settings, detector efficiency, sample dimensions etc. Count/monitor is often something on the order of $10^{-4}$, and it's not uncommon to just scale all data by e.g. $10^4$, so that the numbers are on the order of unity to make the plots nicer. The point is that the scale does not matter, it's arbitrary.

For absolute normalisation, a standard (typically vanadium or a phonon from the sample) needs to be measured and careful calibration must be carried out to account for all of the above scaling. Absolute units for inelastic scattering is usually mbarn/sr/meV/f.u.

barn is the unit for scattering cross section
sr is sterradian, i.e. the result is per solid angle
meV you know; the result is per energy (for elastic scattering/diffraction, absolute units is mbarn/sr/f.u.)
f.u. is formula unit. Some people do per atom instead.

[rozyczko]

ModelComponent:

temperature is a Parameter, which is shared between components.

Then maybe we also need a model collection, which will have temperature as an attribute for all the components in it?
Then we can possibly do calculations on the whole collection, rather than a single component, if this is something you need in your workflow.

[henrikjacobsenfys]

The SampleModel is a kind of model collection, as it consists of multiple components.

[rozyczko]

ModelComponent:

free_all: frees all Parametersof this component, except temperature

Making temperature not adhere to the common way of changing flags is unusual. If you want to underline the fact that temperature is being fitted on very rare occassions, maybe it's better to assign a special flag/attribute to the component, so users have to be specific. E.g.

component.fit_temperature = True

[henrikjacobsenfys]

Yes, this is exactly what I'll do. Temperature should only ever be fitted if the thermometer is broken (has happened to me more than once), or if there's bad thermal contact to the sample (which has also happened to me). So fitting temperature is not needed in most experiments, but will be needed at least a few times by most experimenters.

[rozyczko]

GaussianComponent:

amplitude: amplitude of the Gaussian, defaults to None

I htink it is dangerous to have a badly behaving default.
If I do
my_gaussian = GaussianComponent()
I probably expect some generic gaussian with unity area or amplitude rather than something with only Nones.

[henrikjacobsenfys]

Yes, good point. I'll change this. It's tied to the discussion in #6 about how to describe a peak: via its amplitude $Amp$
$Amp \cdot \exp(-(x/\sigma)^2)$, its area $A$:
$A /\sqrt{2\pi \sigma^2}\exp(-(x/\sigma)^2)$,
or its height $H$:
$H \cdot F/max(F)$, where $F=R*DBF\cdot 1/\sqrt{2\pi \sigma^2}\exp(-(x/\sigma)^2)$, where DBF is the detailed balancing factor described in #6

I think I'll stick to just $A$, as that is the physically meaningful parameter. I may make a helper method to estimate it based on the height of the peak.