Add Data and Function Background types (#276)

* Staging commit * Copied the Standard Layers example to 'miscelaneous'. Renamed it 'backgroundTypes' - this will be the example folder for this * Modified projectClass to allow data backgrounds * Data background implemented. Needs testing * Implemented all background types for standard layers for evaluation (does not compile) * Adds data and custom file names as fields to project "toStruct" struct * Fixes build and tests * Adds code to expand all contrast data to six columns * Applies the "function" background type to all calculation and model types * Moves "applyBackgroundFunction" into the contrast loop as "constructBackground" * Modifies "constructBackground" to work for unequal data and simulation sizes * Modifies "constructBackground" to work from simulation limits * Removes background parameters from "applyBackgroundCorrection", using background array instead * Refactors results struct to replace "backgroundParams" with "backgrounds" * Removes variable definitions unnecessary for compilation * Adds routine "makeSimulationRange" * Fixes array size discrepancies and "shiftedDat" bug * Adds cell array "contrastBackgroundTypes" * Converts "contrastBackgroundParams" to a cell array * Adds "source" column to background and resolution tables * Renames "addDataBackgroundToContrastData" as "insertDataBackgroundIntoContrastData" * Adds code to raise error if no contrasts are defined * Disables support for non-matlab background functions * Moves reduction of "shiftedData" to three columns to "applyBackgroundCorrection" * Reinstates checkIndices for backgroundParams * Tidies up code * Adds test for insert data routine and includes new examples in tests * Tidies up code and removes unecessary routines * Addresses review comments * Addresses further review comments --------- Co-authored-by: arwelHughes <[email protected]> Co-authored-by: arwel <[email protected]>
RascalSoftware · Nov 27, 2024 · 35f0252 · 35f0252
1 parent c944de7
commit 35f0252
Show file tree

Hide file tree

Showing 99 changed files with 1,981 additions and 1,465 deletions.
diff --git a/API/RATMain.m b/API/RATMain.m
@@ -10,40 +10,45 @@
 result = makeEmptyResultStruct(problemStruct.numberOfContrasts, length(problemStruct.fitParams), domains);
 bayesResults = makeEmptyBayesResultsStruct(problemStruct.numberOfContrasts, domains, controls.nChains);
 
-% Decide what we are doing....
-switch controls.procedure
-    case coderEnums.procedures.Calculate % Just a single reflectivity calculation
+if problemStruct.numberOfContrasts > 0
+    switch controls.procedure
+        case coderEnums.procedures.Calculate % Just a single reflectivity calculation
+            controls.calcSldDuringFit = true;
+            result = reflectivityCalculation(problemStruct,problemCells,problemLimits,controls);
+        case coderEnums.procedures.Simplex
+            if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
+                triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning simplex\n\n'));
+            end
+            [problemStruct,result] = runSimplex(problemStruct,problemCells,problemLimits,controls);
+        case coderEnums.procedures.DE
+            if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
+                triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning Differential Evolution\n\n'));
+            end
+            [problemStruct,result] = runDE(problemStruct,problemCells,problemLimits,controls);
+        case coderEnums.procedures.NS
+            if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
+                triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning Nested Sampler\n\n'));
+            end            
+            [problemStruct,result,bayesResults] = runNestedSampler(problemStruct,problemCells,problemLimits,controls,priors);   
+        case coderEnums.procedures.Dream
+            if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
+                triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning DREAM\n\n'));
+            end
+            [problemStruct,result,bayesResults] = runDREAM(problemStruct,problemCells,problemLimits,controls,priors);
+        otherwise
+            error('The procedure "%s" is not supported. The procedure must be one of "%s"', controls.procedure, strjoin(fieldnames(coderEnums.procedures), '", "'));
+    end
+
+    % Then just do a final calculation to fill in SLD if necessary
+    % (i.e. if calcSLD is no for fit)
+    if ~controls.calcSldDuringFit
+        controls.calcSldDuringFit = true;
+        controls.procedure = coderEnums.procedures.Calculate;
         result = reflectivityCalculation(problemStruct,problemCells,problemLimits,controls);
-    case coderEnums.procedures.Simplex
-        if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
-            triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning simplex\n\n'));
-        end
-        [problemStruct,result] = runSimplex(problemStruct,problemCells,problemLimits,controls);
-    case coderEnums.procedures.DE
-        if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
-            triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning Differential Evolution\n\n'));
-        end
-        [problemStruct,result] = runDE(problemStruct,problemCells,problemLimits,controls);
-    case coderEnums.procedures.NS
-        if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
-            triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning Nested Sampler\n\n'));
-        end            
-        [problemStruct,result,bayesResults] = runNestedSampler(problemStruct,problemCells,problemLimits,controls,priors);   
-    case coderEnums.procedures.Dream
-        if ~strcmpi(controls.display, coderEnums.displayOptions.Off)
-            triggerEvent(coderEnums.eventTypes.Message, sprintf('\nRunning DREAM\n\n'));
-        end
-        [problemStruct,result,bayesResults] = runDREAM(problemStruct,problemCells,problemLimits,controls,priors);
-    otherwise
-        error('The procedure "%s" is not supported. The procedure must be one of "%s"', controls.procedure, strjoin(fieldnames(coderEnums.procedures), '", "'));
-end
+    end
 
-% Then just do a final calculation to fill in SLD if necessary
-% (i.e. if calcSLD is no for fit)
-if ~controls.calcSldDuringFit
-    controls.calcSldDuringFit = true;
-    controls.procedure = coderEnums.procedures.Calculate;
-    result = reflectivityCalculation(problemStruct,problemCells,problemLimits,controls);
+else
+    error("RAT cannot proceed without at least one contrast defined in the project")
 end
 
 end
diff --git a/API/checkIndices.m b/API/checkIndices.m
@@ -1,13 +1,30 @@
-function checkIndices(problemStruct)
+function checkIndices(problemStruct, customFiles)
     % Make sure that the indices provided lie within the bounds of the
     % corresponding array.
 
     numBackgroundParams = length(problemStruct.backgroundParams);
+    numCustomFiles = length(customFiles);
     for i = 1:length(problemStruct.contrastBackgroundParams)
-        index = problemStruct.contrastBackgroundParams(i);
-        if (index < 1 && index ~= -1) || index > numBackgroundParams
-            throw(exceptions.indexOutOfRange(sprintf('contrastBackgroundParams(%i) is %i, which is outside the range of backgroundParams', i, index)));
+
+        indices = problemStruct.contrastBackgroundParams{i};
+
+        if length(indices) > 1
+            % The first index is a custom file, the rest are background parameters
+            index = indices(1);
+            if index < 1 || index > numCustomFiles
+                throw(exceptions.indexOutOfRange(sprintf('contrastBackgroundParams{%i}(1) is %i, which is outside the range of customFiles', i, index)));
+            end
+            indices = indices(2:end);
         end
+
+        % All of these indices are background parameters
+        for j = 1:length(indices)
+            index = indices(j);
+            if index < 1 || index > numBackgroundParams
+                throw(exceptions.indexOutOfRange(sprintf('contrastBackgroundParams{%i}(%i) is %i, which is outside the range of backgroundParams', i, j, index)));
+            end
+        end
+
     end
 
     numQzshifts = length(problemStruct.qzshifts);

diff --git a/API/insertDataBackgroundIntoContrastData.m b/API/insertDataBackgroundIntoContrastData.m
@@ -0,0 +1,25 @@
+function contrastData = insertDataBackgroundIntoContrastData(contrastData,backgroundData)
+% Deal with a Data background. The data and errors in this case are
+% inserted into columns 5 and 6 of the relevant datafile.
+%
+% Currently, we throw an error if the qz column (column 1) of the two
+% datafiles are different. Eventually this will be replaced by an
+% interpolation to make it more general.
+
+% Get the arrays from the cells
+dataArray = contrastData{:};
+backgroundArray = backgroundData{:};
+
+% Check that we have the same q
+if ~isequal(dataArray(:,1), backgroundArray(:,1))
+    throw(exceptions.invalidValue("q points must be equal for Data and Background Data"));
+end
+
+% Insert background data into columns 5 and 6 of contrast data
+dataArray(:,5) = backgroundArray(:,2);
+dataArray(:,6) = backgroundArray(:,3);
+
+% Package as a cell array for output...
+contrastData = {dataArray};
+
+end
diff --git a/API/makeEmptyResultStruct.m b/API/makeEmptyResultStruct.m
@@ -2,8 +2,9 @@
     % A function to make an empty container to hold the results of
     % reflectivity calculations. The struct has the following format:
     %
-    % nPar = number of fitted parameters
+    % nParams = number of fitted parameters
     % nContrasts = number of contrasts
+    % nDomains = number of domains - 1 for non-polarised, 2 for domains
     %
     % result = 
     % 
@@ -12,9 +13,10 @@
     %       reflectivity: [nContrastsx1 cell]
     %         simulation: [nContrastsx1 cell]
     %        shiftedData: [nContrastsx1 cell]
-    %          layerSlds: [nContrastsx1 cell]
-    %        sldProfiles: [nContrastsx1 cell]
-    %    resampledLayers: [nContrastsx1 cell]
+    %        backgrounds: [nContrastsx1 cell]
+    %          layerSlds: [nContrastsxnDomains cell]
+    %        sldProfiles: [nContrastsxnDomains cell]
+    %    resampledLayers: [nContrastsxnDomains cell]
     % calculationResults: [1x1 struct]
     %     contrastParams: [1x1 struct]
     %          fitParams: [1xnParams double]
@@ -36,8 +38,6 @@
     % --------------------------------------------------------------------
     % (2) result.contrastParams
 
-    backgroundParams = zeros(nContrasts,1);
-    coder.varsize('backgroundParams',[maxArraySize 1],[1 0]);
     scalefactors = zeros(nContrasts,1);
     coder.varsize('scalefactors',[maxArraySize 1],[1 0]);
     bulkIn = zeros(nContrasts,1);
@@ -51,8 +51,7 @@
     resample = zeros(1, nContrasts);
     coder.varsize('resample',[1 maxArraySize],[0 1]);
 
-    contrastParams = struct('backgroundParams', backgroundParams, ...
-                            'scalefactors', scalefactors, ...
+    contrastParams = struct('scalefactors', scalefactors, ...
                             'bulkIn', bulkIn, ...
                             'bulkOut', bulkOut, ...
                             'resolutionParams', resolutionParams, ...
@@ -64,27 +63,33 @@
 
     reflectivity = cell(nContrasts,1);
     refCell = ones(2,2);
-    coder.varsize('refCell',[10000 2],[1 0]);
     for i = 1:nContrasts
         reflectivity{i} = refCell;
     end
+    coder.varsize('reflectivity{:}',[10000 2],[1 0]);
 
     simulation = cell(nContrasts,1);
     simCell = ones(2,2);
-    coder.varsize('simCell',[10000 2],[1 0]);
     for i = 1:nContrasts
         simulation{i} = simCell;
     end
+    coder.varsize('simulation{:}',[10000 2],[1 0]);
 
     shiftedData = cell(nContrasts,1);
     shiftCell = ones(2,3);
-    coder.varsize('shiftCell',[10000 3],[1 0]);
     for i = 1:nContrasts
         shiftedData{i} = shiftCell;
     end
+    coder.varsize('shiftedData{:}',[10000 3],[1 0]);
+
+    backgrounds = cell(nContrasts,1);
+    backgroundCell = ones(2,3);
+    for i = 1:nContrasts
+        backgrounds{i} = backgroundCell;
+    end
+    coder.varsize('backgrounds{:}',[10000 3],[1 0]);
 
     layerSldCell = ones(2,3);
-    coder.varsize('layerSldCell',[10000 6],[1 1]);    
     if domains
         layerSlds = cell(nContrasts,2);
         for i = 1:nContrasts
@@ -97,9 +102,9 @@
             layerSlds{i} = layerSldCell;
         end
     end
+    coder.varsize('layerSlds{:}',[10000 6],[1 1]);
 
     sldProfileCell = ones(2,2);
-    coder.varsize('sldProfileCell',[10000 2],[1 0]);
     if domains
         sldProfiles = cell(nContrasts,2);
         for i = 1:nContrasts
@@ -108,14 +113,13 @@
         end
     else
         sldProfiles = cell(nContrasts,1);
-
         for i = 1:nContrasts
             sldProfiles{i} = sldProfileCell;
         end
-    end   
+    end
+    coder.varsize('sldProfiles{:}',[10000 3],[1 1]);
 
     resampledLayersCell = ones(2,3);
-    coder.varsize('resampledLayersCell',[10000 3],[1 0]);
     if domains
         resampledLayers = cell(nContrasts,2);
         for i = 1:nContrasts
@@ -128,20 +132,21 @@
             resampledLayers{i} = resampledLayersCell;
         end
     end
+    coder.varsize('resampledLayers{:}',[10000 4],[1 1]);
 
     fitParams = zeros(nParams,1);
     coder.varsize('fitParams',[maxArraySize 1],[1 0]);
 
     fitNames = cell(nParams,1);
-    fitNamesChar = '';
-    coder.varsize('fitNamesChar',[1 maxArraySize],[0 1]);
     for i = 1:nParams
-        fitNames{i} = fitNamesChar;
+        fitNames{i} = '';
     end
+    coder.varsize('fitNames{:}',[1 maxArraySize],[0 1]);
 
     result = struct('reflectivity', {reflectivity}, ...
                     'simulation', {simulation}, ...
                     'shiftedData', {shiftedData}, ...
+                    'backgrounds', {backgrounds}, ...
                     'layerSlds', {layerSlds}, ...
                     'sldProfiles', {sldProfiles}, ...
                     'resampledLayers', {resampledLayers}, ...