diff --git a/DataTree/AcquiredData/DataFiles/DataFilesClass.m b/DataTree/AcquiredData/DataFiles/DataFilesClass.m
index 8bd1bf6d..6b38e8d9 100644
--- a/DataTree/AcquiredData/DataFiles/DataFilesClass.m
+++ b/DataTree/AcquiredData/DataFiles/DataFilesClass.m
@@ -16,6 +16,8 @@
     properties (Access = private)
         lookupTable
         excludedFolders
+        savedFilename
+        changed
     end
     
     methods
@@ -30,6 +32,8 @@
             obj.err = -1;
             obj.errmsg = {};
             obj.dirFormats = struct('type',0, 'choices',{{}});
+            obj.lookupTable = [];
+            obj.changed = false;
 
             global logger
             global cfg
@@ -38,7 +42,6 @@
             cfg    = InitConfig(cfg);
             
             obj.logger = logger;
-            obj.lookupTable = [];
             obj.excludedFolders = {};
             
             skipconfigfile = false;
@@ -76,7 +79,7 @@
             obj.rootdir = filesepStandard(obj.rootdir,'full');           
             
             % Configuration parameters
-            obj.config = struct('RegressionTestActive','',  'AskToFixNameConflicts',1, 'DerivedDataDir','');
+            obj.config = struct('RegressionTestActive','',  'AskToFixNameConflicts',1, 'DerivedDataDir','', 'DerivedDataRootDir','');
             if skipconfigfile==false
                 str = cfg.GetValue('Regression Test Active');
                 if strcmp(str,'true')
@@ -110,108 +113,287 @@
                 [obj.rootdir, 'imagerecon'];
                 [obj.rootdir, 'anatomical'];
                 };
+            obj.config.DerivedDataRootDir = [obj.rootdir, obj.config.DerivedDataDir, '/', f, '/'];
+            obj.savedFilename = [obj.config.DerivedDataRootDir, 'DatasetFiles', '_', obj.filetype, '.mat'];
             if nargin==0
                 return;
             end
             
             obj.err = 0;
-            
-            obj.InitFolderFromats();
+            obj.InitFolderFormats();
             obj.GetDataSet();
+            obj.SaveDataSet();
+        end
+        
+        
+        
+        % -----------------------------------------------------------------------------------
+        function Compare(obj)
+            if obj.IsEmpty()
+                return
+            end
+            obj2 = DataFilesClass();
+            obj2.rootdir = obj.rootdir;
+            obj2.filetype = obj.filetype;
+            obj2.dirFormats = obj.dirFormats;
+            obj2.GetDataSet('skipsaveddata')
+            
+            %%%% Create list of all the saved names for sorting later
+            filesSaved = cell(length(obj.files),1);
+            filesSavedErr = cell(length(obj.filesErr),1);
+            for ii = 1:length(obj.files)
+                filesSaved{ii} = obj.files(ii).name;
+            end
+            for ii = 1:length(obj.filesErr)
+                filesSavedErr{ii} = obj.filesErr(ii).name;
+            end
+
+            %%%% Move any files from saved error list that have been updated to nonerror list 
+            N = length(obj.filesErr);
+            for ii = N:-1:1
+                for jj = 1:length(obj2.files)
+                    if strcmp(obj.filesErr(ii).name, obj2.files(jj).name)
+                        if obj.filesErr(ii).datenum ~= obj2.files(jj).datenum
+                            obj.filesErr(ii) = [];
+                            filesSavedErr{ii} = [];
+                            obj.files(end+1) = obj2.files(jj);
+                            if ~obj2.files(jj).isdir
+                                obj.nfiles = obj.nfiles+1;
+                            end
+                            filesSaved{end+1} = obj2.files(jj).name;        %#ok<AGROW>
+                            obj.logger.Write('%d. Mismatch in dates for saved data in %s. Will re-validate this file.\n', jj, obj.filesErr(ii).name);
+                            obj.changed = true;                            
+                        end
+                        break
+                    end
+                end
+            end
+            
+                        
+            %%%% Now check saved object nonerror list
+            N = length(obj.files);
+            for ii = N:-1:1
+	            found = false;            
+                for jj = 1:length(obj2.files)
+                    if strcmp(obj.files(ii).name, obj2.files(jj).name)
+                        found = true;
+                        if obj.files(ii).datenum ~= obj2.files(jj).datenum
+                            obj.files(ii) = obj2.files(jj).copy();
+                            obj.logger.Write('%d. Mismatch in dates for saved data in %s. Will re-validate this file.\n', jj, obj.files(ii).name);
+                            obj.changed = true;
+                        end
+                        break
+                    end
+                end
+                if ~found
+                    obj.logger.Write('%d. Deleting file %s from saved data.\n', ii, obj.files(ii).name);
+                    if ~obj.files(ii).isdir
+                        obj.nfiles = obj.nfiles-1;
+                    end
+                    obj.files(ii) = [];
+                    filesSaved(ii) = [];                    
+                    obj.changed = true;
+                end     
+            end
+            
+            % We don't need to check if files were added if there's no
+            % change and the number of current files found is same as saved
+            % number of files. Exist early. 
+            if obj.changed == false
+                if (length(obj.files) + length(obj.filesErr)) == length(obj2.files)
+                    return
+                end
+            end
+            
+            
+            %%%% Find which files were added to current dataset and add
+            %%%% them to saved object
+            for ii = 1:length(obj2.files)
+                found = false;
+                for jj = 1:length(obj.files)
+                    if strcmp(obj2.files(ii).name, obj.files(jj).name)
+                        found = true;
+                        break
+                    end
+                end
+                
+                % If the current file filesCurr(ii) is not found in saved dataset then it's a new file. 
+                % Add new file from current dataset to saved dataset
+                if ~found
+                    obj.files(end+1) = obj2.files(ii);
+                    if ~obj2.files(ii).isdir
+                        obj.nfiles = obj.nfiles+1;
+                    end
+                    filesSaved{end+1} = obj2.files(ii).name; %#ok<AGROW>
+                    obj.logger.Write('%d. Adding new file %s to saved data.\n', jj, obj2.files(ii).name);
+                    obj.changed = true;
+                end
+            end
+
+            
+            % Sort non-error files
+            if obj.changed
+                [~, order] = sort(filesSaved);
+                obj.files = obj.files(order);
+                obj.logger.Write('\n');
+            end
+            
+        end
+        
+
+ 
+        % -----------------------------------------------------------------------------------
+        function SaveDataSet(obj)
+            if ~obj.changed
+                return
+            end
+            obj.changed = false;
+            obj.logger.Write('Saving DataFilesClass object in %s\n', obj.savedFilename);
+            if ~exist(obj.config.DerivedDataRootDir, 'dir')
+                mkdir(obj.config.DerivedDataRootDir)
+            end
+            save(obj.savedFilename, '-mat', 'obj');
+        end
+        
+
+ 
+        % -----------------------------------------------------------------------------------
+        function b = LoadDataSet(obj)
+            b = false;
+            if ~exist(obj.savedFilename, 'file')
+                return
+            end
+            dataset = load(obj.savedFilename);
+            dataset.obj.Compare();
+            if isempty(dataset.obj.files) && isempty(dataset.obj.filesErr)
+                return
+            end
+            obj.logger.Write('Loaded saved DataFilesClass object:    %s\n', obj.savedFilename);
+            obj.Copy(dataset.obj);
+
+            obj.ErrorCheck()
+            obj.ErrorCheckName()
+            b = true;
         end
         
         
         
         % -----------------------------------------------------------------------------------
-        function GetDataSet(obj)
+        function GetDataSet(obj, skipsaveddata)
+            if ~exist('skipsaveddata','var')
+                skipsaveddata = '';
+            end
+            if ~strcmp(skipsaveddata, 'skipsaveddata')
+                if obj.LoadDataSet()
+                    return
+                end
+            end
+
+            obj.changed = true;
+
+            mode = '';            
+            if obj.dirFormats.type == 0
+                iFormats = 1:length(obj.dirFormats.choices);
+            else
+                iFormats = obj.dirFormats.type;
+                mode = 'noerrorcheck';
+            end
             if exist(obj.rootdir, 'dir')~=7
                 error('Invalid subject folder: ''%s''', obj.rootdir);
             end
-            
-            for ii = 1:length(obj.dirFormats.choices)
+            for ii = iFormats
                 obj.InitLookupTable();
                 obj.FindDataSet(ii);
                 
+                if strcmp(mode,'noerrorcheck')
+                    continue
+                end
+                if isempty(obj.files)
+                    continue
+                end                
+                
                 % Remove any files that cannot pass the basic test of loading
                 % its data
                 obj.ErrorCheck();
                 if ~isempty(obj.files)
                     break
                 end
+            end            
+            if strcmp(mode,'noerrorcheck')
+                return
             end
-            
             obj.ErrorCheckFinal();
         end
 
         
         
         % -----------------------------------------------------------------------------------
-        function InitFolderFromats(obj)
+        function InitFolderFormats(obj)
             obj.dirFormats.choices = {
                 
                 %%%% 1. Flat #1
                 {
-                ['*_run*.', obj.filetype];
+                    ['*_run*.', obj.filetype];
                 }
 
                 %%%% 2. Flat #2
                 {
-                ['*.', obj.filetype];
+                    ['*.', obj.filetype];
                 }
                 
                 %%%% 3. Subjects 
                 {
-                '*';
-                ['*.', obj.filetype];
+                    '*';
+                    ['*.', obj.filetype];
                 }
 
                 %%%% 4. BIDS #1,    sub-<label>[_ses-<label>][_task-<label>][_run-<index>]_nirs.snirf
                 {                
-                'sub-*';
-                'ses-*';
-                ['nirs/sub-*_run-*_nirs.', obj.filetype];
+                    'sub-*';
+                    'ses-*';
+                    ['nirs/sub-*_run-*_nirs.', obj.filetype];
                 }
                                 
                 %%%% 5. BIDS #2 
                 {
-                'sub-*';
-                ['nirs/sub-*_run-*_nirs.', obj.filetype];
+                    'sub-*';
+                    ['nirs/sub-*_run-*_nirs.', obj.filetype];
                 }
                 
                 %%%% 6. BIDS #3
                 {
-                'sub-*';
-                ['nirs/sub-*_*_nirs.', obj.filetype];
+                    'sub-*';
+                    ['nirs/sub-*_*_nirs.', obj.filetype];
                 }
                 
                 %%%% 7. BIDS #4 
                 {
-                '*';
-                ['nirs/sub-*_*_nirs.', obj.filetype];
+                    '*';
+                    ['nirs/sub-*_*_nirs.', obj.filetype];
                 }
                 
                 %%%% 8. BIDS folder structure
                 {
-                'sub-*';
-                'ses-*';
-                ['nirs/sub-*_run-*_nirs.', obj.filetype];
+                    'sub-*';
+                    'ses-*';
+                    ['nirs/sub-*_run-*_nirs.', obj.filetype];
                 }
                                
                 %%%% 9. BIDS-like folder structure without file naming restrictions
                 {
-                'sub-*';
-                'ses-*';
-                ['nirs/*.', obj.filetype];
+                    'sub-*';
+                    'ses-*';
+                    ['nirs/*.', obj.filetype];
                 }
                                 
                 %%%% 10. BIDS-like folder structure without nirs sub-folder
                 {
-                'sub-*';
-                'ses-*';
-                ['*.', obj.filetype];
+                    'sub-*';
+                    'ses-*';
+                    ['*.', obj.filetype];
                 }
-                                
-                };
+                
+            };
         end
         
 
@@ -254,7 +436,12 @@ function FindDataSet(obj, iFormat, iPattern, parentdir)
                     continue;
                 end
                 if dirs(ii).IsFile()
-                    if strcmp(pattern, '*')
+                    % If the pattern that found this file has no extension
+                    % that means this pattern is meant ONLY for folders. 
+                    % Therefore any file matches should be skipped. 
+                    % Bug Fix - JD, Jun 20, 2023
+                    [~,~,ext] = fileparts(pattern);
+                    if isempty(ext)  
                         continue
                     end
                     if includes(dirs(ii).name, obj.filetype)
@@ -290,7 +477,7 @@ function AddParentDirs(obj, dirname)
                     continue;
                 end
                 obj.files(end+1) = FileClass([obj.rootdir, '/', pathrel2], obj.rootdir);
-                obj.AddLookupTable(obj.files(end).name)
+                obj.AddLookupTable(pathrel2)
             end
         end
         
@@ -298,9 +485,9 @@ function AddParentDirs(obj, dirname)
         
         % ----------------------------------------------------
         function AddFile(obj, file)
-            obj.files(end+1) = file;
+            obj.files(end+1) = FileClass();
+            obj.files(end).Add(file);
             obj.nfiles = obj.nfiles+1;
-            obj.AddLookupTable(obj.files(end).filename)
         end
         
         
@@ -308,6 +495,9 @@ function AddFile(obj, file)
         % ----------------------------------------------------
         function ErrorCheckName(obj)
             for ii = length(obj.files):-1:1
+                if ~obj.files(ii).IsValidated()
+                    continue;
+                end
                 if obj.files(ii).ErrorCheckName()<0
                     q = obj.AskToFixNameConflicts(ii);                    
                     if q == 1
@@ -411,7 +601,7 @@ function ErrorCheckFinal(obj)
             msg{3} = sprintf('it may cause incorrect results in processing. Do you want to rename this %s?', filetype);
             errmsg = [msg{:}];            
         end
-        
+
 
         
         % -------------------------------------------------------
@@ -494,6 +684,10 @@ function ErrorCheck(obj)
                 if obj.files(ii).isdir
                     continue;
                 end
+                if obj.files(ii).IsValidated()
+                    dataflag = true;
+                    continue;
+                end
                 filename = [obj.rootdir, obj.files(ii).name]; %#ok<NASGU>
                 eval( sprintf('o = %s(filename);', constructor) );
                 if  o.GetError() < 0
@@ -513,12 +707,15 @@ function ErrorCheck(obj)
                 end
                 if ~isempty(o.GetErrorMsg())
                     obj.files(ii).SetError(o.GetErrorMsg());
+                else
+                    obj.files(ii).SetValid();                    
                 end
                 hwait = waitbar_improved(ii/length(obj.files), hwait, msg);
             end
             
             if dataflag==false
                 obj.files = FileClass.empty();
+                obj.nfiles = 0;
             else
                 for jj = 1:length(errorIdxs)
                     obj.filesErr(end+1) = obj.files(errorIdxs(jj)).copy;
@@ -527,6 +724,9 @@ function ErrorCheck(obj)
                 obj.nfiles = obj.nfiles - length(errorIdxs);
             end
             close(hwait);
+
+            obj.logger.Write('\n');
+            
         end
         
         
@@ -539,6 +739,7 @@ function ErrorCheck(obj)
             err = obj.err;
         end
                 
+        
         % ----------------------------------------------------------
         function PrintFolderStructure(obj, options)
             if ~exist('options','var')
@@ -569,14 +770,56 @@ function PrintFolderStructure(obj, options)
                     if optionExists(options, 'numbered')
                         n = k*stepsize+stepsize+j;
                         obj.logger.Write(sprintf('%d.%s%s\n', ii, blanks(n), obj.files(ii).filename));
+                        c = c+1;
                     else
                         n = k*stepsize+stepsize;
-                        obj.logger.Write(sprintf('%s%s\n', blanks(n), obj.files(ii).filename));
+                        obj.logger.Write(sprintf('%d.%s%s\n', ii, blanks(n), obj.files(ii).filename));
                     end
                 end
             end
             obj.logger.Write('\n');
         end
+
+        
+        
+        % ----------------------------------------------------------
+        function Copy(obj, obj2)
+            obj.files       = obj2.files.copy();
+            obj.filesErr    = obj2.filesErr.copy();
+            obj.filetype    = obj2.filetype;
+            obj.nfiles      = obj2.nfiles;
+            obj.err         = obj2.err;
+            obj.errmsg      = obj2.errmsg;
+            obj.dirFormats  = obj2.dirFormats;
+            obj.changed     = obj2.changed;
+        end
+        
+        
+        
+        % ----------------------------------------------------------
+        function b = eq(obj, obj2)
+            b = false;
+            if length(obj.files) ~= length(obj2.files)
+                return
+            end
+            if length(obj.filesErr) ~= length(obj2.filesErr)
+                return
+            end
+            for ii = 1:length(obj.files)
+                if obj.files ~= obj2.files
+                    return
+                end
+            end
+            for ii = 1:length(obj.filesErr)
+                if obj.filesErr ~= obj2.filesErr
+                    return
+                end
+            end
+            if obj.filetype ~= obj2.filetype
+                return
+            end
+            b = true;
+        end
         
     end
         
@@ -608,6 +851,18 @@ function AddLookupTable(obj, str)
             b = obj.lookupTable(string2hash(str, n));
         end
         
+
         
+        % ----------------------------------------------------------
+        function Print(obj)
+            for ii = 1:length(obj.files)
+                if obj.files(ii).isdir
+                    continue
+                end
+                fprintf('%s\n', obj.files(ii).name);                
+            end
+        end
+               
     end
+                        
 end
diff --git a/DataTree/AcquiredData/DataFiles/FileClass.m b/DataTree/AcquiredData/DataFiles/FileClass.m
index 66d0f8dc..ea2e1f23 100644
--- a/DataTree/AcquiredData/DataFiles/FileClass.m
+++ b/DataTree/AcquiredData/DataFiles/FileClass.m
@@ -20,6 +20,7 @@
         err
         logger
         errmsg
+        errcodeUnvalidated
     end
     
     methods
@@ -40,8 +41,8 @@
             obj.rootdir    = '';
             obj.errmsg     = '';          % Assume file is not loadable
             obj.logger     = InitLogger(logger);            
-            obj.errmsg    = '';
             
+            obj.errcodeUnvalidated = -9999;
             
             if nargin==0
                 return;
@@ -88,7 +89,13 @@ function Add(obj, obj2)
             obj.filename     = obj2.name;
             obj.name         = getPathRelative(rootpath, obj2.rootdir);
             obj.rootdir 	 = obj2.rootdir;
-            obj.err          = 0;          % Set error to NO ERROR            
+            obj.date         = obj2.date;
+            obj.datenum      = datestr2datenum(obj.date);
+            if obj.IsFile()
+                obj.err          = obj.errcodeUnvalidated;          % Set error to unvalidated
+            else
+                obj.err          = 0;
+            end
         end
         
         
@@ -434,7 +441,7 @@ function Loaded(obj)
             if isempty(obj.name)
                 return;
             end
-            if obj.err ~= 0
+            if (obj.err ~= 0)  &&  (obj.IsValidated)
                 return;
             end
             b = false;            
@@ -443,19 +450,22 @@ function Loaded(obj)
         
         % ----------------------------------------------------
         function err = ErrorCheckName(obj)
+            err = 0;
             [p1,f1] = fileparts(obj.name);
             [p2,f2] = fileparts(filesepStandard(obj.rootdir,'nameonly:file'));
             [~,f3]  = fileparts(p2);
             if strcmp(f1, p1)
-                obj.err = -1;
+                err = -1;
             end
             if strcmp(f1, f2)
-                obj.err = -2;
+                err = -2;
             end
             if strcmp(f1, f3)
-                obj.err = -3;
+                err = -3;
+            end
+            if err ~= 0
+                obj.err = err;
             end
-            err = obj.err;
         end
         
         
@@ -545,6 +555,7 @@ function NameConflictFixed(obj)
         % -----------------------------------------------------
         function SetError(obj, errmsg)
             obj.errmsg = errmsg;
+            obj.err = -1;
         end
         
         
@@ -553,6 +564,62 @@ function SetError(obj, errmsg)
             msg = obj.errmsg;
         end
         
+        
+        % -----------------------------------------------------
+        function SetValid(obj)
+            obj.err = 0;
+        end
+
+        
+        % -----------------------------------------------------
+        function SetErrorUnvalidated(obj)
+            obj.err = obj.errcodeUnvalidated;
+        end
+
+        
+        % -----------------------------------------------------
+        function b = IsValidated(obj)
+            b = false;
+            if obj.err == obj.errcodeUnvalidated
+                return 
+            end
+            b = true;
+        end
+        
+        
+        
+        % -----------------------------------------------------
+        function b = IsUnValidated(obj)
+            b = true;
+            if obj.err == obj.errcodeUnvalidated
+                return 
+            end
+            b = false;
+        end
+
+        
+        
+        % ----------------------------------------------------------
+        function b = eq(obj, obj2)
+            b = false;
+            if ~strcmp(obj.name, obj2.name)
+                return;
+            end
+            if ~strcmp(obj.date, obj2.date)
+                return;
+            end
+            if obj.isdir ~= obj2.isdir
+                return;
+            end
+            if obj.bytes ~= obj2.bytes
+                return;
+            end
+            if obj.datenum ~= obj2.datenum
+                return;
+            end
+            b = true;
+        end
+        
     end
   
 end
\ No newline at end of file
diff --git a/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_DatasetLoad.m b/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_DatasetLoad.m
index b1cec309..3958250c 100644
--- a/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_DatasetLoad.m
+++ b/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_DatasetLoad.m
@@ -8,8 +8,20 @@
 end
 
 try
