Merge pull request #9 from xdf-modules/fix_robust_fit

Fix robust_fit issue with large timestamps

Author: cboulay

load_xdf.m

@@ -73,6 +73,10 @@
 %                                                     ClockOffset difference that is at least this
 %                                                     many seconds away from the median. (default: 1)
 %
+%                'WinsorThreshold' : Error threshold beyond which clock offset deviations are considered
+%                                    outliers, and therefore penalized less severely (linear instead of quadratic).
+%                                    (default: 0.0001)
+%
 %                'ClockResetMaxJitter' : Maximum tolerable jitter (in seconds of error) for clock
 %                                        reset handling. (default: 5)
 %
@@ -365,7 +369,7 @@
             temp(id).time_series = {};
             temp(id).time_stamps = {};
             temp(id).clock_times = [];
-            temp(id).clock_values = [];
+            temp(id).offset_values = [];
             if temp(id).srate > 0
                 temp(id).sampling_interval = 1/temp(id).srate;
             else
@@ -391,7 +395,7 @@
                 % read [CollectionTime]
                 temp(id).clock_times(end+1) = fread(f,1,'double');
                 % read [OffsetValue]
-                temp(id).clock_values(end+1) = fread(f,1,'double');
+                temp(id).offset_values(end+1) = fread(f,1,'double');
             catch e
                 % an error occurred (perhaps a chopped-off file): emit a
                 % warning and scan forward to the next recognized chunk
@@ -431,7 +435,7 @@
         if ~isempty(temp(k).time_stamps)
             try
                 clock_times = temp(k).clock_times;
-                clock_values = temp(k).clock_values;
+                offset_values = temp(k).offset_values;
                 if isempty(clock_times)
                     error('No clock offset values present.'); end
             catch
@@ -447,7 +451,7 @@
                 % importer should be able to deal with recordings where the computer that served a stream
                 % was restarted or hot-swapped during an ongoing recording, or the clock was reset otherwise
                 time_diff = diff(clock_times);
-                value_diff = abs(diff(clock_values));
+                value_diff = abs(diff(offset_values));
                 % points where a glitch in the timing of successive clock measurements happened
                 time_glitch = (time_diff < 0 | (((time_diff - median(time_diff)) ./ median(abs(time_diff-median(time_diff)))) > opts.ClockResetThresholdStds & ...
                     ((time_diff - median(time_diff)) > opts.ClockResetThresholdSeconds)));
@@ -472,15 +476,18 @@
                 ranges = {[1,length(clock_times)]};
             end
 
-            % calculate clock offset mappings for each data range
+            % Calculate clock offset mappings for each data range
             mappings = {};
             for r=1:length(ranges)
                 idx = ranges{r};
                 if idx(1) ~= idx(2)
-                    % to accomodate the Winsorizing threshold (in seconds) we rescale the data (robust_fit sets it to 1 unit)
-                    mappings{r} = robust_fit([ones(idx(2)-idx(1)+1,1) clock_times(idx(1):idx(2))']/opts.WinsorThreshold, clock_values(idx(1):idx(2))'/opts.WinsorThreshold);
+                    Ax = clock_times(idx(1):idx(2))' / opts.WinsorThreshold;
+                    y = offset_values(idx(1):idx(2))'  / opts.WinsorThreshold;
+                    fit_params = robust_fit([ones(size(Ax)) Ax], y);
+                    fit_params(1) = fit_params(1)*opts.WinsorThreshold;
+                    mappings{r} = fit_params;
                 else
-                    mappings{r} = [clock_values(idx(1)) 0]; % just one measurement
+                    mappings{r} = [offset_values(idx(1)) 0]; % just one measurement
                 end
             end
 
@@ -781,7 +788,7 @@ function close_file(f,filename)
 %     minimize 1/2*sum(huber(A*x - y))
 %
 % Based on the ADMM Matlab codes also found at:
-%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
+%   https://web.stanford.edu/~boyd/papers/admm_distr_stats.html
 %
 %                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
 %                                2013-03-04
@@ -790,15 +797,21 @@ function close_file(f,filename)
     rho = 1; end
 if ~exist('iters','var')
     iters = 1000; end
+
+x_offset = min(A(:, 2));
+A(:, 2) = A(:, 2) - x_offset;
 Aty = A'*y;
-L = sparse(chol(A'*A,'lower')); U = L';
+L = chol( A'*A, 'lower' );
+L = sparse(L);
+U = sparse(L');
 z = zeros(size(y)); u = z;
 for k = 1:iters
     x = U \ (L \ (Aty + A'*(z - u)));
     d = A*x - y + u;
     z = rho/(1+rho)*d + 1/(1+rho)*max(0,(1-(1+1/rho)./abs(d))).*d;
     u = d - z;
 end
+x(1) = x(1) - x(2)*x_offset;
 end