Global models

[stestoll]

Today, I tried to globally fit two (generated) signals that stem from the same sample, differing only in the modulation depth. Distance distribution and spin concentration are identical. Here is the script.

clc

% Generate distance distribution
n = 201;
r = linspace(2,8,n);
P = dd_gauss(r,[5 0.2]);

% Generate two signals with different modulation depth
conc = 200;
sig = 0.01;
lam1 = 0.25;
lam2 = 0.5;
t = linspace(0,3,n);
B1 = bg_hom3d(t,conc,lam1);
B2 = bg_hom3d(t,conc,lam2);
V1 = dipolarsignal(t,r,P,lam1,B1,'NoiseLevel',sig);
V2 = dipolarsignal(t,r,P,lam2,B2,'NoiseLevel',sig);

% Fit signals individually and globally
dd = @dd_gauss;
bg = @bg_hom3d;
fitsignal(V1,t,r,dd,bg);
fitsignal(V2,t,r,dd,bg);
fitsignal({V1,V2},{t,t},r,dd,bg);

In the simultaneous fit of the two signals, the spin concentration occurs as a parameter twice, once for each signal. It should occur only once.

Clearly, our current setup for global fitting with fitsignal etc. and the background models are not able to take this into account.

This brings up a bigger issue: DeerLab cannot handle arbitrary global models. Currently, we have hard-coded some basis global models into fitsignal, fitparamodel, etc., but it appears to me that there must be a more general and elegant approach. The current way has two issues: (a) it creates complicated code, (b) it limits the models that can be run.

Is there a better way to provide infrastructure for global models?

[luisfabib]

Let's start from what defines a global model:

1) They are models that are fitted to multiple datasets to yield multiple or single distributions, background or dipolar signals. A global model is defined from multiple local models. 2) They use a set of parameters some of which are local and some are global. Some specific parameters from these local models are shared between them.

My proposal is to have a function globalmodel which takes multiple local models (from the built-in selection available in DeerLab) and allows the user to specify which of the parameters in said models are to be global.

The function writes a new model (in the same style as the other built-in ones) which now takes a mixed set of local and global parameters, which are then internally distributed to the corresponding local models. These models then return a cell array with all the fitted distributions, backgrounds or signals.

Here is a working example of this concept (uses the tables format for the info of the parametric models introduced in abaf21e

The two examples show how a global model for a background fit of two signals could be achieved (solving the example in the OP of this issue) and how a parametric two-Gauss titration model could be defined.

The specifics on whether and how this will be implemented will follow. However, since the interface of these models are almost unchanged with respect to the "standard" ones, the implementation would only require some minor adjustments and would allow their use with functions like fitsignal.

%% Example 1: global background model

clf,clc,clear

t1 = linspace(0,6,300);
t2 = linspace(0,3,300);

bg_global = globalmodel({@bg_hom3d,@bg_hom3d},'Spin concentration');

bg_global()

par0 = [150];
B = bg_global({t1,t2},par0);

B1 = B{1};
B2 = B{2};

plot(t1,B1,t2,B2+0.5)
axis tight, grid on
xlabel('t [us]')
ylabel('B(t)')

%% Example 2: global titration model

clf,clc,clear

r = linspace(0,6,300);

globalpar{1} = 'Center of 1st Gaussian';
globalpar{2} = 'Center of 2nd Gaussian';
globalpar{3} = 'FWHM of 1st Gaussian';
globalpar{4} = 'FWHM of 2nd Gaussian';

dd_global = globalmodel(repmat({@dd_gauss2},5,1),globalpar);

dd_global()

par0 = [2 4 0.5 0.5 linspace(0,1,5)];
P = dd_global(r,par0);

P1 = P{1};
P2 = P{2};
P3 = P{3};
P4 = P{4};
P5 = P{5};

plot(r,P1,r,P2+1,r,P3+2,r,P4+3,r,P5+4)
axis tight, grid on
xlabel('r [nm]')
ylabel('P(r)')

function globalModelFcn = globalmodel(models,globalpar)

if ~iscell(models) || numel(models)<2
    error('The 1st input must be a cell array of multiple function handles.')
end

% Detemine number of models to be mixed
nModels = numel(models);

% Containers
idx = 0;
Info = [];
alreadyUsed(1:numel(globalpar)) = false;

% Loop over local models
for i = 1:nModels

    % Get info for current local model
    info = models{i}();
    nparam = height(info);

    % Containers
    duplicates = [];
    paridx_ = zeros(1,nparam);

    % Loop over model parameters
    for j = 1:nparam

        paramName = info.Parameter{j};

        %Check if parameter is global or local
        if any(strcmp(paramName,globalpar))
            n = find(strcmp(paramName,globalpar));
            %If global variable has already been found...
            if alreadyUsed(n)
                %... just continue with the next parameter
                duplicates(end+1) = j;
                paridx_(j) = globidx(n);
                continue;
            end
            %... if not, then store it 
            type = 'Global';
            alreadyUsed(n) = true;
            globidx(n) = j;
        else
            % If local parameter, just format with info of corresponding model 
            type = sprintf('Local %i',i);
        end

        % Store reference to parameter index in local model
        idx = idx+1;
        paridx_(j) = idx;

        % Adapt parameter index
        info.Index(j) = idx;
        % Adapt parameter name
        info.Parameter{j} = sprintf('%s: %s',type,paramName);
    end

    % Store parameter indices for model i
    paridx{i} = paridx_;

    %Delete rows correponding to global variables already used
    info(duplicates,:) = [];

    % Update global model info table 
    Info = [Info; info];
end

% Sort table such that global parameters are first listed
[Info,idxs] = sortrows(Info,2);
Info.Index = (1:1:height(Info)).';

% Mixed model function handle
%-------------------------------------------------------------------------------
globalModelFcn = @globalModel;

% Function to allow request of information structure or model values
    function output = globalModel(varargin)

        % If no parameters are given, return global model info
        if nargin==0
            output = Info;
            return
        end
        % Otherwise just two parameters allowed
        if nargin~=2
            error('Model requires two input arguments.')
        end
        %Parse inputs
        x = varargin{1};
        params = varargin{2};
        % Validation
        if ~iscell(x)
           x = {x}; 
        end
        if numel(x)==1
           x = repmat(x,nModels,1); 
        end
        if numel(x)~= nModels 
            error('The number of axes for this global model cannot exceed %i.',nModels)
        end

        % Invert the sort done above
        params(idxs) = params;

        % Evaluate all the local models
        y = cell(1,nModels);
        for k = 1:nModels
            y{k} = models{k}(x{k},params(paridx{k}));
        end
        output = y;

    end

end