-    dsetid = H5D.open(gid, name);
 
+    if ~H5L.exists(gid, name, 'H5P_DEFAULT')
+        switch(class(val0))
+            case 'char'
+                val = '';
+            case 'cell'
+                val = {};
+            otherwise
+                val = [];
+        end
+        return
+    end
+
+    dsetid = H5D.open(gid, name);
     % NOTE: HDF5 stores contiguous muti-dimensional arrays in row-major order.
     % Matlab stores them in row-major order. We want to transpose the loaded data 
     % it back to Matlab's column-major storage order and thus get back the 
@@ -29,3 +41,5 @@
             val = [];
     end
 end
+
+
diff --git a/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_Transpose.m b/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_Transpose.m
index 034735be..b1a66613 100644
--- a/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_Transpose.m
+++ b/DataTree/AcquiredData/DataFiles/Hdf5/HDF5_Transpose.m
@@ -1,16 +1,16 @@
-function val = HDF5_Transpose(val, options)
+function val = HDF5_Transpose(val, option)
 
-if ~exist('options','var')
-    options = '';
+if ~exist('option','var')
+    option = '';
 end
 
 % Matlab stores contiguous muti-dimensional arrays in column-major order.
 % HDF5 stores them in row-major order. We want to transpose the data to agree
 % with the file format's storage order.
 if (~isrow(val) && ~iscolumn(val))              || ...      % Matrices
-     ~isempty(findstr('multidim', options))     || ...      % Force multi-dimensional even if vector
-     ~isempty(findstr('2D', options))           || ...      % Force 2D even if vector
-     ~isempty(findstr('3D', options))                       % Force 3D even if vector
+     ~isempty(findstr('multidim', option))     || ...      % Force multi-dimensional even if vector
+     ~isempty(findstr('2D', option))           || ...      % Force 2D even if vector
+     ~isempty(findstr('3D', option))                       % Force 3D even if vector
     
     val = permute(val, ndims(val):-1:1);
     
diff --git a/DataTree/AcquiredData/DataFiles/Hdf5/hdf5write_safe.m b/DataTree/AcquiredData/DataFiles/Hdf5/hdf5write_safe.m
index 61eba6f0..dbd3a09a 100644
--- a/DataTree/AcquiredData/DataFiles/Hdf5/hdf5write_safe.m
+++ b/DataTree/AcquiredData/DataFiles/Hdf5/hdf5write_safe.m
@@ -1,18 +1,18 @@
 function err = hdf5write_safe(fileobj, name, val, options)
 err = -1;
-if isempty(val)
+if ~exist('val','var')
     return;
 end
 if ~exist('options','var')
     options = '';
 end
 
-force_scalar = false;
 force_array = false;
+force_vector = false;
 if any(strcmp(options, 'array'))
     force_array = true;
-elseif any(strcmp(options, 'scalar'))
-    force_scalar = true;
+elseif any(strcmp(options, 'vector'))
+    force_vector = true;
 end
 
 % Identify type of val and use SNIRF v1.1-compliant write function
@@ -32,25 +32,21 @@
 
 % Create dataset
 if iscell(val) || isstring(val)
-    if length(val) > 1 && ~force_scalar || force_array  % Returns true for single strings, believe it or not
+    if (length(val) > 1) || force_array  % Returns true for single strings, believe it or not
         write_string_array(fid, name, val);
     else
         write_string(fid, name, val);
     end
 elseif ischar(val)
     write_string(fid, name, val);
-elseif isfloat(val)
-    if length(val) > 1 && ~force_scalar || force_array
+elseif isfloat(val) || isinteger(val)
+    if force_vector
+        write_numeric_vector(fid, name, val);
+    elseif (length(val) > 1) || force_array
         write_numeric_array(fid, name, val);
     else
         write_numeric(fid, name, val);
     end
-elseif isinteger(val)
-    if length(val) > 1 && ~force_scalar || force_array
-        write_numeric_array(fid, name, val);  % As of now, no integer arrays exist
-    else
-        write_integer(fid, name, val);
-    end
 else
     warning(['An unrecognized variable was saved to ', name])
 end
@@ -82,38 +78,63 @@
 
 
 
+
 % -----------------------------------------------------------------
 function err = write_numeric(fid, name, val)
