function Y=bridge(X,shift)
% BRIDGE	Missing-data interpolation
% Fills in any holes (NaN)
%
% Y = bridge(X,shift)
%
% X	= input vector, or array with time along first dimension (columns) 
% shift	= how far away to collect basis-data for mean bridge (default is
%         lenght of gap). Recommended when prominent and stationary
%         periodicity is known (like seasonality).
%
% To fill missing timepoints when Y is inhomogeneously sampled, put Y onto
% homogeneous grid first, to make the gaps show as NaNs. 

%Time-stamp:<Last updated on 06/09/14 at 11:30:22 by even@nersc.no>
%File:</home/even/matlab/evenmat/bridge.m>

if nargin<2|isempty(shift), shift=[];	end

D=size(X);
if D(2)>1 & D(1)==1, X=X(:); D=fliplr(D); end

if length(D)>2 | any(D(3:end)>1)	% if array
  cols=prod(D(2:end));
  Y=reshape(X,D(1),cols);
  Y=matbridge(Y',shift)';
  Y=reshape(Y,D);
else
  Y=matbridge(X',shift)';
end

% Plotting:
if nargout==0 & length(D)==2 & D(2)==1
  fig bridge 4;clf;
  h=plot(1:D(1),Y,'k.--',1:D(1),X,'k-o');
  legend(...%h([0:2]*N+1),...
	 'interpolated series','data input',0);
  grid on; zoom xon;
  scanplot;
  title('Interpolation by BRIDGE');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% KERNEL %%%%%%%%%%%%%%
function Y=matbridge(Y,shift)	% PS! Works along ROWS!
ydr=isvec(Y);
if ydr==1, Y=Y'; end			% flip col-vector into row

cpx=~isreal(Y);
if cpx, Y=[real(Y);imag(Y)]; end	% fill components separately

[M,N]=size(Y);

for i=1:M				% LOOP THE ROWS
  jn=find(isnan(Y(i,:)));		% indices of nans
  if isempty(jn), break; end
  jh=[find(diff(jn)>1) length(jn)];	% indices of holes in jn (endpoints)
  for j=1:length(jh)			% LOOP THE HOLES
    if j==1				% first hole
      jj=(jn(1)-1):(jn(jh(1))+1);	% indices of hole in Y w/bridgeheads
    else				% the rest
      jj=(jn(jh(j-1)+1)-1):(jn(jh(j))+1);
    end
    if isempty(shift),	sh=length(jj)-1;% automatic shift is bridgelength-1
    else		sh=shift;
    end
    %
    if any(jj>N)				% Gap at end of Y
      Y(i,jj(2:end-1))=Y(i,jj(2:end-1)-sh);	% copy series before
    elseif any(jj<1)				% Gap at start of Y
      Y(i,jj(2:end-1))=Y(i,jj(2:end-1)+sh);	% copy series after
    else					% Find mean-bridge:	
      left=jj-sh; right=jj+sh;
      if any(right>N)				% Not enough data after gap
	%[Y(i,jj(end):end) nans(1,sh)];		% add some nans
	%ym=nanmean([Y(i,left);ans(1:1+sh)]);	% nanmean uses data before
	[Y(i,:) nans(1,sh)];			% add enough nans
	ym=nanmean([Y(i,left);ans(right)]);	% nanmean uses data before
      elseif any(left<1)			% Not enough data before gap
	%[nans(1,sh) Y(i,1:jj(1))];		% ...
	%ym=nanmean([ans(end-sh:end);Y(i,right)]);	
	[nans(1,sh) Y(i,:)];			% ...
	ym=nanmean([ans(sh+1+left);Y(i,right)]);	
      else					% Away from ends of Y = OK
	ym=nanmean([Y(i,left);Y(i,right)]);	
      end
      n=length(jj);
      yml=interp1q([1;n],[ym(1);ym(n)],[1:n]')';   % mean-bridge's line
      yy=ym-yml;				   % deviation from it
      jbh=[jj(1) jj(end)];			   % indices of bridgeheads
      Yb=yy(2:end-1)+interp1q(jbh',Y(i,jbh)',jj(2:end-1)')';% Y-bridge =
						   % hole's line + deviation
      %
      neig=max(left(1)-4*sh,1):min(right(end)+4*sh,N);% Find neighbourhood ...
      2*nanstd(Y(i,neig));			      % range of ...
      yr=nanmean(Y(i,neig))+[-ans ans];		      % variability*
%      if any(find(Yb<yr(1)|Yb>yr(2))), keyboard;end %%%%TEST
      Yb(find(Yb<yr(1)))=yr(1);			% Restrain any ...
      Yb(find(Yb>yr(2)))=yr(2);			% generated spikes*
      Y(i,jj(2:end-1))=Yb;
    end
  end
end

if cpx, Y=complex(Y(1:M/2,:),Y((M/2+1):end,:));  end % into complex again

if ydr==1, Y=Y'; end		% flipback original vector-direction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ^^^KERNEL^^^ %%%%%%%%%%%


% test while in keyboard:
% global t
% cage([t(neig(1)) t(neig(end))],[yr(1) yr(2)]);