-tid = H5T.copy('H5T_NATIVE_DOUBLE');
+warning off;  % Suppress the int truncation warning
+hdftype = convert_to_hdf5_type(val);
+tid = H5T.copy(hdftype);
 sid = H5S.create('H5S_SCALAR');
 dsid = H5D.create(fid, name, tid, sid, 'H5P_DEFAULT');
 H5D.write(dsid, tid, 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', val);
 err = 0;
+warning on;
+
 
 
 
 % -----------------------------------------------------------------
-function err = write_integer(fid, name, val)
-warning off;  % Suppress the int truncation warning
-tid = H5T.copy('H5T_NATIVE_ULONG');
-sid = H5S.create('H5S_SCALAR');
-dsid = H5D.create(fid, name, tid, sid, 'H5P_DEFAULT');
-H5D.write(dsid, tid, 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', int32(val));
+function err = write_numeric_array(fid, name, data)
 err = 0;
-warning on;
+data = HDF5_Transpose(data, '2D');
+sizedata = size(data);
+n = sizedata;
+
+tid = -1;
+sid = -1;
+gid = -1;
+dsid = -1;
+
+maxdims = n;
+
+hdftype = convert_to_hdf5_type(data);
+
+try
+    
+    sid = H5S.create_simple(numel(n), fliplr(n), fliplr(maxdims));
+    gid = HDF5_CreateGroup(fid, fileparts(name));
+    dsid = H5D.create(gid, name, hdftype, sid, 'H5P_DEFAULT');
+    H5D.write(dsid, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', data);
+    
+catch
+    
+    % Clean up; Close everything
+    cleanUp(tid, sid, gid, dsid);
+    err = -1;
+    return;
+    
+end
+cleanUp(tid, sid, gid, dsid);
+
 
 
 
 % -----------------------------------------------------------------
-function err = write_numeric_array(fid, name, data)
+function err = write_numeric_vector(fid, name, data)
 err = 0;
-data = HDF5_Transpose(data);
-sizedata = size(data);
-if sizedata(1) == 1 || sizedata(2) == 1
-    n = length(data);
-else
-    n = sizedata;
-end
+data = HDF5_Transpose(data, '2D');
+sizedata = length(data);
+n = sizedata;
 
 tid = -1;
 sid = -1;
@@ -121,11 +142,12 @@
 dsid = -1;
 
 maxdims = n;
+hdftype = convert_to_hdf5_type(data);
 try
     
-    sid = H5S.create_simple(numel(n), fliplr(n), fliplr(maxdims));   
+    sid = H5S.create_simple(numel(n), fliplr(n), fliplr(maxdims));
     gid = HDF5_CreateGroup(fid, fileparts(name));
-    dsid = H5D.create(gid, name, 'H5T_NATIVE_DOUBLE', sid, 'H5P_DEFAULT');
+    dsid = H5D.create(gid, name, hdftype, sid, 'H5P_DEFAULT');
     H5D.write(dsid, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', data);
     
 catch
@@ -140,6 +162,7 @@
 
 
 
+
 % ------------------------------------------------------
 function cleanUp(tid, sid, gid, dsid)
 if ~isnumeric(tid)
@@ -156,3 +179,30 @@ function cleanUp(tid, sid, gid, dsid)
 end
 
 
+
+
+% -------------------------------------------------------
+function hdftype = convert_to_hdf5_type(val)
+hdftype = '';
+switch class(val)
+    case 'int8'
+        hdftype = 'H5T_STD_8LE';
+    case 'uint8'
+        hdftype = 'H5T_STD_U8LE';
+    case 'int16'
+        hdftype = 'H5T_STD_16LE';
+    case 'uint16'
+        hdftype = 'H5T_STD_U16LE';
+    case 'int32'
+        hdftype = 'H5T_STD_32LE';
+    case 'uint32'
+        hdftype = 'H5T_STD_U32LE';
+    case 'int64'
+        hdftype = 'H5T_STD_64LE';
+    case 'uint64'
+        hdftype = 'H5T_STD_U64LE';
+    case 'double'
+        hdftype = 'H5T_NATIVE_DOUBLE';
+end
+
+
diff --git a/DataTree/AcquiredData/Nirs/NirsClass.m b/DataTree/AcquiredData/Nirs/NirsClass.m
index f9aa6b43..a6156a17 100644
--- a/DataTree/AcquiredData/Nirs/NirsClass.m
+++ b/DataTree/AcquiredData/Nirs/NirsClass.m
@@ -76,7 +76,7 @@
             if ~ischar(varargin{1}) && nargin==1
                 obj.ErrorCheck();
                 return
-            end
+            end           
             
             if nargin==2
                 obj.SetDataStorageScheme(varargin{2});
@@ -213,7 +213,7 @@ function SortData(obj)
         % ---------------------------------------------------------
         function err = LoadTime(obj, fdata, err)
             if nargin == 1
-                fdata = load(obj.GetFilename(),'-mat', 't');
+                load(obj.GetFilename(),'-mat', 't');
                 err = 0;
                 return
             elseif nargin == 2
@@ -238,8 +238,12 @@ function SortData(obj)
         
         % ---------------------------------------------------------
         function err = LoadProbeMeasList(obj, fdata, err)
-            if isproperty(fdata,'SD')
-                obj.CopyProbe(fdata.SD);
+            SD = fdata;
+            while isproperty(SD,'SD')
+                SD = SD.SD;
+            end
+            if isproperty(SD,'SrcPos')
+                obj.CopyProbe(SD);
                 if isempty(obj.SD) && err == 0
                     err = -4;
                 end
@@ -361,10 +365,15 @@ function SaveMat(obj, fname, options)
             end
             [~, k1] = sortrows(obj.SD.MeasList);
             [~, k2] = sortrows(obj2.SD.MeasList);
-            if ~all(obj.SD.MeasList(k1,:) == obj2.SD.MeasList(k2,:))
+            
+            % Compare meas list dimensions first - if that's not equal then the probes aren't equal
+            if ~all( size(obj.SD.MeasList(k1,:)) == size(obj2.SD.MeasList(k2,:)) )
                 return;
-            end            
+            end
 
+            if ~all(obj.SD.MeasList(k1,:) == obj2.SD.MeasList(k2,:))
+                return;
+            end
             
             for kk = 1:length(fields)
                 for jj = 1:length(fields{kk})
@@ -412,6 +421,54 @@ function SaveMat(obj, fname, options)
         
         
         
+        % -------------------------------------------------------
+        function b = ProbeSimilar(obj, obj2)
+            b = false;
+                                              
+            % MeasList
+            [~, k1] = sortrows(obj.SD.MeasList);
+            [~, k2] = sortrows(obj2.SD.MeasList);
+            if ~all( size(obj.SD.MeasList) == size(obj2.SD.MeasList) )
+                return;
+            end
+            if ~all(obj.SD.MeasList(k1,:) == obj2.SD.MeasList(k2,:))
+                return;
+            end
+            
+            if ~all( size(obj.SD.SrcPos) == size(obj2.SD.SrcPos) )
+                return;
+            end
+            if ~all(obj.SD.SrcPos == obj2.SD.SrcPos)
+                return;
+            end
+            
+            if ~all( size(obj.SD.SrcPos3D) == size(obj2.SD.SrcPos3D) )
+                return;
+            end
+            if ~all( obj.SD.SrcPos3D == obj2.SD.SrcPos3D )
+                return;
+            end
+                        
+            if ~all( size(obj.SD.DetPos) == size(obj2.SD.DetPos) )
+                return;
+            end
+            if ~all( obj.SD.DetPos == obj2.SD.DetPos )
+                return;
+            end
+            
+            if ~all( size(obj.SD.DetPos3D) == size(obj2.SD.DetPos3D) )
+                return;
+            end
+            if ~all( obj.SD.DetPos3D == obj2.SD.DetPos3D )
+                return;
+            end
+            
+            
+            b = true;
+        end
+        
+        
+        
         % --------------------------------------------------------------------
         function b = EqualStim(obj, obj2)
             b = false;
@@ -460,8 +517,7 @@ function SaveMat(obj, fname, options)
             end
             if ~all(obj.aux(:) == obj2.aux(:))
                 return;
-            end
-            
+            end            
             B = true;
         end
         
@@ -612,7 +668,7 @@ function SetDataTimeSeries(obj, val)
                     ml(ii).detectorIndex = obj.SD.MeasList(ii,2);
                     ml(ii).dataTypeIndex = 0;
                     ml(ii).wavelengthIndex = obj.SD.MeasList(ii,4);
-            end
+                end
             end
         end
         
@@ -1261,7 +1317,7 @@ function FixProbeSpatialUnit(obj)
         
         
         % ----------------------------------------------------------------------------------
-        function CopyProbe(obj, SD)             %#ok<INUSD>
+        function CopyProbe(obj, SD)
             fields = propnames(obj.SD);
             for ii = 1:length(fields)
                 if eval( sprintf('isfield(SD, ''%s'')', fields{ii}) )
@@ -1339,7 +1395,7 @@ function CopyProbe(obj, SD)             %#ok<INUSD>
                                     
             % MesListAct
             if size(obj.SD.MeasListAct,1) < size(obj.SD.MeasList,1)
-                d = size(obj.SD.MeasListAct,1) - size(obj.SD.MeasList,1);
+                d = size(obj.SD.MeasListAct,1) - size(obj.SD.MeasList,1); %#ok<*PROPLC>
                 if d < 0
                     obj.SD.MeasListAct(end+1:end+abs(d)) = ones(abs(d),1);
                 elseif d>1
@@ -1376,33 +1432,121 @@ function CopyStruct(obj, s)
                 end
             end
         end
-
+        
         
         
         % ----------------------------------------------------------------------------------
-        function ConvertSnirfProbe(obj, snirf)
-            obj.SD.Lambda = snirf.probe.wavelengths;
-            obj.SD.SrcPos = snirf.probe.sourcePos2D;
-            obj.SD.DetPos = snirf.probe.detectorPos2D;
-            obj.SD.SrcPos3D = snirf.probe.sourcePos3D;
-            obj.SD.DetPos3D = snirf.probe.detectorPos3D;
-            obj.SD.MeasList = snirf.GetMeasList();
-            obj.SD.SpatialUnit = snirf.GetLengthUnit();
-            if length(snirf.probe.landmarkLabels) == size(snirf.probe.landmarkPos3D,1)
-                obj.SD.Landmarks3D.labels   = snirf.probe.landmarkLabels;
-            	obj.SD.Landmarks3D.pos      = snirf.probe.landmarkPos3D;
-            end
-            if length(snirf.probe.landmarkLabels) == size(snirf.probe.landmarkPos2D,1)
-            	obj.SD.Landmarks2D.labels   = snirf.probe.landmarkLabels;
-                obj.SD.Landmarks2D.pos      = snirf.probe.landmarkPos2D;
-        	end
-            if     ~isempty(obj.SD.Landmarks3D.labels)
-                obj.SD.Landmarks.pos        = obj.SD.Landmarks3D.pos;
-                obj.SD.Landmarks.labels     = obj.SD.Landmarks3D.labels;
-            elseif ~isempty(obj.SD.Landmarks2D.labels)
-                obj.SD.Landmarks.pos        = obj.SD.Landmarks2D.pos;
-                obj.SD.Landmarks.labels     = obj.SD.Landmarks2D.labels;
+        function b = IsProbeFlat(obj)
+            if all(size(obj.SD.SrcPos) == size(obj.SD.SrcPos3D))
+                if ~all(obj.SD.SrcPos == obj.SD.SrcPos3D)
+                    b = false;
+                    return;
+                end
+            end
+            if all(size(obj.SD.DetPos) == size(obj.SD.DetPos3D))
+                if ~all(obj.SD.DetPos == obj.SD.DetPos3D)
+                    b = false;
+                    return
+                end
+            end
+            optpos = [obj.SD.SrcPos3D; obj.SD.DetPos3D];
+            ncoord = size(optpos, 2);
+            for ii = 1:ncoord
+                if length(unique(optpos(:,ii))) == 1
+                    b = true;
+                    return
+                end
+            end
+            if isempty(optpos)
+                b = true;
+                return
             end                
+            optpos = [obj.SD.SrcPos; obj.SD.DetPos];
+            ncoord = size(optpos, 2);
+            flag = zeros(1, ncoord);
+            for ii = 1:ncoord
+                if length(unique(optpos(:,ii))) > 1
+                    flag(ii) = 1;
+                end
+            end
+            if sum(flag) == ncoord
+                b = false;
+                return
+            end
+            b = true;
+        end
+        
+        
+        
+        % ----------------------------------------------------------------------------------
+        function b = Are3DLandmarksFlat(obj)
+            b = true;
+            if isempty(obj.SD.Landmarks3D.pos)
+                return;
+            end
+            ncoord = size(obj.SD.Landmarks3D.pos, 2);
+            for ii = 1:ncoord
+                if length(unique(obj.SD.Landmarks3D.pos(:,ii)))==1
+                    return;
+                end
+            end
+            b = false;
+        end
+        
+
+        
+        % ----------------------------------------------------------------------------------
+        function b = AreLandmarksFlat(obj)
+            b = true;
+            if isempty(obj.SD.Landmarks.pos)
+                return;
+            end
+            ncoord = size(obj.SD.Landmarks.pos, 2);
+            for ii = 1:ncoord
+                if length(unique(obj.SD.Landmarks.pos(:,ii)))==1
+                    return;
+                end
+            end
+            b = false;
+        end
+        
+
+        
+        % ----------------------------------------------------------------------------------
+        function err = ConvertSnirfProbe(obj, snirf)
+            err = 0;
+            if isempty(snirf)
+                return
+            end
+            if isempty(snirf.probe)
+                return
+            end
+            try
+	            obj.SD.Lambda = snirf.probe.wavelengths;
+	            obj.SD.SrcPos = snirf.probe.sourcePos2D;
+	            obj.SD.DetPos = snirf.probe.detectorPos2D;                
+	            obj.SD.SrcPos3D = snirf.probe.sourcePos3D;
+	            obj.SD.DetPos3D = snirf.probe.detectorPos3D;
+	            obj.SD.MeasList = snirf.GetMeasList();
+	            obj.SD.SpatialUnit = snirf.GetLengthUnit();
+	            if length(snirf.probe.landmarkLabels) == size(snirf.probe.landmarkPos3D,1)
+	                obj.SD.Landmarks3D.labels   = snirf.probe.landmarkLabels;
+	            	obj.SD.Landmarks3D.pos      = snirf.probe.landmarkPos3D;
+	            end
+                if length(snirf.probe.landmarkLabels) == size(snirf.probe.landmarkPos2D,1)
+                    obj.SD.Landmarks2D.labels   = snirf.probe.landmarkLabels;
+                    obj.SD.Landmarks2D.pos      = snirf.probe.landmarkPos2D;
+                end
+                if     ~isempty(obj.SD.Landmarks3D.labels)
+                    obj.SD.Landmarks.pos        = obj.SD.Landmarks3D.pos;
+                    obj.SD.Landmarks.labels     = obj.SD.Landmarks3D.labels;
+                elseif ~isempty(obj.SD.Landmarks2D.labels)
+                    obj.SD.Landmarks.pos        = obj.SD.Landmarks2D.pos;
+                    obj.SD.Landmarks.labels     = obj.SD.Landmarks2D.labels;
+                end
+            catch
+                err = -1;
+            end
         end
         
         
@@ -1441,10 +1585,7 @@ function ConvertSnirfStim(obj, snirf)
         
         % ----------------------------------------------------------------------------------
         function ConvertSnirfAux(obj, snirf)
-            obj.aux = zeros(length(obj.t), length(snirf.aux));
-            for ii = 1:length(snirf.aux)
-                obj.aux(:,ii) = snirf.aux(ii).dataTimeSeries;
-            end
+            obj.aux = snirf.GetAuxDataMatrix();
         end
         
         
@@ -1530,7 +1671,19 @@ function ErrorCheck(obj)
                     end
                 end
             end
-            
+            if obj.IsProbeFlat()
+                obj.SD.SrcPos3D = [];
+                obj.SD.DetPos3D = [];
+            end
+            if obj.Are3DLandmarksFlat() && ~obj.AreLandmarksFlat()
+                obj.SD.Landmarks3D.pos      = obj.SD.Landmarks.pos;
+                obj.SD.Landmarks3D.labels   = obj.SD.Landmarks.labels;
+            end
+            for ii = 1:length(obj.SD.AnchorList)
+                if ischar(obj.SD.AnchorList{ii,1})
+                    obj.SD.AnchorList{ii,1} = str2num(obj.SD.AnchorList{ii,1});
+                end
+            end
         end
         
         
diff --git a/DataTree/AcquiredData/Snirf/AuxClass.m b/DataTree/AcquiredData/Snirf/AuxClass.m
index ab812d24..95bbb24f 100644
--- a/DataTree/AcquiredData/Snirf/AuxClass.m
+++ b/DataTree/AcquiredData/Snirf/AuxClass.m
@@ -144,8 +144,8 @@
             
             hdf5write_safe(fid, [location, '/name'], obj.name);
             hdf5write_safe(fid, [location, '/dataTimeSeries'], obj.dataTimeSeries, 'array');
-            hdf5write_safe(fid, [location, '/time'], obj.time, 'array');
-            hdf5write_safe(fid, [location, '/timeOffset'], obj.timeOffset, 'array');
+            hdf5write_safe(fid, [location, '/time'], obj.time, 'vector');
+            hdf5write_safe(fid, [location, '/timeOffset'], obj.timeOffset, 'vector');
         end
         
         
@@ -179,16 +179,14 @@ function SetDataTimeSeries(obj, val)
         % ----------------------------------------------------------------------------------
         function Copy(obj, obj2)
             if isempty(obj)
-                obj = DataClass();
+                obj = AuxClass();
             end
-            if ~isa(obj2, 'DataClass')
+            if ~isa(obj2, 'AuxClass')
                 return;
             end
-            for ii=1:length(obj2.measurementList)
-                obj.measurementList(ii) = obj2.measurementList(ii).copy();      % shallow copy ok because MeasListClass has no handle properties
-            end
-            obj.dataTimeSeries = obj2.dataTimeSeries;
-            obj.time = obj2.time;
+            obj.dataTimeSeries  = obj2.dataTimeSeries;
+            obj.time            = obj2.time;
+            obj.timeOffset      = obj2.timeOffset;
         end
         
         
@@ -198,6 +196,12 @@ function Copy(obj, obj2)
             if ~strcmp(obj.name, obj2.name)
                 return;
             end
+            if ~all(size(obj.dataTimeSeries)==size(obj2.dataTimeSeries))
+                return;
+            end
+            if ~all(size(obj.time)==size(obj2.time))
+                return;
+            end
             if ~all(obj.dataTimeSeries(:)==obj2.dataTimeSeries(:))
                 return;
             end
diff --git a/DataTree/AcquiredData/Snirf/DataClass.m b/DataTree/AcquiredData/Snirf/DataClass.m
index fa5a3344..1e39d8e3 100644
--- a/DataTree/AcquiredData/Snirf/DataClass.m
+++ b/DataTree/AcquiredData/Snirf/DataClass.m
@@ -242,7 +242,7 @@
             end
             
             hdf5write_safe(fid, [location, '/dataTimeSeries'], obj.dataTimeSeries, 'array');
-            hdf5write_safe(fid, [location, '/time'], obj.time, 'array');
+            hdf5write_safe(fid, [location, '/time'], obj.time, 'vector');
             
             for ii = 1:length(obj.measurementList)
                 obj.measurementList(ii).SaveHdf5(fid, [location, '/measurementList', num2str(ii)]);
@@ -280,10 +280,19 @@
                     return;
                 end
             end
+            if all(obj.dataTimeSeries==0)
+                return
+            end
             b = true;
         end
         
         
+        % -------------------------------------------------------
+        function b = IsValid(obj)
+            b = obj.IsDataValid();
+        end
+        
+        
         % ----------------------------------------------------------------------
         function err = ErrorCheck(obj, err, params)
             if ~exist('params','var')
@@ -1067,7 +1076,6 @@ function Copy(obj, obj2)
                 obj = DataClass();
             end
             if isempty(obj2)
-                obj = DataClass();
                 return;
             end
             if ~isa(obj2, 'DataClass')
@@ -1086,7 +1094,6 @@ function CopyMeasurementList(obj, obj2)
                 obj = DataClass();
             end
             if isempty(obj2)
-                obj = DataClass();
                 return;
             end
             if ~isa(obj2, 'DataClass')
diff --git a/DataTree/AcquiredData/Snirf/MeasListClass.m b/DataTree/AcquiredData/Snirf/MeasListClass.m
index 2c80f071..c785ee82 100644
--- a/DataTree/AcquiredData/Snirf/MeasListClass.m
+++ b/DataTree/AcquiredData/Snirf/MeasListClass.m
@@ -5,12 +5,17 @@
         sourceIndex
         detectorIndex
         wavelengthIndex
+        wavelengthActual
+        wavelengthEmissionActual
         dataType
+        dataUnit        
         dataTypeLabel
         dataTypeIndex   % Used for condition when dataType=99999 ("Processed") and dataTypeLabel='HRF...'
         sourcePower
         detectorGain
         moduleIndex
+        sourceModuleIndex
+        detectorModuleIndex
     end
     
        
@@ -40,15 +45,20 @@
             
             % Fields which are part of the SNIRF spec which are loaded and saved 
             % from/to SNIRF files
-            obj.sourceIndex      = 0;
-            obj.detectorIndex    = 0;
-            obj.wavelengthIndex  = 0;
-            obj.dataType         = 0;
-            obj.dataTypeLabel    = '';
-            obj.dataTypeIndex    = 0;
-            obj.sourcePower      = 0;
-            obj.detectorGain     = 0;
-            obj.moduleIndex      = 0;
+            obj.sourceIndex                 = 0;
+            obj.detectorIndex               = 0;
+            obj.wavelengthIndex             = 0;
+            obj.wavelengthActual            = 0;
+            obj.wavelengthEmissionActual    = 0;
+            obj.dataType                    = 0;
+            obj.dataUnit                    = '';
+            obj.dataTypeLabel               = '';
+            obj.dataTypeIndex               = 0;
+            obj.sourcePower                 = 0;
+            obj.detectorGain                = 0;
+            obj.moduleIndex                 = 0;
+            obj.sourceModuleIndex           = 0;
+            obj.detectorModuleIndex         = 0;
             
             dataTypeValues = DataTypeValues();
 
@@ -114,16 +124,20 @@
                 [gid, fid] = HDF5_GroupOpen(fileobj, location);
                 
                 % Load datasets
-                obj.sourceIndex     = HDF5_DatasetLoad(gid, 'sourceIndex');
-                obj.detectorIndex   = HDF5_DatasetLoad(gid, 'detectorIndex');
-                obj.wavelengthIndex = HDF5_DatasetLoad(gid, 'wavelengthIndex');
-                obj.dataType        = HDF5_DatasetLoad(gid, 'dataType');
-                obj.dataTypeIndex   = HDF5_DatasetLoad(gid, 'dataTypeIndex');
-                obj.dataTypeLabel   = HDF5_DatasetLoad(gid, 'dataTypeLabel', obj.dataTypeLabel);
-                obj.detectorIndex   = HDF5_DatasetLoad(gid, 'detectorIndex');
-                obj.sourcePower     = HDF5_DatasetLoad(gid, 'sourcePower');
-                obj.sourcePower     = HDF5_DatasetLoad(gid, 'sourcePower');
-                obj.moduleIndex     = HDF5_DatasetLoad(gid, 'moduleIndex');
+                obj.sourceIndex                 = HDF5_DatasetLoad(gid, 'sourceIndex');
+                obj.detectorIndex               = HDF5_DatasetLoad(gid, 'detectorIndex');
+                obj.wavelengthIndex             = HDF5_DatasetLoad(gid, 'wavelengthIndex');
+                obj.wavelengthActual            = HDF5_DatasetLoad(gid, 'wavelengthActual');
+                obj.wavelengthEmissionActual    = HDF5_DatasetLoad(gid, 'wavelengthEmissionActual');
+                obj.dataType                    = HDF5_DatasetLoad(gid, 'dataType');
+                obj.dataUnit                    = HDF5_DatasetLoad(gid, 'dataUnit', obj.dataUnit);
+                obj.dataTypeIndex               = HDF5_DatasetLoad(gid, 'dataTypeIndex');
+                obj.dataTypeLabel  				= HDF5_DatasetLoad(gid, 'dataTypeLabel', obj.dataTypeLabel);
+                obj.sourcePower     			= HDF5_DatasetLoad(gid, 'sourcePower');
+                obj.moduleIndex     			= HDF5_DatasetLoad(gid, 'moduleIndex');
+                obj.detectorGain                = HDF5_DatasetLoad(gid, 'detectorGain');
+                obj.sourceModuleIndex           = HDF5_DatasetLoad(gid, 'sourceModuleIndex');
+                obj.detectorModuleIndex         = HDF5_DatasetLoad(gid, 'detectorModuleIndex');
                 
                 HDF5_GroupClose(fileobj, gid, fid);
             catch
@@ -168,15 +182,41 @@
                 return;
             end
             
+            % Save required fields
             hdf5write_safe(fid, [location, '/sourceIndex'], uint64(obj.sourceIndex));
             hdf5write_safe(fid, [location, '/detectorIndex'], uint64(obj.detectorIndex));
             hdf5write_safe(fid, [location, '/wavelengthIndex'], uint64(obj.wavelengthIndex));
             hdf5write_safe(fid, [location, '/dataType'], uint64(obj.dataType));
-            hdf5write_safe(fid, [location, '/dataTypeLabel'], obj.dataTypeLabel);
             hdf5write_safe(fid, [location, '/dataTypeIndex'], uint64(obj.dataTypeIndex));
-            hdf5write_safe(fid, [location, '/sourcePower'], obj.sourcePower);
-            hdf5write_safe(fid, [location, '/detectorGain'], obj.detectorGain);
-            hdf5write_safe(fid, [location, '/moduleIndex'], uint64(obj.moduleIndex));
+
+            % Save optional fields
+            if ~isempty(obj.wavelengthActual)
+                hdf5write_safe(fid, [location, '/wavelengthActual'], obj.wavelengthActual);
+            end
+            if ~isempty(obj.wavelengthEmissionActual)
+                hdf5write_safe(fid, [location, '/wavelengthEmissionActual'], obj.wavelengthEmissionActual);
+            end
+            if ~isempty(obj.dataUnit)
+                hdf5write_safe(fid, [location, '/dataUnit'], obj.dataUnit);
+            end
+            if ~isempty(obj.dataTypeLabel)
+                hdf5write_safe(fid, [location, '/dataTypeLabel'], obj.dataTypeLabel);
+            end
+            if ~isempty(obj.sourcePower)
+                hdf5write_safe(fid, [location, '/sourcePower'], obj.sourcePower);
+            end
+            if ~isempty(obj.moduleIndex)
+                hdf5write_safe(fid, [location, '/moduleIndex'], uint64(obj.moduleIndex));
+            end
+            if ~isempty(obj.detectorGain)
+                hdf5write_safe(fid, [location, '/detectorGain'], obj.detectorGain);
+            end
+            if ~isempty(obj.sourceModuleIndex)
+                hdf5write_safe(fid, [location, '/sourceModuleIndex'], uint64(obj.sourceModuleIndex));
+            end
+            if ~isempty(obj.detectorModuleIndex)
+                hdf5write_safe(fid, [location, '/detectorModuleIndex'], uint64(obj.detectorModuleIndex));
+            end
         end
 
                 
diff --git a/DataTree/AcquiredData/Snirf/ProbeClass.m b/DataTree/AcquiredData/Snirf/ProbeClass.m
index 8992acfe..66880366 100644
--- a/DataTree/AcquiredData/Snirf/ProbeClass.m
+++ b/DataTree/AcquiredData/Snirf/ProbeClass.m
@@ -3,21 +3,24 @@
     properties
         wavelengths
         wavelengthsEmission
+        sourceLabels
         sourcePos2D
+        sourcePos3D
+        detectorLabels
         detectorPos2D
+        detectorPos3D
+        landmarkLabels
         landmarkPos2D
         landmarkPos3D
-        sourcePos3D
-        detectorPos3D
         frequencies
         timeDelays
         timeDelayWidths
         momentOrders
         correlationTimeDelays
         correlationTimeDelayWidths
-        sourceLabels
-        detectorLabels
-        landmarkLabels
+        coordinateSystem
+        coordinateSystemDescription
+        useLocalIndex
     end
     
     properties (Access = private)
@@ -33,6 +36,24 @@
             % Set class properties not part of the SNIRF format
             obj.SetFileFormat('hdf5');
             obj.scaling = 1;            
+
+            obj.wavelengths          = [];
+            obj.wavelengthsEmission  = [];
+            obj.sourcePos2D  = [];
+            obj.detectorPos2D  = [];
+            obj.sourcePos3D  = [];
+            obj.detectorPos3D  = [];
+            obj.frequencies  = [];
+            obj.timeDelays  = [];
+            obj.timeDelayWidths  = [];
+            obj.momentOrders = [];
+            obj.correlationTimeDelays = [];
+            obj.correlationTimeDelayWidths = [];
+            obj.sourceLabels = {};
+            obj.detectorLabels = {};
+            obj.coordinateSystem = '';
+            obj.coordinateSystemDescription = '';
+            obj.useLocalIndex = 0;
             
             % Set SNIRF fomat properties
             if nargin>0
@@ -71,37 +92,16 @@
                         obj.landmarkPos3D = SD.DummyPos;
                         obj.landmarkLabels = SD.Landmarks3D.labels;
                     end
-                    obj.frequencies  = 1;
-                    obj.timeDelays  = 0;
-                    obj.timeDelayWidths  = 0;
-                    obj.momentOrders = [];
-                    obj.correlationTimeDelays = 0;
-                    obj.correlationTimeDelayWidths = 0;
-                    for ii=1:size(SD.SrcPos)
+                    for ii = 1:size(SD.SrcPos)
                         obj.sourceLabels{ii} = ['S',num2str(ii)];
                     end
-                    for ii=1:size(SD.DetPos)
+                    for ii = 1:size(SD.DetPos)
                         obj.detectorLabels{ii} = ['D',num2str(ii)];
                     end
                 elseif ischar(varargin{1})
                     obj.SetFilename(varargin{1});
                     obj.Load(varargin{1});
                 end
-            else
-                obj.wavelengths          = [];
-                obj.wavelengthsEmission  = [];
-                obj.sourcePos2D  = [];
-                obj.detectorPos2D  = [];
-                obj.sourcePos3D  = [];
-                obj.detectorPos3D  = [];
-                obj.frequencies  = 1;
-                obj.timeDelays  = 0;
-                obj.timeDelayWidths  = 0;
-                obj.momentOrders = [];
-                obj.correlationTimeDelays = 0;
-                obj.correlationTimeDelayWidths = 0;
-                obj.sourceLabels = {};
-                obj.detectorLabels = {};
             end
         end
 
@@ -150,8 +150,40 @@ function BackwardCompatibility(obj)
 
         
         
+        % -------------------------------------------------------
+        function ScaleProjection2D(obj)
+            if size(obj.sourcePos2D,1)
+                dm_src2d = distmatrix([obj.sourcePos2D, zeros(size(obj.sourcePos2D,1),1)]);
+            else
+                dm_src2d = distmatrix(obj.sourcePos2D, zeros);
+            end
+            dm_src3d = distmatrix(obj.sourcePos3D);
+            exp = round(log10(mean(dm_src3d(dm_src3d(:)>0)))) - round(log10(mean(dm_src2d(dm_src2d(:)>0))));
+            scaleFactor = 10^exp;
+            obj.sourcePos2D = obj.sourcePos2D*scaleFactor;
+
+            if size(obj.detectorPos2D,1)
+                dm_det2d = distmatrix([obj.detectorPos2D, zeros(size(obj.detectorPos2D,1),1)]);
+            else
+                dm_det2d = distmatrix(obj.detectorPos2D, zeros);
+            end
+            dm_det3d = distmatrix(obj.detectorPos3D);
+            exp = round(log10(mean(dm_det3d(dm_det3d(:)>0)))) - round(log10(mean(dm_det2d(dm_det2d(:)>0))));
+            scaleFactor = 10^exp;
+            obj.detectorPos2D = obj.detectorPos2D*scaleFactor;
+        end
+            
+        
+        
+        
         % -------------------------------------------------------
         function Project_3D_to_2D(obj)             
+            if all(obj.sourcePos2D(:)==0)
+                obj.sourcePos2D = [];
+            end
+            if all(obj.detectorPos2D(:)==0)
+                obj.detectorPos2D = [];
+            end
             if isempty(obj.sourcePos2D) && isempty(obj.detectorPos2D)
                 if isempty(obj.landmarkPos3D) || ~obj.isValidLandmarkLabels()
                     nSource = size(obj.sourcePos3D,1);
@@ -220,6 +252,7 @@ function Project_3D_to_2D(obj)
                         %
                         obj.sourcePos2D = convert_optodepos_to_circlular_2D_pos(obj.sourcePos3D, T, norm_factor);
                         obj.detectorPos2D = convert_optodepos_to_circlular_2D_pos(obj.detectorPos3D, T, norm_factor);
+                        obj.ScaleProjection2D();
                     end
                 end
             end
@@ -356,23 +389,44 @@ function Project_3D_to_2D(obj)
             obj.landmarkPos3D  = obj.landmarkPos3D / obj.scaling;
             
             % Now save
-            hdf5write_safe(fid, [location, '/wavelengths'], obj.wavelengths, 'array');
-            hdf5write_safe(fid, [location, '/wavelengthsEmission'], obj.wavelengthsEmission, 'array');
+            hdf5write_safe(fid, [location, '/wavelengths'], obj.wavelengths, 'vector');
             hdf5write_safe(fid, [location, '/sourcePos2D'], obj.sourcePos2D, 'array');
             hdf5write_safe(fid, [location, '/detectorPos2D'], obj.detectorPos2D, 'array');
             hdf5write_safe(fid, [location, '/landmarkPos2D'], obj.landmarkPos2D, 'array');
             hdf5write_safe(fid, [location, '/sourcePos3D'], obj.sourcePos3D, 'array');
             hdf5write_safe(fid, [location, '/detectorPos3D'], obj.detectorPos3D, 'array');
             hdf5write_safe(fid, [location, '/landmarkPos3D'], obj.landmarkPos3D, 'array');
-            hdf5write_safe(fid, [location, '/frequencies'], obj.frequencies, 'array');
-            hdf5write_safe(fid, [location, '/timeDelays'], obj.timeDelays, 'array');
-            hdf5write_safe(fid, [location, '/timeDelayWidths'], obj.timeDelayWidths, 'array');
-            hdf5write_safe(fid, [location, '/momentOrders'], obj.momentOrders, 'array');
-            hdf5write_safe(fid, [location, '/correlationTimeDelays'], obj.correlationTimeDelays, 'array');
-            hdf5write_safe(fid, [location, '/correlationTimeDelayWidths'], obj.correlationTimeDelayWidths, 'array');
-            hdf5write_safe(fid, [location, '/sourceLabels'], obj.sourceLabels, 'array');
-            hdf5write_safe(fid, [location, '/detectorLabels'], obj.detectorLabels, 'array');
-            hdf5write_safe(fid, [location, '/landmarkLabels'], obj.landmarkLabels, 'array');
+            
+            if ~isempty(obj.sourceLabels)
+                hdf5write_safe(fid, [location, '/sourceLabels'], obj.sourceLabels, 'array');
+            end
+            if ~isempty(obj.detectorLabels)
+                hdf5write_safe(fid, [location, '/detectorLabels'], obj.detectorLabels, 'array');
+            end
+            if ~isempty(obj.landmarkLabels)
+                hdf5write_safe(fid, [location, '/landmarkLabels'], obj.landmarkLabels, 'array');
+            end
+            if ~isempty(obj.wavelengthsEmission)
+                hdf5write_safe(fid, [location, '/wavelengthsEmission'], obj.wavelengthsEmission, 'vector');
+            end
+            if ~isempty(obj.frequencies)
+                hdf5write_safe(fid, [location, '/frequencies'], obj.frequencies, 'vector');
+            end
+            if ~isempty(obj.timeDelays)
+                hdf5write_safe(fid, [location, '/timeDelays'], obj.timeDelays, 'vector');
+            end
+            if ~isempty(obj.timeDelayWidths)
+                hdf5write_safe(fid, [location, '/timeDelayWidths'], obj.timeDelayWidths, 'vector');
+            end
+            if ~isempty(obj.momentOrders)
+                hdf5write_safe(fid, [location, '/momentOrders'], obj.momentOrders, 'vector');
+            end
+            if ~isempty(obj.correlationTimeDelays)
+                hdf5write_safe(fid, [location, '/correlationTimeDelays'], obj.correlationTimeDelays, 'vector');
+            end
+            if ~isempty(obj.correlationTimeDelayWidths)
+                hdf5write_safe(fid, [location, '/correlationTimeDelayWidths'], obj.correlationTimeDelayWidths, 'vector');
+            end
         end
         
         
@@ -541,6 +595,16 @@ function Project_3D_to_2D(obj)
                     return;
                 end
             end
+            if ~isempty(obj.sourcePos2D)
+                if size(obj.sourcePos2D,2) < 2
+                    return;
+                end
+            end
+            if ~isempty(obj.sourcePos3D)
+                if size(obj.sourcePos3D,2) < 3
+                    return;
+                end
+            end
             if iscolumn(obj.detectorPos2D)
                 return;
             end
@@ -549,6 +613,16 @@ function Project_3D_to_2D(obj)
                     return;
                 end
             end
+            if ~isempty(obj.detectorPos2D)
+                if size(obj.detectorPos2D,2) < 2
+                    return;
+                end
+            end
+            if ~isempty(obj.detectorPos3D)
+                if size(obj.detectorPos3D,2) < 3
+                    return;
+                end
+            end
             b = true;
         end
         
diff --git a/DataTree/AcquiredData/Snirf/SnirfClass.m b/DataTree/AcquiredData/Snirf/SnirfClass.m
index 679d730d..1f53ced7 100644
--- a/DataTree/AcquiredData/Snirf/SnirfClass.m
+++ b/DataTree/AcquiredData/Snirf/SnirfClass.m
@@ -1178,7 +1178,11 @@ function SetMetaDataTags(obj, val)
                     p = length(obj2.GetTime());
                     q = length(obj.aux(ii).GetTime());
                 end
-                datamat(:,ii) = resample(obj.aux(ii).GetDataTimeSeries(), p, q);
+                m = resamplee(obj.aux(ii).GetDataTimeSeries(), p, q);
+                if isrow(m)
+                    m = m';
+                end
+                datamat = [datamat, m]; %#ok<AGROW>
             end
         end
         
@@ -1254,6 +1258,26 @@ function SetStims_MatInput(obj, s, t)
         end
         
         
+        
+        % ---------------------------------------------------------
+        function s = EditStates_MatInput(obj, s, t)
+            if nargin<2
+                return
+            end
+            if isempty(t)
+                return;
+            end
+            errmargin = mean(diff(t));
+            for ii = 1:size(s,2)
+                tidxs = find(s(:,ii)~=0);
+                for jj = 1:length(tidxs)
+                    obj.stim(ii).EditState(t(tidxs(jj)), s(tidxs(jj),ii), errmargin);
+                end
+            end
+        end
+        
+
+        
         % ---------------------------------------------------------
         function s = GetStims(obj, t)
             % Returns a .nirs style stim signal. Stim state marks are
diff --git a/DataTree/AcquiredData/Snirf/StimClass.m b/DataTree/AcquiredData/Snirf/StimClass.m
index 0af23906..72fae473 100644
--- a/DataTree/AcquiredData/Snirf/StimClass.m
+++ b/DataTree/AcquiredData/Snirf/StimClass.m
@@ -228,7 +228,9 @@ function Copy(obj, obj2)
             end
             
             % Now check contents
-            if ~all( abs(obj.data(:)-obj2.data(:)) < (obj.errmargin/10) )
+            d1 = sortrows(obj.data);
+            d2 = sortrows(obj2.data);
+            if ~all( abs(d1(:)-d2(:)) < (obj.errmargin/10) )
                 return;
             end
             
@@ -487,14 +489,18 @@ function EditAmplitude(obj, tPts, amp)
 
         
         % ----------------------------------------------------------------------------------
-        function EditState(obj, tPts, state)
+        function EditState(obj, tPts, state, errmargin)
+            if ~exist('state','var')
+                state = 1;
+            end
             if isempty(obj.data)
                 return;
             end
-            if ~exist('state','var')
-                state = 1;
+            if exist('errmargin','var')
+                k = obj.GetIncludedStimIdxs(tPts, errmargin);
+            else
+                k = obj.GetIncludedStimIdxs(tPts);
             end
-            k = GetIncludedStimIdxs(obj, tPts);
             obj.states(k,2) = state;
         end
         
@@ -529,12 +535,14 @@ function SetTpts(obj, tPts, k)
         
         
         % -------------------------------------------------------
-        function k = GetIncludedStimIdxs(obj, tPts)
+        function k = GetIncludedStimIdxs(obj, tPts, errmargin)
             k = [];
-            if length(tPts)<2
-                errmargin = obj.errmargin;
-            else
-                errmargin = mean(diff(tPts));
+            if ~exist('errmargin','var')
+                if length(tPts)<2
+                    errmargin = obj.errmargin;
+                else
+                    errmargin = mean(diff(tPts));
+                end
             end
             for ii = 1:size(obj.data,1)
                 if ~isempty(find( abs(obj.data(ii,1)-tPts) < errmargin))
@@ -555,7 +563,7 @@ function SetDuration(obj, duration, tPts)
             if ~exist('tPts','var')
                 tPts = obj.data(:,1);
             end
-            k = GetIncludedStimIdxs(obj, tPts);
+            k = obj.GetIncludedStimIdxs(tPts);
             obj.data(k,2) = duration;
         end
 
diff --git a/DataTree/AppSettings.cfg b/DataTree/AppSettings.cfg
index f7c804fc..6c0cdd72 100644
--- a/DataTree/AppSettings.cfg
+++ b/DataTree/AppSettings.cfg
@@ -38,7 +38,7 @@ No
 No
 
 % Quiet Mode # On, Off
-Off
+On
 
 % Replace TSV File Tabs with Spaces # auto-replace, ask me, don't replace and don't ask me
 ask me
diff --git a/DataTree/DataTreeClass.m b/DataTree/DataTreeClass.m
index c26d464a..98dbc1a0 100644
--- a/DataTree/DataTreeClass.m
+++ b/DataTree/DataTreeClass.m
@@ -788,7 +788,6 @@ function Save(obj, hwait)
             
             t1 = tic;
             for ii = 1:length(obj.groups)
-                obj.logger.Write('Saving group %d in %s\n', ii, [obj.groups(ii).pathOutputAlt, obj.groups(ii).GetFilename()]);
                 obj.groups(ii).Save(hwait);
             end
             obj.logger.Write('Completed saving processing results for all groups in %0.3f seconds.\n', toc(t1));
diff --git a/DataTree/Examples/EditStimExample.m b/DataTree/Examples/EditStimExample.m
index 8fbdefb8..c17c6ed8 100644
--- a/DataTree/Examples/EditStimExample.m
+++ b/DataTree/Examples/EditStimExample.m
@@ -24,7 +24,7 @@ function EditStimExample(dataSetDir)
 %
 %   Examples:
 %       cd <local_path>/DataTree
-%       EditStimExample('./Examples/SubjDataSample')
+%       EditStimExample('./Examples/SampleData')
 %
 
 
@@ -33,7 +33,7 @@ function EditStimExample(dataSetDir)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if ~exist('dataSetDir','var')
     f = which('EditStimExample');
-    dataSetDir = [fileparts(f), '/SubjDataSample']; 
+    dataSetDir = [fileparts(f), '/SampleData']; 
 end
 
 
diff --git a/DataTree/GroupClass.m b/DataTree/GroupClass.m
index 26cc682d..42061a5f 100644
--- a/DataTree/GroupClass.m
+++ b/DataTree/GroupClass.m
@@ -308,9 +308,7 @@ function Add(obj, subj, sess, run)
                 if ~exist(fnameFull, 'file')
                     procStreamDefault.SaveConfigFile(fnameFull, o.type);
                 end
-                
-                obj.logger.Write('Attempting to load %s-level proc stream from %s\n', o.type, fnameFull);
-                
+                               
                 % Load file to the first empty procStream in the dataTree at each processing level
                 err = o.LoadProcStreamConfigFile(fnameFull);
 
@@ -329,7 +327,7 @@ function Add(obj, subj, sess, run)
                 % disseminate it to all subjects and all runs in this group
                 else
                     
-                    obj.logger.Write('Loading proc stream from %s\n', fnameFull);
+                    obj.logger.Write('Loading %s-level proc stream from %s\n', o.type, fnameFull);
                     g.CopyFcalls(o.procStream, o.type);
                     
                 end
@@ -694,7 +692,7 @@ function Save(obj, hwait)
                 return
             end
             
-            obj.logger.Write('Saving processed data in %s\n', [obj.path, obj.outputDirname, obj.outputFilename]);
+            obj.logger.Write('Saving group tree processing stream input in  %s\n', [obj.path, obj.outputDirname, obj.outputFilename]);
             
             if ishandle(hwait)
                 obj.logger.Write('Auto-saving processing results ...\n', obj.logger.ProgressBar(), hwait);
diff --git a/DataTree/ProcStream/ProcResultClass.m b/DataTree/ProcStream/ProcResultClass.m
index 375430ce..4a43d329 100644
--- a/DataTree/ProcStream/ProcResultClass.m
+++ b/DataTree/ProcStream/ProcResultClass.m
@@ -283,7 +283,7 @@ function SaveInit(obj, pathname, filename)
             end
             output = ProcResultClass();
             
-            obj.logger.Write(sprintf('Initializing empty derived data output file:   %s', [filesepStandard(pathname), filename, '.mat']));
+            obj.logger.WriteFileOnly(sprintf('Initializing empty derived data output file:   %s', [filesepStandard(pathname), filename, '.mat']), 661, 529, 23); 
             save([filesepStandard(pathname), filename, '.mat'], '-mat', 'output');
         end
         
diff --git a/DataTree/ProcStream/ProcStreamClass.m b/DataTree/ProcStream/ProcStreamClass.m
index c5d7bc57..b570ee60 100644
--- a/DataTree/ProcStream/ProcStreamClass.m
+++ b/DataTree/ProcStream/ProcStreamClass.m
@@ -410,7 +410,7 @@ function ExportProcStream(obj, filename, fcalls)
             [p,f] = fileparts(temp); 
             fname = [filesepStandard(p), f, '_processing.json'];
             if strcmpi(cfg.GetValue('Export Processing Stream Functions'), 'yes')
-                logger.Write('Saving processing stream  %s:\n', fname);
+                logger.Write('Saving processing stream function call chain in %s:\n', fname);
                 appname = sprintf('%s', getNamespace());
                 vernum  = sprintf('v%s', getVernum(appname));
                 dt      = sprintf('%s', char(datetime(datetime, 'Format','MMMM d, yyyy,   HH:mm:ss')));
diff --git a/DataTree/TreeNodeClass.m b/DataTree/TreeNodeClass.m
index 9794abdc..25c4978b 100644
--- a/DataTree/TreeNodeClass.m
+++ b/DataTree/TreeNodeClass.m
@@ -2186,6 +2186,17 @@ function ApplyParamEditToAll(obj, iFcall, iParam, val)
 
         
         
+        % ------------------------------------------------------------
+        function msg = GetErrorMsg(obj)
+            msg = '';
+            if isempty(obj.acquired)
+                return;
+            end
+            msg = obj.acquired.GetErrorMsg();
+        end
+
+        
+        
         % ------------------------------------------------------------
         function Print(obj, indent)
             obj.logger.Write('%s%s\n', blanks(indent), obj.procStream.output.SetFilename([obj.path, obj.GetOutputFilename()]) );
diff --git a/DataTree/Version.txt b/DataTree/Version.txt
index 4761f0e9..87967a79 100644
--- a/DataTree/Version.txt
+++ b/DataTree/Version.txt
@@ -1 +1 @@
-1.15.2
\ No newline at end of file
+1.21.1
\ No newline at end of file
diff --git a/Example pipelines/README.txt b/Example pipelines/README.txt
new file mode 100644
index 00000000..93397127
--- /dev/null
+++ b/Example pipelines/README.txt	
@@ -0,0 +1,7 @@
+
+To reproduce the results in the livescripts, simple_pipeline.mlx and glm_pipeline.mlx, select the file in "Brite Fingertapping 5-19-2020.snirf"
+
+	<Homer3 root folder> / SampleData / 	
+	
+That is the best data file for illustrating motion correction. 
+
diff --git a/Example pipelines/glm_pipeline.m b/Example pipelines/glm_pipeline.m
new file mode 100644
index 00000000..3f303140
--- /dev/null
+++ b/Example pipelines/glm_pipeline.m	
@@ -0,0 +1,153 @@
+%% GLM Analysis in Homer3 Toolbox
+% This livescript will demonstrate how to perform GLM analysis using Homer3 
+% toolbox. The following analysis steps will be performed:
+%% 
+%% Loading SNIRF Data Files
+% First, the SNIRF data will be loaded from <https://fnirs.org/resources/software/snirf/ 
+% Shared Near Infrared Spectroscopy Format> (SNIRF), which was designed to facilitate 
+% sharing and analysing NIRS data. A detailed description of the format can be 
+% found <https://github.com/fNIRS/snirf/blob/master/snirf_specification.md#snirf-data-format-summary 
+% here>. We will first select a measurement file. 
+
+clear all; close all;
+
+% Select data file to process
+p = fileparts(which('simple_pipeline.m'));
+while 1
+    [fname, pname] = uigetfile([p, '/../SampleData/', '*.snirf'], 'Select snirf file(s)', 'MultiSelect', 'on');
+    if fname == 0
+        q = menu('You have not selected a data file. Click OK to continue or CANCEL to quit', ...
+            {'OK','CANCEL'});
+        if q==2
+            return
+        end
+    else
+        break;
+    end
+end
+
+% Load raw intensity signals and metadata into data_nirs 
+if iscell(fname) 
+    data_nirs = cell(1, numel(fname));
+    for fi = 1:numel(fname)
+        data_nirs{fi} = SnirfClass([pname, fname{fi}]); 
+    end
+else
+    data_nirs = {SnirfClass([pname, fname])}; 
+end
+
+
+%
+% We will now proceed to work with just one measurement data, which we will 
+% call _acquired. _
+acquired = data_nirs{1};
+
+% The _acquired _Snirf data object contains both raw intensity signals for each 
+% wavelength and metadata, such as stimuli information, probes source/detector 
+% geometry, measurement date and time, etc. We will now work on the intensity 
+% signals and, for that, we save a copy of it in the _intensity_nirs_ variable. 
+% It's important to note that just assigning "acquired.data" data structure to 
+% _intensity_nirs _ variable will only point to the data in _acquired _data object. 
+% This means that any further modifications that we do on _intensity_nirs_ will 
+% also be done on the data contained in _acquired. _To prevent this from happening, 
+% we create a copy of it as follows:
+intensity_nirs = acquired.data.copy; % raw intensity signals 
+%% Convert Intensities to Optical Densities
+% We will proceed to convert these signals into optical density (OD) signals. 
+% These signals represent the attenuation of the transmitted light $(I_0)$in tissue, 
+% as measured by the receivers $(I$), for each wavelength. 
+% 
+% $$OD(t, \lambda)$$ $$=-log_{10}$$$$\left(\frac{I(t,\lambda}{I_0(t,\lambda)}\right)}$$
+% 
+% The computation in Homer3 is as follows:
+
+dod = hmrR_Intensity2OD(intensity_nirs);
+%% Bandpass Filtering
+% In a common analysis pipeline, fNIRS data is filtered such that the frequency 
+% of the hemodynamic response remains and frequencies of other signal components 
+% are filtered out. Although there are different opinions on the best filter cutoffs 
+% for estimating the hemodynamic response, most researchers filter out frequencies 
+% above 0.1 Hz (<https://www.frontiersin.org/articles/10.3389/fnhum.2018.00505/full 
+% Pinti et al., 2018>) in the healthy adult population. We will additionally filter 
+% frequencies below 0.01 Hz, and we will do so by means of a bandpass filter.
+
+dodFilt = hmrR_BandpassFilt(dod, 0.00, 0.5);
+% Optical densities to concentrations
+% Once our data has been filtered, we compute the concentration changes in O2Hb 
+% and HHb by means of the Modified Lambert-Beer Law <https://iopscience.iop.org/article/10.1088/0031-9155/33/12/008 
+% (Delpy et al., 1988)>.  This conversion requires the user to insert ppf (partial 
+% pathlength factor) values corresponding to each wavelength. These are age and 
+% wavelength specific. Here we use 6.06.
+
+dhbFilt = hmrR_OD2Conc(dodFilt, acquired.probe, [1, 1]);
+%% *GLM analysis*
+% After performing basic data preprocessing steps, we will proceed with the 
+% GLM analysis. GLM stands for General Linear Model. In GLM analysis, a reference 
+% function, which consists of a linear combination of regressors, is fitted to 
+% the time course of each channel. The design matrix for the GLM reference function 
+% is created by convolving the boxcar function (boxcar time course acquires values 
+% of 1 when the condition is _on _and 0 when the condition is_ off_) with the 
+% hemodynamic response function. 
+% 
+% We will use _hmrR_GLM _function to run the GLM. Below, we set the parameters 
+% for the function. Please refer to Homer3 function documentation to see an explixit 
+% explanation of these parameters 
+
+trange = [-2 15];
+glmSolveMethod = 1; % Ordinary least squares for inverting design matrix
+idxBasis = 2; % Modified gamma function for modeling HRF
+paramsBasis = [0.1 3.0 1.8 3.0]; % Parameters to model the shape of HRF (recommended values)
+%%
+% In addition, various other regressors can be added to the model. For example, 
+% Homer3 allows regressing Short Separation Channels (SSC), model drifts and auxiliary 
+% signals. To enable these regressors, we will specify some more parameters. Specifically, 
+% the drift order will be set to 3. Since the SSC were included in the set-up 
+% of an example measurement used here, we will set the SSC threshold to 1 cm distance.
+
+driftOrder = 3; % How many drift regressors to use (polynomial order)
+rhoSD_ssThresh = 10; % Threshold for SSC channels, in cm
+flagNuisanceRMethod = 0; % Regression is performed with the nearest SS channel
+
+% The specified parameters can now be used to run the GLM.
+
+stim = acquired.stim;
+probe = acquired.probe;
+[data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks, hmrstats] = ...
+    hmrR_GLM(dhbFilt, stim, probe, [], [], [], [], trange, glmSolveMethod, idxBasis, paramsBasis,...
+    rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, []);
+%% Plot Final Results in a Graph
+% Displaying NIRS-data is typically done by plotting both the average changes 
+% of oxy- and deoxyhemoglobin in time. Here, we plot the resulting average * on 
+% HbO and HbR signals
+% Extract plot data from output of last function 
+
+% These are plot figure positioning parameters
+px      = .01;
+py      = .30;
+sx      = .30;
+sy      = .40;
+xstep   = .15;
+
+% Now plot the GLM processed block average 
+tAvg = data_yavg.time;
+yAvg = data_yavg.GetDataTimeSeries('reshape');
+nChannels = size(yAvg,3);
+channels = 1:nChannels;
+hbPlotNames = {'HbO ', 'HbR ', 'HbT '};
+hbType   = [1, 2, 3];
+for iHb = 1:length(hbType)
+    figure('menubar','none', 'numbertitle','off', 'name',[hbPlotNames{iHb}, 'Concentration Block Average'], ...
+        'units','normalized','position',[px, py, sx, sy]);
+    px = px + xstep;
+    for iCh = channels
+        plot(tAvg, yAvg(:, hbType(iHb), channels(iCh)));
+        hold on;
+    end
+    xlim([min(tAvg), max(tAvg)]);
+    hold off;
+    pause(1)
+end
+xlim([min(tAvg), max(tAvg)]);
+hold off;
+
+
diff --git a/Example pipelines/glm_pipeline.mlx b/Example pipelines/glm_pipeline.mlx
new file mode 100644
index 00000000..0e348896
Binary files /dev/null and b/Example pipelines/glm_pipeline.mlx differ
diff --git a/Example pipelines/simple_pipeline.m b/Example pipelines/simple_pipeline.m
new file mode 100644
index 00000000..5b219f79
--- /dev/null
+++ b/Example pipelines/simple_pipeline.m	
@@ -0,0 +1,259 @@
+%% 
+% This livescript will introduce basic principles of Homer3 and guide you through 
+% the implementation of a simple pipeline without the Homer3 GUI, in a tutorial 
+% style. You will learn how to:
+%% Loading SNIRF Data Files
+% We will load the files from the <https://fnirs.org/resources/software/snirf/ 
+% Shared Near Infrared Spectroscopy Format> (SNIRF), which was designed to facilitate 
+% sharing and analysing NIRS data. A detailed description about the format can 
+% be found <https://github.com/fNIRS/snirf/blob/master/snirf_specification.md#snirf-data-format-summary 
+% here>. We can select one or more measurement files in the dialog box. Here, 
+% we chose to load all three measurement files in our experiment into _data_nirs._
+
+clear all; close all;
+
+% Select data file to process
+p = fileparts(which('simple_pipeline.m'));
+while 1
+    [fname, pname] = uigetfile([p, '/../SampleData/', '*.snirf'], 'Select snirf file(s)', 'MultiSelect', 'on');
+    if fname == 0
+        q = menu('You have not selected a data file. Click OK to continue or CANCEL to quit', ...
+            {'OK','CANCEL'});
+        if q==2
+            return
+        end
+    else
+        break;
+    end
+end
+
+% Load raw intensity signals and metadata into data_nirs 
+if iscell(fname) 
+    data_nirs = cell(1, numel(fname));
+    for fi = 1:numel(fname)
+        data_nirs{fi} = SnirfClass([pname, fname{fi}]); 
+    end
+else
+    data_nirs = {SnirfClass([pname, fname])}; 
+end
+
+
+%
+% We will now proceed to work with just one measurement data, which we will 
+% call _acquired. _
+acquired = data_nirs{1};
+
+% The _acquired _Snirf data object contains both raw intensity signals for each 
+% wavelength and metadata, such as stimuli information, probes source/detector 
+% geometry, measurement date and time, etc. We will now work on the intensity 
+% signals and, for that, we save a copy of it in the _intensity_nirs_ variable. 
+% It's important to note that just assigning "acquired.data" data structure to 
+% _intensity_nirs _ variable will only point to the data in _acquired _data object. 
+% This means that any further modifications that we do on _intensity_nirs_ will 
+% also be done on the data contained in _acquired. _To prevent this from happening, 
+% we create a copy of it as follows:
+intensity_nirs = acquired.data.copy; % raw intensity signals
+%% Convert Intensities to Optical Densities
+% We will proceed to convert these signals into optical density (OD) signals. 
+% These signals represent the attenuation of the transmitted light $(I_0)$in tissue, 
+% as measured by the receivers $(I$), for each wavelength. 
+% 
+% $$OD(t, \lambda)$$ $$=-log_{10}$$$$\left(\frac{I(t,\lambda}{I_0(t,\lambda)}\right)}$$
+% 
+% The computation in Homer3 is as follows:
+
+dod = hmrR_Intensity2OD(intensity_nirs);
+%% Detect Motion Artifacts
+% Parts of the data can be contaminated by artefacts, which could be systemic 
+% artefacts (e.g. heartbeat,respiration), or be caused by extrinsic factors (e.g. 
+% subject movement) (<https://www.spiedigitallibrary.org/journals/neurophotonics/volume-3/issue-3/031405/False-positives-and-false-negatives-in-functional-near-infrared-spectroscopy/10.1117/1.NPh.3.3.031405.full?SSO=1 
+% Tachtsidis&Scholkmann, 2016>). Both these types of artefacts need to be handled 
+% separately. In this tutorial, we focus on the latter.
+% 
+% Motion artefacts detection is often done by visual inspection of non-physiological 
+% data patterns or automated analysis. Homer3 provides the hmrR_MotionArtifactByChannel 
+% function, which detects outliers using standard deviation and amplitude thresholds. 
+% 
+% We will set these thresholds to be 50 and 5, respectively, and choose a 1 
+% second sliding window in which the channel-wise amplitude and standard deviation 
+% will be evaluated to detect outliers. If any active data channel exhibits a 
+% signal change greater than STDEVthresh or AMPthresh, then a segment of data 
+% around that time point is marked as a motion artifact. The channel-wise motion 
+% artifacts are recorded in the output matrix _tIncAutoCh, _where 1s represent 
+% motion-free data and 0s represent motion contaminated data. Time points in the 
+% recording where at least one channel has a motion artifact are indicated in 
+% a similar manner by the vector _tIncAuto_.
+
+% PARAMETERS:
+tMotion = 1.0;
+tMask = 1.0;
+STDEVthresh = 10;
+AMPthresh = 1.00;
+[tIncAuto, tIncAutoCh] = hmrR_MotionArtifactByChannel(dod, acquired.probe, [], [], [], ...
+    tMotion, tMask, STDEVthresh, AMPthresh);
+%% Correct Motion Artifacts
+% Having already detected the motion artifacts in our data, we will now proceed 
+% to artefact handling. Artefact handling can be done in two ways: either, the 
+% data is removed from further analysis, or it is interpolated e.g. using cubic 
+% spline interpolation algorithm (<https://pubmed.ncbi.nlm.nih.gov/20308772/ Scholkmann 
+% et al., 2010>). We opted for the latter, which is imlpemented in homer as follows:
+
+dodMC = hmrR_MotionCorrectSpline(dod.copy, [], tIncAutoCh, 0.99);       % Motion corrected ODs
+
+
+%
+% We now plot the resulting motion-corrected signals.
+t = dod.time;
+y = dod.dataTimeSeries;
+yMC = dodMC.dataTimeSeries;
+
+% Plot figure positioning parameters:   
+px      = .01;
+py      = .30;
+sx      = .30;
+sy      = .40;
+xstep   = .15;
+
+% Channels to plot 
+channels = 1:size(y, 2);
+figure('menubar','none', 'numbertitle','off', 'name','OD Before and After Motion Correction', ...
+    'units','normalized','position',[px, py, sx, sy]);
+px = px + xstep;
+%% 
+% 
+
+ax1 = subplot(211); hold on
+for iCh = 1:length(channels)
+    plot(t, detrend(y(:, iCh)));
+    xlabel('time (s)')
+end
+title('Original raw ODs')
+xlim([t(1), t(end)]);
+
+ax2 = subplot(212); hold on
+for iCh = 1:length(channels)
+    plot(t, detrend(yMC(:, iCh)));
+    xlabel('time (s)')
+end
+title('Motion corrected ODs')
+xlim([t(1), t(end)]);
+
+linkaxes([ax1, ax2])
+pause(1);
+%% Convert Optical Densities to Hemoglobin Concentration Changes
+% Once our data has been corrected for motion, we compute the concentration 
+% changes in O2Hb and HHb by means of the Modified Lambert-Beer Law <https://iopscience.iop.org/article/10.1088/0031-9155/33/12/008 
+% (Delpy et al., 1988)>. This conversion requires the user to insert ppf (partial 
+% pathlength factor) values corresponding to each wavelength. These are age and 
+% wavelength specific. Here we use 6.06.
+
+dhb  = hmrR_OD2Conc(dodMC, acquired.probe, [1,1]);
+
+
+%
+% Data visualization:
+
+% Extract plot data from output of last function 
+t = dhb.time;
+y = dhb.GetDataTimeSeries('reshape');
+
+% fs = 1/diff(dhb.time(1:2)); % sampling frequency
+tstart = acquired.stim(1).data(1,1);
+tend = tstart + 30;
+
+% Hb type - {1 -> HbO,  2 -> HbR,  3 -> HbT}
+hbType   = [1, 2];
+
+% channels to plot 
+subplotDims = [5,4];                            % Subplot dimensions in terms of rows and columns
+nSubplots = subplotDims(1) * subplotDims(2);    % number of subplots = number of subplot rows X number of subplot columns
+nChannels = size(y,3);
+if nSubplots > nChannels
+    channels = 1:nChannels;
+else
+    channels = 1:nSubplots;
+end
+color_map = {'r', 'b', 'green'};
+figure('menubar','none', 'numbertitle','off', 'name','Hb0/HbR Concentration Time Course', ...
+    'units','normalized','position',[px, py, sx*1.3, sy]);
+px = px + xstep;
+
+for iHb = 1:length(hbType)
+    axes_array = [];
+    for iCh = channels
+        axes_array = [axes_array, subplot(5,4,iCh)]; hold on
+        plot(t, detrend(y(:, hbType(iHb), channels(iCh))),...
+            color_map{iHb}, 'LineWidth', 1.5);
+        xlabel('time (s)')
+        xlim([tstart, tend]);
+        hold on
+    end
+    pause(1)
+end
+linkaxes(axes_array, 'x')
+%% Bandpass Filtering
+% In a common analysis pipeline, fNIRS data is filtered such that the frequency 
+% of the hemodynamic response remains and frequencies of other signal components 
+% are filtered out. Although there are different opinions on the best filter cutoffs 
+% for estimating the hemodynamic response, most researchers filter out frequencies 
+% above 0.1 Hz (<https://www.frontiersin.org/articles/10.3389/fnhum.2018.00505/full 
+% Pinti et al., 2018>) in the healthy adult population. We will additionally filter 
+% frequencies bellow 0.01 Hz, and we will do so bu means of a bandpass filter.
+
+dhbFilt = hmrR_BandpassFilt(dhb, 0.01, 0.5);
+%% Average Over Trials
+% Averaging of data is usually first done across trials per condition, and, 
+% if available, over different sessions (also called runs), and then over subjects 
+% or per subject groups. In this tutorial, we have analysed data from a single 
+% subject performing the same task in multiple sessions, and have focused in just 
+% one of these sessions.
+% 
+% We are interested in performing a block average over trials belonging to our 
+% fingertapping task. To do so, we first have a look at the stimulus information 
+% contained out _acquired _Snirf data structure.
+
+for ii = 1:length(acquired.stim)
+    fprintf('Event type:  %s\n', acquired.stim(ii).name)
+end
+
+
+% To perform the block average, the recorded data is segmented into those time 
+% periods of experimental interest. The segmented period includes a few seconds 
+% of time before the task onset and few seconds of time after the task offset. 
+% We select stim(1), corresponding to the fingertapping task, and choose a prestimulus 
+% and poststimulus time of 5 and 30 seconds, respectively. The short period before 
+% the onset of the experimental task is usually used as a baseline. This allows 
+% making inference between a baseline period (e.g. rest) and a task period (e.g. 
+% finger tapping). 
+% 
+% Having set all of these parameters, we then proceed to average over trials 
+% per channel as follows:
+tPre = -5;
+tPost = 30;
+dhbAvg = hmrR_BlockAvg(dhbFilt, acquired.stim(1), [tPre, tPost]);
+
+%% Plot Final Results in a Graph
+% Displaying NIRS-data is typically done by plotting both the average changes 
+% of oxy- and deoxyhemoglobin in time. Here, we plot the resulting average * on 
+% HbO and HbR signals
+% Extract plot data from output of last function 
+tAvg = dhbAvg.time;
+yAvg = dhbAvg.GetDataTimeSeries('reshape');
+channels = 1:nChannels;
+hbPlotNames = {'HbO ', 'HbR ', 'HbT '};
+hbType   = [1, 2, 3];
+for iHb = 1:length(hbType)
+    figure('menubar','none', 'numbertitle','off', 'name',[hbPlotNames{iHb}, 'Concentration Block Average'], ...
+        'units','normalized','position',[px, py, sx, sy]);
+    px = px + xstep;
+    for iCh = channels
+        plot(tAvg, yAvg(:, hbType(iHb), channels(iCh)));
+        hold on;
+    end
+    xlim([min(tAvg), max(tAvg)]);
+    hold off;
+    pause(1)
+end
+xlim([min(tAvg), max(tAvg)]);
+hold off;
+
diff --git a/Example pipelines/simple_pipeline.mlx b/Example pipelines/simple_pipeline.mlx
new file mode 100644
index 00000000..5decc6f9
Binary files /dev/null and b/Example pipelines/simple_pipeline.mlx differ
diff --git a/FuncRegistry/FuncRegClass.m b/FuncRegistry/FuncRegClass.m
index 9ff080ae..50166b0c 100644
--- a/FuncRegistry/FuncRegClass.m
+++ b/FuncRegistry/FuncRegClass.m
@@ -403,15 +403,14 @@ function Copy(obj, obj2)
         % ----------------------------------------------------------------------------------
         function lastdt = DateLastModified(obj)
             try
-                lastdt = datetime(1970,01,01);
-                for ii = 1:length(obj.entries)
-                    for jj = 1:length(obj.userfuncdir)
-                        file = dir([obj.userfuncdir{jj}, '/', obj.entries(ii).GetName(), '.m']);
-                        if ~isempty(file)
-                            if  datetime(file.date,'locale','system') > lastdt
-                                lastdt =  datetime(file.date,'locale','system');
+                lastdt = 0;
+                for jj = 1:length(obj.userfuncdir)
+                    files = dir([obj.userfuncdir{jj}, '*.m']);
+                    for ii = 1:length(files)
+                        if ~isempty(files(ii))
+                            if  files(ii).datenum > lastdt
+                                lastdt =  files(ii).datenum;
                             end
-                            break;
                         end
                     end
                 end
diff --git a/FuncRegistry/RegistriesClass.m b/FuncRegistry/RegistriesClass.m
index 8c87320f..2fccca1a 100644
--- a/FuncRegistry/RegistriesClass.m
+++ b/FuncRegistry/RegistriesClass.m
@@ -385,7 +385,7 @@ function Import(obj, funcpath)
             regfile = dir([obj.userfuncdir{1}, 'Registry.mat']);
             for ii = 1:length(obj.funcReg)
                 try
-                    if obj.funcReg(ii).DateLastModified() > datetime(regfile.date,'locale','system')
+                    if obj.funcReg(ii).DateLastModified() > regfile.datenum
                         return;
                     end
                 catch
diff --git a/FuncRegistry/UserFunctions/Archive/hmrR_BandpassFilt_Nirs.m b/FuncRegistry/UserFunctions/Archive/hmrR_BandpassFilt_Nirs.m
index 81bd5ea0..38a07eb6 100644
--- a/FuncRegistry/UserFunctions/Archive/hmrR_BandpassFilt_Nirs.m
+++ b/FuncRegistry/UserFunctions/Archive/hmrR_BandpassFilt_Nirs.m
@@ -42,26 +42,19 @@
 %[fa,fb]=butter(FilterOrder,lpf*2/fs);
 if FilterType==1 | FilterType==5
     [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low');
-elseif FilterType==4
-%    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low',Filter_Rp,Filter_Rs);
-else
-%    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low',Filter_Rp);
 end
 ylpf=filtfilt(fb,fa,y);
 
+y2 = ylpf;
+
 % high pass filter
+if hpf ~= 0
 FilterType = 1;
 FilterOrder = 5;
 if FilterType==1 | FilterType==5
     [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high');
-elseif FilterType==4
-%    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high',Filter_Rp,Filter_Rs);
-else
-%    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high',Filter_Rp);
 end
-
-if FilterType~=5
-    y2=filtfilt(fb,fa,ylpf); 
-else
-    y2 = ylpf;
+y2=filtfilt(fb,fa,ylpf); 
 end
+
+
diff --git a/FuncRegistry/UserFunctions/Archive/hmrR_MotionCorrectSpline_Nirs.m b/FuncRegistry/UserFunctions/Archive/hmrR_MotionCorrectSpline_Nirs.m
index 22036a06..9ac000dd 100644
--- a/FuncRegistry/UserFunctions/Archive/hmrR_MotionCorrectSpline_Nirs.m
+++ b/FuncRegistry/UserFunctions/Archive/hmrR_MotionCorrectSpline_Nirs.m
@@ -55,6 +55,8 @@
 ml = SD.MeasList;
 mlAct = SD.MeasListAct; % prune bad channels
 
+mlAct = [mlAct; mlAct]; % hardcoded fix to include two wavelengths
+
 lstAct = find(mlAct==1);
 dodSpline = dod;
 t = t(:);  % needs to be a column vector
diff --git a/FuncRegistry/UserFunctions/hmrR_MotionArtifact.m b/FuncRegistry/UserFunctions/hmrR_MotionArtifact.m
index fbe40967..183c9b8f 100644
--- a/FuncRegistry/UserFunctions/hmrR_MotionArtifact.m
+++ b/FuncRegistry/UserFunctions/hmrR_MotionArtifact.m
@@ -99,9 +99,7 @@
     % Get list of active channels 
     mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);
     mlActAuto{iBlk} = mlAct_Initialize(mlActAuto{iBlk}, MeasList);
-    lstAct1 = mlAct_Matrix2IndexList(mlActAuto{iBlk}, MeasList);
-    lstAct2 = mlAct_Matrix2IndexList(mlActMan{iBlk}, MeasList);
-    lstAct = unique([lstAct1(:)', lstAct2(:)']);
+    lstAct = mlAct_CombineIndexLists( mlActMan{iBlk}, mlActAuto{iBlk}, MeasList );
     
     if length(fs)~=1
         fs = 1/(fs(2)-fs(1));
diff --git a/FuncRegistry/UserFunctions/hmrR_MotionArtifactByChannel.m b/FuncRegistry/UserFunctions/hmrR_MotionArtifactByChannel.m
index 3f72fc8f..a848cf3e 100644
--- a/FuncRegistry/UserFunctions/hmrR_MotionArtifactByChannel.m
+++ b/FuncRegistry/UserFunctions/hmrR_MotionArtifactByChannel.m
@@ -108,9 +108,7 @@
     % get list of active channels
     mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);
     mlActAuto{iBlk} = mlAct_Initialize(mlActAuto{iBlk}, MeasList);
-    lstAct1 = mlAct_Matrix2IndexList(mlActAuto{iBlk}, MeasList);
-    lstAct2 = mlAct_Matrix2IndexList(mlActMan{iBlk}, MeasList);
-    lstAct = unique([lstAct1(:)', lstAct2(:)']);
+    lstAct = mlAct_CombineIndexLists( mlActMan{iBlk}, mlActAuto{iBlk}, MeasList );
     
     % set artifact buffer for tMask seconds on each side of spike
     art_buffer = round(tMask*fs); % time in seconds times sample rate    
diff --git a/FuncRegistry/UserFunctions/hmrR_MotionCorrectWavelet.m b/FuncRegistry/UserFunctions/hmrR_MotionCorrectWavelet.m
index 11823a14..f1d44376 100644
--- a/FuncRegistry/UserFunctions/hmrR_MotionCorrectWavelet.m
+++ b/FuncRegistry/UserFunctions/hmrR_MotionCorrectWavelet.m
@@ -70,9 +70,7 @@
     
     mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);
     mlActAuto{iBlk} = mlAct_Initialize(mlActAuto{iBlk}, MeasList);
-    lstAct1 = mlAct_Matrix2IndexList(mlActAuto{iBlk}, MeasList);
-    lstAct2 = mlAct_Matrix2IndexList(mlActMan{iBlk}, MeasList);
-    lstAct = unique([lstAct1(:)', lstAct2(:)']);
+    lstAct = mlAct_CombineIndexLists( mlActMan{iBlk}, mlActAuto{iBlk}, MeasList );
     
     SignalLength = size(dod,1); % #time points of original signal
     N = ceil(log2(SignalLength)); % #of levels for the wavelet decomposition
diff --git a/FuncRegistry/UserFunctions/hmrR_PruneChannels.m b/FuncRegistry/UserFunctions/hmrR_PruneChannels.m
index 4fa31f52..203dcb39 100644
--- a/FuncRegistry/UserFunctions/hmrR_PruneChannels.m
+++ b/FuncRegistry/UserFunctions/hmrR_PruneChannels.m
@@ -30,8 +30,8 @@
 % Prune_Channels: mlActAuto = hmrR_PruneChannels(data, probe, mlActMan, tIncMan, dRange, SNRthresh, SDrange)
 %
 % PARAMETERS:
-% dRange: [1e4, 1e7]
-% SNRthresh: 2
+% dRange: [1.0e4, 1.0e7]
+% SNRthresh: 2.0
 % SDrange: [0.0, 45.0]
 %
 % TO DO:
@@ -127,6 +127,22 @@
     
     chanList = chanList(:) & mlActMan{iBlk}(:,3);
     
+    nChannelsFirstWavelength = length(chanList)/2;
+
+    firstWavelengthChannel_idx = chanList(1:nChannelsFirstWavelength);
+    secondWavelengthChannel_idx = chanList(nChannelsFirstWavelength+1:end);
+    
+    %set secondWavelengthChannel_idx to zero if firstWavelengthChannel_idx
+    %is zero and vice versa
+    secondWavelengthChannel_idx(~firstWavelengthChannel_idx) = 0; 
+    firstWavelengthChannel_idx(~secondWavelengthChannel_idx) = 0;
+    
+    % combine both
+    chanList = [firstWavelengthChannel_idx; secondWavelengthChannel_idx];
+    
+    % sanity check - total_num_channels should be number of channel pairs
+    total_num_channels_check = sum(chanList(1:nChannelsFirstWavelength) == chanList(nChannelsFirstWavelength+1:end));
+    
     % update MeasListAct  
     mlActAuto{iBlk} = mlAct_Initialize(chanList, MeasList);
 end
diff --git a/FuncRegistry/UserFunctions/hmrR_StimRejection.m b/FuncRegistry/UserFunctions/hmrR_StimRejection.m
index 6125aa55..338e419b 100644
--- a/FuncRegistry/UserFunctions/hmrR_StimRejection.m
+++ b/FuncRegistry/UserFunctions/hmrR_StimRejection.m
@@ -53,7 +53,7 @@
     s = snirf.GetStims(t);
     
     dt = (t(end)-t(1))/length(t);
-    tRangeIdx = [floor(tRange(1)/dt):ceil(tRange(2)/dt)];
+    tRangeIdx = floor(tRange(1)/dt):ceil(tRange(2)/dt);
     
     smax = max(s,[],2);
     lstS = find(smax==1);
@@ -66,7 +66,7 @@
             s(lstS(iS),:) = -1*abs(s(lstS(iS),:));
         end
     end
-    snirf.SetStims_MatInput(s, t);
+    snirf.EditStates_MatInput(s, t);
 end
 stim = snirf.stim;
 
diff --git a/FuncRegistry/UserFunctions/hmrR_TaskNorm.m b/FuncRegistry/UserFunctions/hmrR_TaskNorm.m
new file mode 100644
index 00000000..f6321c73
--- /dev/null
+++ b/FuncRegistry/UserFunctions/hmrR_TaskNorm.m
@@ -0,0 +1,272 @@
+% SYNTAX:
+% [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+%
+% UI NAME:
+% Task_Normalization
+%
+% DESCRIPTION:
+% Calculate the average of the HRF over the duration of the task.
+% copyright  � 2023, Medea.
+% Authors: Anna Falivene, Charlotte Johnson
+
+%
+% INPUT:
+% dc: SNIRF.data container with the concentration data
+% stim: SNIRF.stim container with the stimulus condition data
+% deltat: defines the range before the onset and after the duration
+% [tPre_onset tPost_duration]. Both positive (or null) values
+% show: parameter to able/disable figure(). if show=0 no figure will
+% appear, if =1 figures with the HBO and HBR concentration for each
+% channel will appear
+%
+% OUTPUT:
+% dcAvg: averaged results
+% dcAvgStd: standard deviation across trials
+% nTrials: the number of trials averaged
+%
+% USAGE OPTIONS:
+% Task_Normalization_on_Concentration_Data: [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+%
+% PARAMETERS:
+% deltat: [0, 0]
+% show: 1
+%
+% PREREQUISITES:
+% Delta_OD_to_Conc: dc = hmrR_OD2Conc( dod, probe, ppf )
+%
+function [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+
+snirf = SnirfClass(dc, stim); %read the concentration data
+
+delta1=deltat(1); %baseline duration
+delta2=deltat(2); %post task duration
+dcAvg= DataClass();
+dcAvgStd = DataClass();
+d = dc.GetDataTimeSeries();
+
+%to compute the mean duration of the stimuli/tasks considered
+m_durations1=[];
+for ic=1:size(stim,2)
+    if size(stim(1, ic).data) ~=0
+
+        s1= stim(:,ic).GetData();
+        m_durations1=[m_durations1;mean(s1(:,2),1)];
+    end
+end
+
+for ic=1:size(stim,2)
+    if size(stim(1, ic).data) ~=0
+
+        s1= stim(:,ic).GetData();
+    else
+        continue
+    end
+    t = snirf.GetTimeCombined();
+    Fs=1/(t(2)-t(1)); %sampling rate
+    nTrials=size(s1,1);
+    namestr=[];
+   
+
+    for i=1:nTrials
+        my_field = string(['Repetition',num2str(i)]);
+         namestr=[namestr;my_field];
+        str.(my_field) = [];
+
+%data segmentation of each repetition
+        T1=find(t>=s1(i,1),1,'first');%finding onset
+        T2=find(t<=s1(i,1)+s1(i,2),1,'last');%finding end of trial (onset+dur)
+        newtime=t(T1:T2);
+        newdata=d(T1:T2,:); 
+        rsdata=(round(mean(m_durations1)))*50;%number of samples you want to have when interpolating/normalizing 
+        normtime=linspace(newtime(1),newtime(end),rsdata);
+        %data interpolation for the task execution phase of the trial
+        norm_data=interp1(newtime, newdata, normtime, 'spline');
+%data interpolation for the baseline and post-task phases of the trial
+        if delta1~=0 | delta2~=0
+
+            %data segmentation of baseline and post task phases
+            T1d=T1-round(delta1*Fs);
+            T2d=T2+round(delta2*Fs);
+            
+            rs1=(delta1*50);%defines numer of samples of delta1
+            rs2=(delta2*50);%defines numer of samples of delta2
+            n1=linspace(t(T1d),t(T1),rs1);%=Generate linearly spaced vector.
+            n2=linspace(t(T2),t(T2d),rs2);% =Generate linearly spaced vector.
+            if delta1~=0
+                d1=interp1(t(T1d:T1), d(T1d:T1,:), n1, 'spline');% interpolation of data in delta1
+            else
+                d1=[];
+            end
+            if delta2~=0
+                d2=interp1(t(T2:T2d), d(T2:T2d,:), n2, 'spline');% interpolation of data in delta2
+            else
+                d2=[];
+            end
+
+            norm_data=[ d1(1:end-1,:); norm_data; d2(2:end,:)];
+
+        end
+
+        str.(namestr(i,:))=norm_data;%saves new data in a structure
+    end
+
+    meanvalues=[]; stdvalues=[]; %forming a matrix where the MEAN and STD values are stored
+
+    ml = dc.GetMeasListSrcDetPairs('reshape'); %opening the measlist useful to calculate mean and stds 
+
+    %compute mean and std over the trials
+    for k=1:(3*size(ml,1))
+        ALLM=[];
+
+        for nt=1:nTrials
+            ALLM=[ALLM str.(namestr(nt,:))(:,k)];
+        end
+        meanvalues=[meanvalues mean(ALLM,2)];
+        stdvalues=[stdvalues std(ALLM,1,2)];
+
+    end
+
+    % '%cycle' vector
+    if delta1~=0 & delta2~=0
+        x1=linspace(0,100,rsdata+rs1+rs2-2);%=Generate linearly spaced vector from 0-100%.
+        disp (['onset: ', num2str(x1(rs1-1)), '%', ' - ', 'end of stimulus: ' ,num2str(x1(rs1-1+rsdata)), '%'])
+    elseif delta1~=0 & delta2==0 %only baseline
+        x1=linspace(0,100,rsdata+rs1-1);
+        disp (['onset: ', num2str(x1(rs1-1)), '%'])
+    elseif delta1==0 & delta2~=0 %no baseline but only resting period after the stim
+        x1=linspace(0,100,rsdata+rs2-1);
+        disp (['end of stimulus: ' ,num2str(x1(rs1-1+rsdata)), '%'])
+    else
+        x1=linspace(0,100,rsdata); %only exercise task
+    end
+
+    avghrf=[]; figs=[];
+    tot_avg=[];
+
+    %based on BlockAVG, Computing mean and std per channel in the different tasks.
+    for ll=1:3:size(meanvalues,2)-2
+
+        yHBO=meanvalues(:,ll); %HbO
+        yHBR=meanvalues(:,ll+1); %hbR
+        yHBT=meanvalues(:,ll+2); %hbt
+        sdHBO=stdvalues(:,ll);
+        sdHBR=stdvalues(:,ll+1);
+
+
+        ALLHBO=[];ALLHBR=[];%creating a matrix for All HbO and HbR data
+
+        for nt=1:nTrials
+%data concentration of each trial
+            ALLHBO=[ALLHBO str.(namestr(nt,:))(:,ll)];
+            ALLHBR=[ALLHBR str.(namestr(nt,:))(:,ll+1)];
+        end
+
+%baseline correction (as in Homer BlockAvg function) can be performed only if delta1 is different from 0
+        if delta1~=0
+            foomHBO = ones(size(yHBO,1),1)*mean(yHBO(1:rs1-1),1);
+            foomHBR = ones(size(yHBR,1),1)*mean(yHBR(1:rs1-1),1);
+            foomHBT = ones(size(yHBT,1),1)*mean(yHBT(1:rs1-1),1);
+            avghrf(:,ll)=yHBO-foomHBO;
+            avghrf(:,ll+1)=yHBR-foomHBR;
+            avghrf(:,ll+2)=yHBT-foomHBT;
+
+        else
+            %no correction
+            avghrf(:,ll)=yHBO;
+            avghrf(:,ll+1)=yHBR;
+            avghrf(:,ll+2)=yHBT;
+        end
+
+%%% Plot section 
+        f=figure();
+        ax1=subplot('Position',[0.35 0.58 0.33 0.34 ]);
+        txt=['Source ',num2str(dc.measurementList(1,ll).sourceIndex),' detector ', num2str(dc.measurementList(1,ll).detectorIndex)];
+
+        patch([x1 fliplr(x1)], [avghrf(:,ll)-sdHBO ;flipud(avghrf(:,ll)+sdHBO)], [0.6350 0.0780 0.1840] , 'FaceAlpha', 0.5 )
+        hold on
+        plot(x1,avghrf(:,ll),'r', 'LineWidth',2)
+
+        patch([x1 fliplr(x1)], [avghrf(:,ll+1)-sdHBR ;flipud(avghrf(:,ll+1)+sdHBR)], [0.3 0.75 0.9], 'FaceAlpha', 0.5)
+        plot(x1,avghrf(:,ll+1),'Color', '#0072BD', 'LineWidth',2)
+
+        hold off
+        lgd=legend('HbO','', 'HbR','');
+        title(txt )
+
+        ax2=subplot(223);
+        plot(x1,ALLHBO, 'LineWidth',1)
+        legend(namestr(:,:))
+        title('HbO_2 in each repetition')
+
+
+        ax3=subplot(224);
+        plot(x1,ALLHBR, 'LineWidth',1)
+        legend(namestr(:,:))
+        title('HbR in each repetition')
+
+        axx=[ax1 ax2 ax3];
+        grid (axx, 'on' )
+        ax1.YGrid='off';
+        xlabel(axx,'% of trial')
+        % ylabel(axx, '[M]')
+    
+        ax1.GridAlpha=0.5;
+        xtickangle([ax1 ax2 ax3],0)
+        
+        if delta1~=0 & delta2~=0
+            xticks([ax1 ax2 ax3],[0 x1(rs1-1) x1(rs1-1+rsdata) 100])
+            xticklabels([ax1 ax2 ax3],{'0','onset', 'end of the stimulus', '100'})
+
+
+        elseif delta1~=0 & delta2==0 %only baseline
+            xticks([ax1 ax2 ax3],[0 x1(rs1-1) 100])
+            xticklabels([ax1 ax2 ax3],{'0','onset', '100'})
+
+        elseif delta1==0 & delta2~=0 %no baseline but only resting period after the stim
+            xticks([ax1 ax2 ax3],[0  x1(end-(rs2-1)) 100])
+            xticklabels([ax1 ax2 ax3],{'0', 'end of the stimulus', '100'})
+
+        else %only exercises
+            xticks([ax1 ax2 ax3],[0  100])
+            xticklabels([ax1 ax2 ax3],{'0','100'})
+
+        end
+     
+        hlink = linkprop([ax1 ax2 ax3],{'YGrid','GridAlpha'});
+     
+
+        figs=[figs; f];
+    end
+    figs=(find(figs~=0));
+
+%saving figures
+    global maingui
+    folder=uigetdir();%decide where figures are stored.
+    namefig=maingui.dataTree.currElem.name; %can be modified according to how the .snirf file were named
+    namefig=namefig(1:end-6);
+    namefig=[namefig, '.fig'];
+    HRF=fullfile(folder,namefig);
+
+    savefig(figs,HRF,'compact');
+
+%add the data into the dcAvg and dcAvgStd Dataclasses
+  
+    for nn=1:size(ml,1)
+        dcAvg.AddChannelHbO(ml(nn,1), ml(nn,2), ic);
+        dcAvg.AddChannelHbR(ml(nn,1), ml(nn,2), ic);
+        dcAvg.AddChannelHbT(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbO(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbR(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbT(ml(nn,1), ml(nn,2), ic);
+        
+    end
+   
+    dcAvg.AppendDataTimeSeries(avghrf);
+    dcAvgStd.AppendDataTimeSeries(stdvalues);
+    dcAvg.SetTime(x1, true);
+    dcAvgStd.SetTime(x1, true);
+end
+
+
+
+
diff --git a/FuncRegistry/UserFunctions/mlAct_CombineIndexLists.m b/FuncRegistry/UserFunctions/mlAct_CombineIndexLists.m
new file mode 100644
index 00000000..bd872904
--- /dev/null
+++ b/FuncRegistry/UserFunctions/mlAct_CombineIndexLists.m
@@ -0,0 +1,9 @@
+function lstAct = mlAct_CombineIndexLists( mlAct1, mlAct2, ml )
+
+lstAct1 = mlAct_Matrix2IndexList(mlAct1, ml);
+lstAct2 = mlAct_Matrix2IndexList(mlAct2, ml);
+lstAct = [lstAct1(:)', lstAct2(:)'];
+
+k = find(ismember(lstAct,lstAct1) & ismember(lstAct,lstAct2)); 
+lstAct = unique(lstAct(k));
+
diff --git a/FuncRegistry/Version.txt b/FuncRegistry/Version.txt
index e21e727f..8e03717d 100644
--- a/FuncRegistry/Version.txt
+++ b/FuncRegistry/Version.txt
@@ -1 +1 @@
-1.4.0
\ No newline at end of file
+1.5.1
\ No newline at end of file
diff --git a/Homer3.m b/Homer3.m
index 6e00e17d..fddf7e9e 100644
--- a/Homer3.m
+++ b/Homer3.m
@@ -46,11 +46,11 @@
 end
 
 PrintSystemInfo(logger, 'Homer3', getArgs(groupDirs, inputFileFormat, unitTest, nargin));
-checkForHomerUpdates();
+checkForUpdates('Homer3');
 gdir = cfg.GetValue('Last Group Folder');
 if isempty(gdir)
     if isdeployed()
-        groupDirs = {[getAppDir(), 'SubjDataSample']};
+        groupDirs = {[getAppDir(), 'SampleData']};
     end
 end
 
diff --git a/Install/createInstallFile.m b/Install/createInstallFile.m
index 178f5ca8..46fc10c6 100644
--- a/Install/createInstallFile.m
+++ b/Install/createInstallFile.m
@@ -14,7 +14,7 @@ function createInstallFile(options)
 
 % Generate a LastCheckForUpdates.dat in case that we  haven't run Homer
 % before creating a build
-checkForHomerUpdates()
+checkForUpdates('Homer3')
 
 if ~exist('options','var') || isempty(options)
     options = 'all';
@@ -48,7 +48,7 @@ function createInstallFile(options)
 mkdir([dirnameInstall, installfilename]);
 mkdir([dirnameInstall, installfilename, '/FuncRegistry']);
 mkdir([dirnameInstall, installfilename, '/FuncRegistry/UserFunctions']);
-mkdir([dirnameInstall, installfilename, '/SubjDataSample']);
+mkdir([dirnameInstall, installfilename, '/SampleData']);
 
 % Generate executables
 if ~strcmp(options, 'nobuild')
@@ -86,8 +86,8 @@ function createInstallFile(options)
     end
 end
 
-if exist([dirnameApp, 'SubjDataSample'],'dir')
-    copyfile([dirnameApp, 'SubjDataSample'], [dirnameInstall, installfilename, '/SubjDataSample']);
+if exist([dirnameApp, 'SampleData'],'dir')
+    copyfile([dirnameApp, 'SampleData'], [dirnameInstall, installfilename, '/SampleData']);
 end
 
 for ii=1:length(platform.createshort_script)
diff --git a/Install/createShortcut.bat b/Install/createShortcut.bat
index 6093aef5..af5e8869 100644
--- a/Install/createShortcut.bat
+++ b/Install/createShortcut.bat
@@ -1,6 +1,10 @@
 
 @echo ====================================
-@echo Inside createShortcut.bat %1 %2 
+if not "%~3"=="" (
+    @echo Inside createShortcut.bat %1 %2 %3
+) else if not "%~2"==""  (
+    @echo Inside createShortcut.bat %1 %2
+)
 
 @SET DESKTOP_REG_ENTRY="HKCU\Software\Microsoft\Windows\CurrentVersion\Explorer\Shell Folders"
 @SET DESKTOP_REG_KEY="Desktop"
@@ -16,8 +20,13 @@
 @echo IF EXIST "%DESKTOP_DIR%"\%2.lnk (del /Q "%DESKTOP_DIR%"\%2.lnk)
 @IF EXIST "%DESKTOP_DIR%"\%2.lnk (del /Q "%DESKTOP_DIR%"\%2.lnk)
 
-@echo cscript "%dirnameSrc%"\createShortcut.vbs %1\%2
-@cscript "%dirnameSrc%"\createShortcut.vbs %1\%2
+if not "%~3"=="" (
+    @echo cscript "%dirnameSrc%"\createShortcut.vbs %1\%2\%3
+    @cscript "%dirnameSrc%"\createShortcut.vbs %1\%2\%3
+) else if not "%~2"==""  (
+    @echo cscript "%dirnameSrc%"\createShortcut.vbs %1\%2
+    @cscript "%dirnameSrc%"\createShortcut.vbs %1\%2
+) 
 
 @echo IF EXIST %1\%2.lnk (move %1\%2.lnk "%DESKTOP_DIR%"\%2.lnk)
 @IF EXIST %1\%2.lnk (move %1\%2.lnk "%DESKTOP_DIR%"\%2.lnk)
diff --git a/Install/createShortcut.vbs b/Install/createShortcut.vbs
index 94bbaafe..d50820c2 100644
--- a/Install/createShortcut.vbs
+++ b/Install/createShortcut.vbs
@@ -1,10 +1,15 @@
 
 Set objArgs = WScript.Arguments
 application = objArgs(0)
+startinfolder = ""
+if objArgs.Count > 1 then
+    startinfolder = objArgs(1)
+end if
+
 
 Set sh = CreateObject("WScript.Shell")
 Set shortcut = sh.CreateShortcut(application + ".lnk")
 
 shortcut.TargetPath = application
-shortcut.WorkingDirectory = ""
+shortcut.WorkingDirectory = startinfolder
 shortcut.Save
diff --git a/Install/setup.m b/Install/setup.m
index f68b70dd..4a7541a8 100644
--- a/Install/setup.m
+++ b/Install/setup.m
@@ -54,14 +54,31 @@ function main()
 nSteps = 100;
 iStep = 1;
 
-fprintf('dirnameSrc = %s\n', dirnameSrc)
-fprintf('dirnameDst = %s\n', dirnameDst)
-
 logger = Logger([dirnameSrc, 'Setup']);
 
+logger.Write('dirnameSrc = %s\n', dirnameSrc)
+logger.Write('dirnameDst = %s\n', dirnameDst)
+logger.Write('\n')
+
+% Create destination folders
+try 
+    mkdir(dirnameDst);
+catch ME
+    msg{1} = sprintf('Error: Could not create installation folder. It might be in use by other applications.\n');
+    msg{2} = sprintf('Try closing and reopening file browsers or any other applications that might be using the\n');
+    msg{3} = sprintf('installation folder and then retry installation.');
+    MenuBox(msg, 'OK');
+    close(h);
+    rethrow(ME)
+end
+
+myCopyFile([dirnameSrc, 'Version.txt'],     dirnameDst);
+logger.Write('\n')
+
 [~, exename] = getAppname();
 
 v = getVernum(exename);
+
 logger.Write('==========================================\n');
 logger.Write('Setup script for %s v%s:\n', exename, v);
 logger.Write('==========================================\n\n');
@@ -78,17 +95,6 @@ function main()
 
 pause(2);
 
-% Create destination folders
-try 
-    mkdir(dirnameDst);
-catch ME
-    msg{1} = sprintf('Error: Could not create installation folder. It might be in use by other applications.\n');
-    msg{2} = sprintf('Try closing and reopening file browsers or any other applications that might be using the\n');
-    msg{3} = sprintf('installation folder and then retry installation.');
-    menu([msg{:}], 'OK');
-    close(h);
-    rethrow(ME)
-end
 
 % Get full paths for source and destination directories
 dirnameSrc = fullpath(dirnameSrc);
@@ -96,16 +102,16 @@ function main()
 
 % Copy files from source folder to destination installation folder
 for ii = 1:length(platform.exename)
-    copyFile([dirnameSrc, platform.exename{ii}], [dirnameDst, platform.exename{ii}]);
+    myCopyFile([dirnameSrc, platform.exename{ii}], [dirnameDst, platform.exename{ii}]);
 end
-copyFile([dirnameSrc, 'db2.mat'],           dirnameDst);
-copyFile([dirnameSrc, 'AppSettings.cfg'],   dirnameDst);
-copyFile([dirnameSrc, 'DataTree'],          [dirnameDst, 'DataTree']);
-copyFile([dirnameSrc, 'FuncRegistry'],      [dirnameDst, 'FuncRegistry']);
-copyFile([dirnameSrc, 'SubjDataSample'],    [dirnameDst, 'SubjDataSample']);
-copyFile([dirnameSrc, 'SDGcolors.csv'],     dirnameDst);
-copyFile([dirnameSrc, 'Version.txt'],     dirnameDst);
-copyFile([dirnameSrc, 'LastCheckForUpdates.dat'],     dirnameDst);
+myCopyFile([dirnameSrc, 'db2.mat'],           dirnameDst);
+myCopyFile([dirnameSrc, 'AppSettings.cfg'],   dirnameDst);
+myCopyFile([dirnameSrc, 'DataTree'],          [dirnameDst, 'DataTree']);
+myCopyFile([dirnameSrc, 'FuncRegistry'],      [dirnameDst, 'FuncRegistry']);
+myCopyFile([dirnameSrc, 'SampleData'],    [dirnameDst, 'SampleData']);
+myCopyFile([dirnameSrc, 'SDGcolors.csv'],     dirnameDst);
+myCopyFile([dirnameSrc, 'Version.txt'],     dirnameDst);
+myCopyFile([dirnameSrc, 'LastCheckForUpdates.dat'],     dirnameDst);
 
 % Create desktop shortcuts to Homer3
 createDesktopShortcuts(dirnameSrc, dirnameDst);
@@ -154,7 +160,7 @@ function main()
             if ispathvalid(dirnameSrc)
                 err = -1;
             end
-            copyFile(dirnameSrc0, dirnameSrc);
+            myCopyFile(dirnameSrc0, dirnameSrc);
         end
         rmdir_safe(sprintf('%sDesktop/Test/', homePageFullPath()));
         
@@ -170,70 +176,6 @@ function main()
 
 
 
-
-
-% -------------------------------------------------------------------
-function copyFile(src, dst, type)
-global h
-global nSteps
-global iStep
-global logger
-
-if ~exist('type', 'var')
-    type = 'file';
-end
-
-try
-    % If src is one of several possible filenames, then src to any one of
-    % the existing files.
-    if iscell(src)
-        for ii=1:length(src)
-            if ~isempty(dir(src{ii}))
-                src = src{ii};
-                break;
-            end
-        end
-    end
-    
-    assert(logical(exist(src, type)));
-    
-    % Check if we need to untar the file 
-    k = findstr(src,'.tar.gz'); %#ok<FSTR>
-    if ~isempty(k)
-        untar(src,fileparts(src));
-        src = src(1:k-1);
-    end
-    
-    % Copy file from source to destination folder
-    
-    logmsg = sprintf('Copying %s to %s\n', src, dst);
-    if isempty(logger)
-        fprintf(logmsg);
-    else
-        logger.Write(logmsg);
-    end
-    copyfile(src, dst);
-
-    if ~isempty(iStep)
-        waitbar(iStep/nSteps, h); iStep = iStep+1;
-    end
-    pause(1);
-catch ME
-    if ishandles(h)
-        close(h);
-    end
-    printStack(ME);
-    if iscell(src)
-        src = src{1};
-    end
-    menu(sprintf('Error: Could not copy %s to installation folder.', src), 'OK');
-    pause(5);
-    rethrow(ME);
-end
-
-
-
-
 % --------------------------------------------------------------
 function deleteShortcuts()
 global platform
@@ -265,7 +207,7 @@ function createDesktopShortcuts(dirnameSrc, dirnameDst)
         cmd = sprintf('call "%s\\createShortcut.bat" "%s" %s.exe', dirnameSrc(1:end-1), dirnameDst, exename);
         system(cmd);
         
-        cmd = sprintf('call "%s\\createShortcut.bat" "%s" SubjDataSample', dirnameSrc(1:end-1), dirnameDst(1:end-1));
+        cmd = sprintf('call "%s\\createShortcut.bat" "%s" SampleData', dirnameSrc(1:end-1), dirnameDst(1:end-1));
         system(cmd);
         
     elseif islinux()
diff --git a/MainGUI/MainGUI.fig b/MainGUI/MainGUI.fig
index decd8961..84d6f9f1 100644
Binary files a/MainGUI/MainGUI.fig and b/MainGUI/MainGUI.fig differ
diff --git a/MainGUI/MainGUI.m b/MainGUI/MainGUI.m
index c4cef348..68b7f5a9 100644
--- a/MainGUI/MainGUI.m
+++ b/MainGUI/MainGUI.m
@@ -357,12 +357,9 @@ function MainGUI_EnableDisablePlotEditMode(handles, val)
 % --------------------------------------------------------------------
 function [nFileSuccess, nFilesWarning, nFilesFailed] = WarningsReport(handles)
 global maingui
-warnings = maingui.dataTree.GetWarningsReport();
 [nFileSuccess, nFilesWarning, nFilesFailed] = maingui.dataTree.GetErrorStats();
-if ~isempty(warnings)
-    set(handles.MainGUI,'visible','on')
-    MessageBox(warnings, 'WARNINGS')
-    %set(handles.listboxGroupTree, 'foregroundcolor',maingui.errcolor)
+if sum([nFilesWarning, nFilesFailed]) > 0
+    set(handles.pushbuttonHideErrors,'visible','on')
 end
 
 
@@ -431,7 +428,7 @@ function DisplayGroupTree(handles)
         set(handles.textStatus, 'foregroundcolor',maingui.errcolor);
         if nFilesFailed > 0
             set(handles.listboxFilesErr, 'visible','on', 'enable','on','value',1, 'string',listboxFilesErr)
-            set(handles.pushbuttonHideErrors, 'visible','on');
+            set(handles.pushbuttonHideErrors, 'visible','on', 'tooltipstring','Show/Hide Error Report Window.');
             if nFileSuccess==0
                 set(handles.pushbuttonHideErrors, 'visible','off');
                 p1 = get(handles.listboxGroupTree, 'position');
@@ -645,6 +642,9 @@ function listboxGroupTree_Callback(hObject0, eventdata, handles)
 msg = sprintf('%s:  %s', maingui.dataTree.filesErr(idx).filename, ...
     maingui.dataTree.filesErr(idx).GetErrorMsg());
 fprintf('%s\n', msg);
+if isempty(maingui.handles)
+    maingui.handles.msgbox = [];
+end
 if ishandle(maingui.handles.msgbox)
     delete(maingui.handles.msgbox);
 end
@@ -884,7 +884,7 @@ function menuItemPlotProbeGUI_Callback(hObject, ~, handles)
 LaunchChildGuiFromMenu('PlotProbeGUI', hObject, GetDatatype(handles), maingui.condition);
 
 % --------------------------------------------------------------------
-function menuItemPlotProbe2_Callback(hObject, ~, handles)
+function menuItemPlotProbe2GUI_Callback(hObject, ~, handles)
 global maingui
 procElem = maingui.dataTree.currElem;
 % Derived data that we want to save in a Snirf file.
@@ -917,7 +917,7 @@ function menuItemPlotProbe2_Callback(hObject, ~, handles)
     maingui.dataTree.LoadCurrElem()
 end
 probe = procElem.acquired.probe;
-stim = procElem.acquired.stim;
+stim = procElem.procStream.input.acquired.stim;
 metaDataTags = procElem.acquired.metaDataTags;
 
 obj = SnirfClass(data, stim, probe, metaDataTags);
@@ -1033,6 +1033,7 @@ function UpdateChildGuis(handles)
     MessageBox('Could not load current processing element. Acquisition files might be outdated or corrupted');
     return;
 end
+set(handles.textFileWarning, 'string',maingui.dataTree.currElem.GetErrorMsg(), 'fontsize',8, 'fontweight','bold');
 
 DisplayAxesSDG(handles);
 hObject = DisplayData(handles, hObject);
@@ -1199,6 +1200,11 @@ function UpdateChildGuis(handles)
 % ----------------------------------------------------------------------------------
 function DisplayStim(handles, hAxes)
 global maingui
+
+if strcmp(get(handles.menuItemHideStims, 'checked'), 'on')
+    return;
+end
+
 dataTree = maingui.dataTree;
 procElem = dataTree.currElem;
 
@@ -1418,14 +1424,14 @@ function Update(varargin)
 
 
 % --------------------------------------------------------------------
-function menuItemResetGroupFolder_Callback(hObject, eventdata, handles)
+function menuItemResetGroupFolder_Callback(~, ~, handles)
 resetGroupFolder();
 DisplayGroupTree(handles);
 
 
 
 % --------------------------------------------------------------------
-function pushbuttonDataPanLeft_Callback(hObject, eventdata, handles)
+function pushbuttonDataPanLeft_Callback(hObject, ~, handles)
 global maingui
 procElem = maingui.dataTree.currElem;
 iCh0     = maingui.axesSDG.iCh;
@@ -1473,7 +1479,7 @@ function pushbuttonDataPanRight_Callback(hObject, eventdata, handles)
 
 
 % --------------------------------------------------------------------
-function pushbuttonDataResetView_Callback(hObject, eventdata, handles)
+function pushbuttonDataResetView_Callback(hObject, ~, handles)
 global maingui
 set(handles.checkboxFixRangeX, 'value',0);
 set(handles.checkboxFixRangeY, 'value',0);
@@ -1704,7 +1710,7 @@ function checkboxExcludeTime_Callback(hObject, ~, handles)
 
 
 % --------------------------------------------------------------------
-function checkboxExcludeStims_Callback(hObject, eventdata, handles)
+function checkboxExcludeStims_Callback(hObject, ~, handles)
 global maingui
 
 hAxesData = maingui.axesData.handles.axes;
@@ -1731,7 +1737,7 @@ function checkboxExcludeStims_Callback(hObject, eventdata, handles)
 
 
 % --------------------------------------------------------------------
-function ExcludeTime_ButtonDownFcn(hObject, eventdata, handles)
+function ExcludeTime_ButtonDownFcn(hObject, ~, handles)
 global maingui
 
 % Make sure the user clicked on the axes and not 
@@ -1765,7 +1771,7 @@ function ExcludeTime_ButtonDownFcn(hObject, eventdata, handles)
 
 
 % --------------------------------------------------------------------
-function ExcludeStims_ButtonDownFcn(hObject, eventdata, handles)
+function ExcludeStims_ButtonDownFcn(hObject, ~, handles)
 global maingui
 
 if ~strcmp(get(hObject,'type'),'axes')
@@ -1954,7 +1960,6 @@ function menuItemPowerSpectrum_Callback(~, ~, ~)
 
 % -------------------------------------------------------------------------------
 function togglebuttonMinimizeGUI_Callback(hObject, ~, handles)
-u0 = get(handles.MainGUI, 'units');
 k = [1.0, 1.0, 0.8, 0.8];
 p0 = get(handles.MainGUI, 'position');
 if strcmp(get(hObject, 'tooltipstring'), 'Minimize GUI Window')
@@ -2065,10 +2070,10 @@ function panProbeCallback(hObject, ~, handles)
 elseif get(hObject,'string')=='>'
     xlim( [xrange(1)+xd/5 xrange(2)+xd/5] );
     maingui.axesSDG.xlim = [xrange(1)+xd/5 xrange(2)+xd/5];
-elseif get(hObject,'string')=='/\'
+elseif strcmp(get(hObject,'string'), '/\')
     ylim( [yrange(1)+yd/5 yrange(2)+yd/5] );
     maingui.axesSDG.ylim = [yrange(1)+yd/5 yrange(2)+yd/5];
-elseif get(hObject,'string')=='\/'
+elseif strcmp(get(hObject,'string'), '\/')
     ylim( [yrange(1)-yd/5 yrange(2)-yd/5] );
     maingui.axesSDG.ylim = [yrange(1)-yd/5 yrange(2)-yd/5];
 end
@@ -2220,6 +2225,16 @@ function menuItemPowerSpectrum_Loglog_Callback(~, ~, ~) %#ok<*DEFNU>
 
 % ---------------------------------------------------------
 function pushbuttonHideErrors_Callback(hObject, ~, handles)
+global maingui
+warnings = maingui.dataTree.GetWarningsReport();
+if strcmpi(get(handles.listboxFilesErr, 'enable'), 'off')
+    if ~isempty(warnings)
+        set(handles.MainGUI,'visible','on')
+        MessageBox(warnings, 'WARNINGS');
+    end
+    return;
+end
+
 pos2 = get(handles.listboxFilesErr, 'position');
 pos1 = get(handles.listboxGroupTree, 'position');
 x = pos2(4);
@@ -2232,6 +2247,10 @@ function pushbuttonHideErrors_Callback(hObject, ~, handles)
     set(hObject, 'string','/\');
     set(handles.listboxGroupTree, 'position', [pos1(1), pos1(2)-x, pos1(3), pos1(4)+x]);
 end
+if ~isempty(warnings) && strcmpi(get(handles.listboxFilesErr, 'visible'), 'off')
+    set(handles.MainGUI,'visible','on')
+    MessageBox(warnings, 'WARNINGS');
+end
 
 
 
@@ -2251,3 +2270,16 @@ function menuItemExportStim_Callback(~, ~, ~)
 function menuItemReloadStim_Callback(hObject, ~, handles)
 global maingui
 maingui.dataTree.currElem.EditStim();
+
+
+
+
+% --------------------------------------------------------------------
+function menuItemHideStims_Callback(hObject, ~, handles)
+if strcmp(get(hObject, 'checked'), 'on')
+    set(hObject, 'checked','off');
+elseif strcmp(get(hObject, 'checked'), 'off')
+    set(hObject, 'checked','on');
+end
+DisplayData(handles, hObject);
+
diff --git a/SampleData/Brite Fingertapping 5-12-2020.snirf b/SampleData/Brite Fingertapping 5-12-2020.snirf
new file mode 100644
index 00000000..c5033459
Binary files /dev/null and b/SampleData/Brite Fingertapping 5-12-2020.snirf differ
diff --git a/SampleData/Brite Fingertapping 5-15-2020.snirf b/SampleData/Brite Fingertapping 5-15-2020.snirf
new file mode 100644
index 00000000..1582f3bd
Binary files /dev/null and b/SampleData/Brite Fingertapping 5-15-2020.snirf differ
diff --git a/SampleData/Brite Fingertapping 5-19-2020.snirf b/SampleData/Brite Fingertapping 5-19-2020.snirf
new file mode 100644
index 00000000..aa91a927
Binary files /dev/null and b/SampleData/Brite Fingertapping 5-19-2020.snirf differ
diff --git a/SubjDataSample/test_process.cfg b/SampleData/test_process.cfg
similarity index 100%
rename from SubjDataSample/test_process.cfg
rename to SampleData/test_process.cfg
diff --git a/SubjDataSample/test.nirs b/SubjDataSample/test.nirs
deleted file mode 100644
index 33d4b205..00000000
Binary files a/SubjDataSample/test.nirs and /dev/null differ
diff --git a/UnitTests/FindUnitTestsFolders.m b/UnitTests/FindUnitTestsFolders.m
index ad4d0d66..b4b9bf07 100644
--- a/UnitTests/FindUnitTestsFolders.m
+++ b/UnitTests/FindUnitTestsFolders.m
@@ -1,6 +1,6 @@
 function groupFolders = FindUnitTestsFolders()
 global logger
-global SCRAMBLE
+global SCRAMBLE_CHANNELS
 
 logger = InitLogger(logger);
 
@@ -22,7 +22,7 @@
     groupFolders{kk} = dirs(ii).name;
     logger.Write('Found unit test folder %s\n', pathfull);
     
-    if SCRAMBLE
+    if SCRAMBLE_CHANNELS
         ScrambleChannelsForGroup(pathfull);
     end
     kk = kk+1;
diff --git a/UnitTests/UnitTests_Init.m b/UnitTests/UnitTests_Init.m
index 968e4024..fc34103e 100644
--- a/UnitTests/UnitTests_Init.m
+++ b/UnitTests/UnitTests_Init.m
@@ -6,7 +6,7 @@ function UnitTests_Init(standalone, start, appname)
 global logger
 global cfg
 global maingui
-global SCRAMBLE
+global SCRAMBLE_CHANNELS
 global ERROR_ODDS_CONST
 
 if ~exist('standalone','var') || isempty(standalone)
@@ -44,7 +44,7 @@ function UnitTests_Init(standalone, start, appname)
     reg = RegistriesClass();
     reg.DeleteSaved();
     
-    SCRAMBLE         = true;
+    SCRAMBLE_CHANNELS  = true;
     ERROR_ODDS_CONST = generateErrorOddsConstant(0);    
 end
 
diff --git a/UnitTests/testDataFilesClass.m b/UnitTests/testDataFilesClass.m
new file mode 100644
index 00000000..3a477b4c
--- /dev/null
+++ b/UnitTests/testDataFilesClass.m
@@ -0,0 +1,133 @@
+function [files0, files1, files2] = testDataFilesClass(options)
+% Syntax:
+%       [files0, files1, files2] = testDataFilesClass()
+%       [files0, files1, files2] = testDataFilesClass(options)
+%
+% Examples:
+%       [files0, files1, files2] = testDataFilesClass()
+%       [files0, files1, files2] = testDataFilesClass('create')
+%       [files0, files1, files2] = testDataFilesClass('create:remove')
+%
+%
+if ~exist('options','var')
+    options = '';
+end
+datasetDir = './group1';
+
+if includes(options,'create')
+    if ~exist(datasetDir,'dir')
+        fprintf('mkdir(''%s'');\n', datasetDir)
+        mkdir(datasetDir);
+    end
+    fprintf('cd(''%s'');\n', datasetDir)
+    cd(datasetDir)
+end
+dirs = dir('sub-*');
+if isempty(dirs)
+    fprintf('generateSimDataset(pwd, 5, 3, 4);\n')
+    if isempty(which('AtlasViewerGUI'))
+        MenuBox('ERROR:   AtlasViewer has not been installed. You must add AtlasViewer to Matlab search paths to generate simulated dataset with generateSimDataset.m')
+        return
+    end
+    generateSimDataset(pwd, 5, 3, 4);
+end
+fprintf('\n')
+pause(2)
+
+
+if exist('./derivatives','dir')
+    fprintf('rmdir(''./derivatives'',''s'')\n');
+    rmdir('./derivatives','s')
+end
+fprintf('tic; files0 = DataFilesClass(pwd, ''snirf''); toc;\n');
+tic; files0 = DataFilesClass(pwd, 'snirf'); toc %#ok<*NASGU>
+fprintf('\n')
+pause(2)
+
+% Update files
+filesUpdated = {
+    './sub-02/ses-01/sub-02_ses-01_run02.snirf';
+    './sub-02/ses-03/sub-02_ses-03_run01.snirf';
+    './sub-03/ses-02/sub-03_ses-02_run02.snirf';
+    './sub-03/ses-02/sub-03_ses-02_run04.snirf';
+    './sub-05/ses-01/sub-05_ses-01_run04.snirf';
+    };
+s = repmat(SnirfClass(),1,length(filesUpdated));
+for ii = 1:length(filesUpdated)
+    fprintf('Update file %s\n', filesUpdated{ii})
+    s(ii) = SnirfClass(filesUpdated{ii});
+    s(ii).Save();
+end
+fprintf('\n')
+pause(2)
+
+
+% Added files
+filesAdded = {
+    './sub-01/ses-02/sub-01_ses-02_run10.snirf';
+    './sub-03/ses-02/sub-04_ses-02_run14.snirf';
+    './sub-04/ses-02/sub-04_ses-02_run11.snirf';
+    './sub-04/ses-03/sub-04_ses-03_run12.snirf';
+    };
+for ii = 1:length(filesAdded)
+    fprintf('Added file %s\n', filesAdded{ii})
+    s = SnirfClass('./sub-01/ses-01/sub-01_ses-01_run01.snirf');
+    s.Save(filesAdded{ii});
+end
+fprintf('\n')
+pause(2)
+
+
+% Deleted files
+filesDeleted = {
+    './sub-02/ses-03/sub-02_ses-03_run02.snirf';
+    './sub-04/ses-02/sub-04_ses-02_run01.snirf';
+    './sub-05/ses-02/sub-05_ses-02_run03.snirf';
+    };
+for ii = 1:length(filesDeleted)
+    fprintf('Delete file %s\n', filesDeleted{ii})
+    delete(filesDeleted{ii})
+end
+fprintf('\n')
+pause(2)
+
+
+
+% Error files
+filesError = {
+    './sub-01/ses-01/sub-02_ses-01_run02.snirf';
+    './sub-04/ses-03/sub-02_ses-03_run01.snirf';
+    './sub-05/ses-02/sub-03_ses-02_run02.snirf';
+    './sub-05/ses-02/sub-03_ses-02_run04.snirf';
+    };
+s = repmat(SnirfClass(),1,length(filesError));
+for ii = 1:length(filesError)
+    fprintf('Generate error file %s\n', filesError{ii})
+    s(ii) = SnirfClass(filesError{ii});
+    s(ii).data.dataTimeSeries(:,1) = [];
+    s(ii).Save();
+end
+fprintf('\n')
+pause(2)
+
+
+
+fprintf('tic; files1 = DataFilesClass(pwd, ''snirf''); toc;\n');
+tic; files1 = DataFilesClass(pwd, 'snirf'); toc
+fprintf('\n')
+pause(2)
+
+fprintf('tic; files2 = DataFilesClass(pwd, ''snirf''); toc;\n');
+tic; files2 = DataFilesClass(pwd, 'snirf'); toc
+fprintf('\n')
+pause(2)
+
+if includes(options,'remove')
+    fprintf('cd(''..'');\n')
+    cd('..')
+    fprintf('rmdir(''%s'',''s'');\n', datasetDir)
+    rmdir(datasetDir,'s')
+    fprintf('\n')
+end
+
+
diff --git a/Utils/Shared/ConfigFileClass.m b/Utils/Shared/ConfigFileClass.m
index 9dedf9ec..b46e9cb1 100644
--- a/Utils/Shared/ConfigFileClass.m
+++ b/Utils/Shared/ConfigFileClass.m
@@ -661,6 +661,9 @@ function FindCfgFiles(obj, rootdir)
                 if strcmp(dirs(ii).name, 'submodules')
                     continue;
                 end
+                if strcmp(dirs(ii).name, 'Install')
+                    continue;
+                end
                 obj.FindCfgFiles([rootdir, dirs(ii).name]);
             end
         end
diff --git a/Utils/Shared/Logger.m b/Utils/Shared/Logger.m
index 71b420ae..97a8fd4e 100644
--- a/Utils/Shared/Logger.m
+++ b/Utils/Shared/Logger.m
@@ -194,6 +194,20 @@ function WriteNoFmt(self, s, options, hwait)
         
                 
         
+        % -------------------------------------------------
+        function WriteFileOnly(self, s, varargin)
+            optionsValueOrig = self.options.value;
+            self.SetFileOnly();
+            if isempty(varargin)
+                self.Write(s, varargin);
+            else
+                self.Write(s, varargin{:});
+            end
+            self.options.value = optionsValueOrig;
+        end
+        
+        
+        
         % -------------------------------------------------
         function ct = CurrTime(self, msg, options, hwait)
             ct = '';
@@ -363,22 +377,22 @@ function SetDebugLevel(self, options)
         
         
         % -------------------------------------------------
-        function val = FileOnly(self)
-            val = self.options.FILE_ONLY;
+        function SetFileOnly(self)
+            self.options.value = self.options.FILE_ONLY;
         end
         
         
         
         % -------------------------------------------------
-        function val = Null(self)
-            val = self.options.NULL;
+        function SetConsoleOnly(self)
+            self.options.value = self.options.CONSOLE_ONLY;
         end
         
         
         
         % -------------------------------------------------
-        function val = ConsoleOnly(self)
-            val = self.options.CONSOLE_ONLY;
+        function val = Null(self)
+            val = self.options.NULL;
         end
         
         
diff --git a/Utils/Shared/Version.txt b/Utils/Shared/Version.txt
index 081af9a1..fe4e75fb 100644
--- a/Utils/Shared/Version.txt
+++ b/Utils/Shared/Version.txt
@@ -1 +1 @@
-1.7.1
\ No newline at end of file
+1.8.3
\ No newline at end of file
diff --git a/Utils/Shared/ceiln.m b/Utils/Shared/ceiln.m
new file mode 100644
index 00000000..e51755fd
--- /dev/null
+++ b/Utils/Shared/ceiln.m
@@ -0,0 +1,8 @@
+function val = ceiln(val0, n)
+if ~exist('n','var')
+    n = 0;
+end
+pw = ceil(log10(abs(val0))); % Find order of magnitude of val.
+res = 10^(pw-n-2); % Resolution to round to.
+valtemp = val0*res;
+val = ceil(valtemp)/res; % < change floor() to ceil(), for ceiling equivalent.
diff --git a/Utils/Shared/checkForUpdates.m b/Utils/Shared/checkForUpdates.m
new file mode 100644
index 00000000..1540fac1
--- /dev/null
+++ b/Utils/Shared/checkForUpdates.m
@@ -0,0 +1,65 @@
+function checkForUpdates(appname)
+global cfg
+
+cfg = InitConfig(cfg);
+
+% If user has Check For Updates enabled
+if (strcmp(cfg.GetValue('Check For Updates'),'on'))
+    
+    % If it has been a week since Homer checked for an update
+    if (datetime - GetLastCheckForUpdates() > duration(168,0,0))
+        url = sprintf('https://openfnirs.org/software/homer/%s/', appname);
+        [s,status] = urlread(url,'timeout',4);        
+        if (~status)
+            % App is offline or server could not be reached
+            fprintf('Server could not be reached to check for updates.')
+            return
+        end
+        SetLastCheckForUpdates();
+        version = regexp(s, 'id="version">(.*?)<\/', 'tokens');
+        if isempty(version)
+            return
+        end
+        web_vrnum = version{1}{1};
+        this_vrnum = getVernum();
+        promptFlag = compareVernum(web_vrnum, this_vrnum);  % If fetched vernum is greater
+
+        % Open a hidden web browser        
+        if (promptFlag)            
+            wb = com.mathworks.mde.webbrowser.WebBrowser.createBrowser;
+            wb.setCurrentLocation(url);
+            p = getParentRecursive(wb);
+            p.setVisible(0);            
+            choice = questdlg(sprintf('Version of your copy of %s is v%s\nA newer version v%s is available.\nWould you like to download?', appname, this_vrnum, web_vrnum), ...
+                             'Update Available', ...
+                             'Yes', ...
+                             'Remind me later', ...
+                             'Don''t show this again', ...
+                             'Remind me later');            
+            if strcmp(choice, 'Yes')
+                web('https://github.com/BUNPC/AtlasViewer/releases');
+                close(wb);
+            elseif strcmp(choice, 'Don''t show this again')
+                web('https://github.com/BUNPC/AtlasViewer');
+                cfg.SetValue('Check For Updates', 'off');
+            end            
+            close(wb);
+        end        
+        pause(1);  % To ensure <script> is run
+        
+    end
+    
+end
+
+
+
+% ------------------------------------------------------------------------------
+function v1_greater = compareVernum(v1, v2)
+v1_greater = false;
+v1num = versionstr2num(v1);
+v2num = versionstr2num(v2);
+if v1num > v2num
+    v1_greater = true;
+end    
+
+
diff --git a/Utils/Shared/distmatrix.m b/Utils/Shared/distmatrix.m
new file mode 100644
index 00000000..d086c143
--- /dev/null
+++ b/Utils/Shared/distmatrix.m
@@ -0,0 +1,49 @@
+% Usage:
+%
+%     d = distmatrix(varargin)
+%
+% Example 1:
+%
+%     d = distmatrix(v);
+%
+% Example 2:
+%
+%     d = distmatrix(v1,v2);
+%
+%
+% Author: Jay Dubb (jdubb@nmr.mgh.harvard.edu)
+% Date:   11/13/2008
+%  
+
+function d = distmatrix(varargin)
+
+d=[];
+if length(varargin)==1
+    
+    v = double(varargin{1});
+    n = size(v,1);
+    d = zeros(n);
+    for i=1:n
+        for j=i+1:n        
+            d(i,j) = ((v(i,1) - v(j,1))^2 + ...
+                      (v(i,2) - v(j,2))^2 + ...
+                      (v(i,3) - v(j,3))^2)^0.5;
+        end     
+    end
+
+elseif length(varargin)==2
+
+    v1 = double(varargin{1});
+    v2 = double(varargin{2});
+    m = size(v1,1);
+    n = size(v2,1);
+    d = zeros(m,n);    
+    for i=1:m
+        for j=1:n
+            d(i,j) = ((v1(i,1) - v2(j,1))^2 + ...
+                      (v1(i,2) - v2(j,2))^2 + ...
+                      (v1(i,3) - v2(j,3))^2)^0.5;
+        end     
+    end    
+
+end
diff --git a/Utils/Shared/floorn.m b/Utils/Shared/floorn.m
new file mode 100644
index 00000000..a7f50ab8
--- /dev/null
+++ b/Utils/Shared/floorn.m
@@ -0,0 +1,9 @@
+function val = floorn(val0, n)
+if ~exist('n','var')
+    n = 0;
+end
+pw = ceil(log10(abs(val0))); % Find order of magnitude of val.
+res = 10^(pw-n-2); % Resolution to round to.
+valtemp = val0*res;
+val = floor(valtemp)/res; % < change floor() to ceil(), for ceiling equivalent.
+
diff --git a/Utils/Shared/getAppDir.m b/Utils/Shared/getAppDir.m
index dbde5abf..e5a1d022 100644
--- a/Utils/Shared/getAppDir.m
+++ b/Utils/Shared/getAppDir.m
@@ -42,6 +42,9 @@
         dirname = [dirnameHome, '/atlasviewer/'];
         cd(currdir);
     end
+    if ~ispathvalid(dirname)
+        mkdir(dirname);
+    end
 else
     dirname = fileparts(which('AtlasViewerGUI.m'));
 end
@@ -72,10 +75,15 @@
         dirname = [dirnameHome, '/homer3/'];
         cd(currdir);
     end
+    if ~ispathvalid(dirname)
+        mkdir(dirname);
+    end
 else
     dirname = fileparts(which('Homer3.m'));
 end
-
+if ~ispathvalid(dirname)
+    dirname = filesepStandard(pwd);
+end
 dirname(dirname=='\') = '/';
 if dirname(end) ~= '/'
     dirname(end+1) = '/';
@@ -100,6 +108,9 @@
         dirname = [dirnameHome, '/datatree/'];
         cd(currdir);
     end
+    if ~ispathvalid(dirname)
+        mkdir(dirname);
+    end
 else
     dirname = fileparts(which('DataTreeClass.m'));
 end
diff --git a/Utils/Shared/pretty_print_struct.m b/Utils/Shared/pretty_print_struct.m
index bb87eb2a..79fd0e4d 100644
--- a/Utils/Shared/pretty_print_struct.m
+++ b/Utils/Shared/pretty_print_struct.m
@@ -1,4 +1,4 @@
-function pretty_print_struct(st, indent, option, logger)
+function pretty_print_struct(st, indent, option)
 spaces = '';
 
 if ~exist('st','var') || isempty(st)
@@ -10,11 +10,6 @@ function pretty_print_struct(st, indent, option, logger)
 if ~exist('option','var') || isempty(option)
     option = 1;
 end
-if ~exist('logger','var')
-    logger = [];
-end
-logger = InitLogger(logger);
-
 if iswholenum(indent)
     spaces = blanks(indent);
 elseif ischar(indent)
@@ -22,6 +17,7 @@ function pretty_print_struct(st, indent, option, logger)
 end
 
 if isstruct(st) || isobject(st)
+    fields = properties(st);    
     s = evalc('disp(st)');
     c = str2cell_fast(s, char(10));
     for ii=1:length(c)
diff --git a/Utils/Shared/resamplee.m b/Utils/Shared/resamplee.m
new file mode 100644
index 00000000..6555e018
--- /dev/null
+++ b/Utils/Shared/resamplee.m
@@ -0,0 +1,10 @@
+function resig = resamplee(sig,upsample,downsample)
+if upsample*downsample<2^31
+    resig = resample(sig,upsample,downsample);
+else
+    sig1half = sig(1:floor(length(sig)/2));
+    sig2half = sig(floor(length(sig)/2):end);
+    resig1half = resamplee(sig1half, floor(upsample/2), length(sig1half));
+    resig2half = resamplee(sig2half, upsample-floor(upsample/2), length(sig2half));
+    resig = [resig1half(:); resig2half(:)];
+end
diff --git a/Utils/Shared/simData/generateSimDataset.m b/Utils/Shared/simData/generateSimDataset.m
index 8d7bf8db..e73847be 100644
--- a/Utils/Shared/simData/generateSimDataset.m
+++ b/Utils/Shared/simData/generateSimDataset.m
@@ -1,28 +1,93 @@
-function nirs = generateSimDataset(dirname, nSubj, nSess, nRuns, options)
+function out = generateSimDataset(dirname, nSubj, nSess, nRuns, options)
+%
+% Syntax:
+%   out = generateSimDataset()
+%   out = generateSimDataset(dirname)
+%   out = generateSimDataset(dirname, nSubj)
+%   out = generateSimDataset(dirname, nSubj, nSess)
+%   out = generateSimDataset(dirname, nSubj, nSess, nRuns)
+%   out = generateSimDataset(dirname, nSubj, nSess, nRuns, options)
+%
+% Description:
+%   options - { {'probe'|'digpts}, 'springs', 'data' }
+%
+% Examples:
+%
+%   %%%%% Change current folder to <root groups test folder>
+%   emptyGroupTestDir = 'c:/users/public/groups/emptyGroupTestDir';
+%   if ~exist(emptyGroupTestDir,'dir'), mkdir(emptyGroupTestDir); end
+%   cd(emptyGroupTestDir);
+%
+%   %%%%% Now generate probe and/or digitized points and/or raw/unprocessed dataTree data.
+%   %%%%% These 9 examples covers all the major AtlasViewer + Homer3 workflows
+%
+%   % Example 1: probe file in SD format with spring registration and
+%   % raw/unprocessed dataTree data with 3 subjects, 3 sessions each, and 3 runs each.
+%   out = generateSimDataset(pwd, 3, 3, 3, 'probe:springs:data');
+%
+%   % Example 2: probe file in SD format with spring registration and
+%   % raw/unprocessed dataTree data with 3 subjects, 3 sessions each with 3 runs each.
+%   out = generateSimDataset();
+%
+%   % Example 3: probe file in SD format with spring registration and
+%   % acquisition data with 3 subjects, 1 session each, and 3 runs each.
+%   out = generateSimDataset(pwd);
+%
+%   % Example 4: probe file in SD format with spring registration and acquisition data
+%   % with 2 subjects with 1 session each, and each session with 5 runs each.
+%   out = generateSimDataset(pwd, 2, 1, 5);
+%
+%   % Example 5: probe file in SD format with NO registration and NO acquisition data
+%   out = generateSimDataset(pwd, [],[],[],'probe');
+%
+%   % Example 6: probe file in SD format with spring registration and NO data
+%   out = generateSimDataset(pwd, [],[],[],'probe:springs');
+%
+%   % Example 7: Unprocessed dataTree data with 3 subjects, 3 sessions, and 3 runs each
+%   % and with acquisition files containing probe with landmark registration.
+%   % with 3 subjects with 1 session each, and each session with 3 runs each.
+%   out = generateSimDataset(pwd, [],[],[],'probe:landmarks:data');
+%
+%   % Example 8: Digpts file with associated probe file in SD format containing measurement list
+%   % and NO acquisition data
+%   out = generateSimDataset(pwd, [],[],[],'digpts');
+%
+%   % Example 9: Digpts file with associated probe file in SD format containing measurement list
+%   % and data with 1 subject with 2 sessions each, and each session with 4 runs each.
+%   out = generateSimDataset(pwd, 1, 2, 4, 'digpts:data');
+%
+%
+global atlasViewer
+global probeFilename
 
-% 
+atlasViewer = [];
+probeFilename = 'newfile1.SD';
 
-nirs = [];
+out = [];
 
 t0 = tic;
 
-
-if ~exist('dirname','var')
-    dirname = filesepStandard(pwd);
+if ~exist('dirname','var') || isempty(dirname)
+    dirname = pwd;
 end
-if ~exist('nSubj','var')
+if ~exist('nSubj','var') || isempty(nSubj)
     nSubj = 3;
 end
-if ~exist('nSess','var')
-    nSess = 3;
+if ~exist('nSess','var') || isempty(nSess)
+    nSess = 1;
 end
-if ~exist('nRuns','var')
+if ~exist('nRuns','var') || isempty(nRuns)
     nRuns = 3;
 end
 if ~exist('options','var')
-    options = 'probe:springs';
+    options = 'probe:springs:data';
 end
 
+dirname = filesepStandard(dirname);
+
+% Reset folder to known state
+resetDataset();
+
 % Create template data
 setNamespace('AtlasViewerGUI');
 dirnameAtlas = getAtlasDir();
@@ -31,7 +96,38 @@
 end
 refpts = initRefpts();
 refpts = getRefpts(refpts, dirnameAtlas);
-SD = genProbeFromRefpts(refpts, 36, 2, options);
+
+if optionExists(options, 'probe')
+    genProbeFromRefpts(refpts, 36, 2, options);
+elseif optionExists(options, 'digpts')
+    genDigptsFromRefpts(refpts, 36, 2);
+end
+
+if optionExists(options, 'data')
+    digpts = initDigpts();
+    digpts = getDigpts(digpts);
+    probe  = initProbe();
+    probe  = getProbe(probe, dirname, digpts);
+    genSimData(probe, dirname, nSubj, nSess, nRuns, t0);
+    subjDirs = dir('sub-*');
+    for ii = 1:length(subjDirs)
+        digptsNew = movePts(digpts, randNearZero(1,3,t0), randNearOne(1,3,t0), randNearZero(1,3,t0));
+        digptsNew.pathname = [filesepStandard(dirname), subjDirs(ii).name];
+        saveDigpts(digptsNew);
+        if ispathvalid([dirname, probeFilename])
+            if optionExists(options, 'springs')
+                copyfile([dirname, probeFilename], [filesepStandard(dirname), subjDirs(ii).name, '/', probeFilename]);
+            end
+        end
+    end
+    delete('./digpts*.txt');
+end
+
+
+
+% ----------------------------------------------
+function genSimData(probe, dirname, nSubj, nSess, nRuns, t0)
+SD = convertProbe2SD(probe);
 nirs = NirsClass(SD);
 
 % Generate stims
@@ -48,27 +144,166 @@
     offset = offset+50;
 end
 
-
-%
+% Generate data
 for iSubj = 1:nSubj
-    sname = sprintf('subj-%d', iSubj);
-    if ispathvalid([dirname, sname])
-        rmdir([dirname, sname], 's')
+    if iSubj < 10
+        iSubjName = sprintf('0%d', iSubj);
+    else
+        iSubjName = sprintf('%d', iSubj);
     end
-    mkdir([dirname, sname]);
-
-    for iRun = 1:nRuns
-        for iM = 1:size(nirs.SD.MeasList,1)
-            [nirs.t, nirs.d(:,iM)] = simulateDataTimeSeries(ntpts, 1, .4, t0);
-        end
-        snirf = SnirfClass(nirs.d, nirs.t, nirs.SD, [], nirs.s, nirs.CondNames);
-        if iRun<10
-            rname = sprintf('%s%s/%s_run0%d.snirf', dirname, sname, sname, iRun);
+    subjName = sprintf('sub-%s', iSubjName);
+    subjNameFull = filesepStandard([dirname, subjName], 'dir:nameonly');
+    if ispathvalid(subjNameFull)
+        rmdir(subjNameFull, 's')
+    end
+    mkdir(subjNameFull);
+    
+    for iSess = 1:nSess
+        if iSess < 10
+            iSessName = sprintf('0%d', iSess);
         else
-            rname = sprintf('%s%s/%s_run%d.snirf', dirname, sname, sname, iRun);
+            iSessName = sprintf('%d', iSess);
+        end
+        sessName = sprintf('ses-%s', iSessName);
+        sessNameFull = filesepStandard([subjNameFull, sessName], 'dir:nameonly');
+        if ispathvalid(sessNameFull)
+            rmdir(sessNameFull, 's')
         end
-        snirf.Save(rname);
-        fprintf('Created run %s\n', rname);
+        mkdir(sessNameFull);
+        
+        for iRun = 1:nRuns
+            for iM = 1:size(nirs.SD.MeasList,1)
+                [nirs.t, nirs.d(:,iM)] = simulateDataTimeSeries(ntpts, 1, .4, t0);
+            end
+            snirf = SnirfClass(nirs.d, nirs.t, nirs.SD, [], nirs.s, nirs.CondNames);
+            if iRun<10
+                runNameFull = sprintf('%s%s_%s_run0%d.snirf', sessNameFull, subjName, sessName, iRun);
+            else
+                runNameFull = sprintf('%s%s_%s_run%d.snirf', sessNameFull, subjName, sessName, iRun);
+            end
+            snirf.Save(runNameFull);
+            fprintf('Created run %s\n', runNameFull);
+        end
+    end
+end
+
+
+% -----------------------------------------------------------------
+function resetDataset(dirname)
+if ~exist('dirname','var')
+    dirname = filesepStandard(pwd);
+end
+
+currdir = filesepStandard(pwd);
+cd(dirname);
+
+files = dir('*');
+for ii = 1:length(files)
+    if strcmpi(files(ii).name, '.')
+        continue
+    end
+    if strcmpi(files(ii).name, '..')
+        continue
+    end
+    if strcmpi(files(ii).name, 'anatomical')
+        continue
+    end
+    if files(ii).isdir
+        rmdir(files(ii).name,'s')
+    else
+        delete(files(ii).name)
+    end
+end
+cd(currdir);
+
+
+
+
+% -----------------------------------------------
+function digpts = movePts(digpts, r, s, t)
+if ~exist('r','var')
+    r = [0,0,0];
+end
+if ~exist('s','var')
+    s = [1,1,1];
+end
+if ~exist('t','var')
+    t = [0,0,0];
+end
+
+alpha = deg2rad(r(1));
+beta  = deg2rad(r(2));
+theta = deg2rad(r(3));
+
+A = [ ...
+    1            0             0    0;
+    0   cos(alpha)   -sin(alpha)    0;
+    0   sin(alpha)    cos(alpha)    0;
+    0            0             0    1;
+    ];
+
+B = [ ...
+    cos(beta)    0     sin(beta)     0;
+    0            1             0     0;
+    -sin(beta)    0     cos(beta)     0;
+    0            0             0     1;
+    ];
+
+C = [ ...
+    cos(theta)   -sin(theta)   0     0;
+    sin(theta)    cos(theta)   0     0;
+    0            0             1     0;
+    0            0             0     1;
+    ];
+
+D = [ ...
+    1            0             0     t(1);
+    0            1             0     t(2);
+    0            0             1     t(3);
+    0            0             0     1;
+    ];
+
+S = [ ...
+    s(1)           0             0     0;
+    0           s(2)           0     0;
+    0            0           s(3)    0;
+    0            0             0     1;
+    ];
+
+T = S*D*C*B*A;
+
+digpts.srcpos = xform_apply(digpts.srcpos, T);
+digpts.detpos = xform_apply(digpts.detpos, T);
+digpts.refpts.pos = xform_apply(digpts.refpts.pos, T);
+
+
+
+% ----------------------------------------------------------------
+function r = randNearOne(w, h, t0)
+generateRandNumSeed(t0);
+r = ones(w,h) + (rand(w,h) - 0.5) / 2;
+
+
+
+% ----------------------------------------------------------------
+function r = randNearZero(w, h, t0)
+generateRandNumSeed(t0);
+r = (rand(w,h) - 0.5) / 2;
+
+
+
+% ---------------------------------------------------
+function generateRandNumSeed(time0)
+if time0 == 0
+    x = uint32(100*rand);
+    rng(x);
+    y = uint32(100*rand);
+    rng(y);
+else
+    s = 0;
+    while s==0
+        s = uint64(1e4*toc(time0));
     end
+    rng(s);
 end
- 
+
diff --git a/Utils/submodules/syncSubmodules.m b/Utils/submodules/syncSubmodules.m
index fa1679fa..e3c018b1 100644
--- a/Utils/submodules/syncSubmodules.m
+++ b/Utils/submodules/syncSubmodules.m
@@ -79,7 +79,6 @@
 end
 fprintf('\n');
 
-cmds{ii,1} = sprintf('cd %s', synctool.repoParentFull); ii = ii+1;
 for jj = 1:size(submodules,1)
     submodulename = getRepos(submodules, jj, options);
     
diff --git a/Version.txt b/Version.txt
index 09e7584a..59d7d104 100644
--- a/Version.txt
+++ b/Version.txt
@@ -1 +1 @@
-1.80.2
\ No newline at end of file
+1.86.0
\ No newline at end of file