function [base,loss_ts,rt_loss_ts,consecloss_ts, consecarr_ts,delay_ts] = dag_ts(depfilename,arrfilename,plotflag)

% [base,loss_ts,rt_loss_ts,consecloss_ts, consecarr_ts,delay_ts] = dag_ts(depfilename,arrfilename, plotflag)
% 
% input dat file format:
%
% column 1 - timestamp, fraction of a second in 2^-32 sec units
% column 2 - timestamp, seconds
% column 3 - seq.number of the packet
%
% outputs:
%
% base:
% column 1,2 - departure timestamp
% column 3,4 - arrival timestamp
% missing packets = -1
% timestamp format: fractional part of the timestamp in 2^-32 sec (columns 1,3)
%                   whole sec part of the timestmap (columns 2,4)
%
% loss_ts: 
% loss time series
% column 1 - dep time in seconds relative to the first packet
% column 2 - loss indicator ( 1 = lost packet )
%
% rt_loss_ts: 
% real time loss series
% delta - size of the discrete timeframe (100*mean_inter_departure_time)
% column 1 - number of loss events in the given discrete timeframe
%
% consecloss_ts: 
% time series of consecutive losses
% column 1 - dep time of the first lost packet in the loss event 
%            in seconds relative to the first packet
% column 2 - number of packets lost
%
% consecarr_ts: 
% time series of consecutive arrivals
% column 1 - dep time of the first lost packet in the arrival event 
%            in seconds relative to the first packet
% column 2 - number of packets
%
% delay_ts: 
% delay time series
% column 1 - dep time of the packet in seconds relative to the first packet
% column 2 - delay in seconds (-1 indicates lost packet)
%

if nargin == 2
   plotflag = 1;
end

load (depfilename)
load (arrfilename)

% test for out of order packets

ooflag_dep =  ~isempty(find(diff(dep(:,3))<0))
ooflag_arr =  ~isempty(find(diff(arr(:,3))<0))

% test for duplicated packets
% delete the duplicates from the data

dpflag_dep =  ~isempty(find(diff(sort(dep(:,3))) == 0))
dpflag_arr =  ~isempty(find(diff(sort(arr(:,3))) == 0))

if (dpflag_dep) 
    dep = sortrows(dep,3);
    dep(find(diff(dep(:,3)) == 0) +1,3) = -1;
    dep = dep(find(dep(:,3) ~= -1),:);
end

if (dpflag_arr) 
    arr = sortrows(arr,3);
    arr(find(diff(arr(:,3)) == 0) +1,3) = -1;
    arr = arr(find(arr(:,3) ~= -1),:);
end

% construct the base data structure

base = zeros(max(dep(:,3))-min(dep(:,3))+1,4) - 1;
base(dep(:,3)-min(dep(:,3))+1,1:2) = dep(:,1:2);
base(arr(:,3)-min(dep(:,3))+1,3:4) = arr(:,1:2);


% loss time series
% column 1 - dep time in seconds relative to the first packet
% column 2 - loss indicator ( 1 = lost packet )

tmp = find(base(:,2) ~= -1);
loss_ts = zeros(size(tmp,1),2);
loss_ts(:,1) = base(tmp,2) - base(tmp(1),2) + (base(tmp,1) - base(tmp(1),1)).*(2^-32);
tmp = find(base(tmp,4) == -1);
loss_ts(tmp,2) = 1;
clear tmp;

% real time loss series
% delta - size of the discrete timeframe
% column 1 - number of loss events in the given discrete timeframe

delta = 100 * mean(diff(loss_ts(2:end,1)));

last_tptr = 0;
tptr = delta;
j = 1;

rt_loss_ts = zeros(floor(loss_ts(end,1)/delta) + 1,1);

while (tptr < loss_ts(end,1))
   rt_loss_ts(j) = sum(loss_ts(find(loss_ts(:,1) < tptr & loss_ts(:,1) >= last_tptr),2));
   j = j+1;
   last_tptr = tptr;
   tptr = tptr + delta;
end

clear tptr last_tptr j;

% time series of consecutive losses
% column 1 - dep time of the first lost packet in the loss event 
%            in seconds relative to the first packet
% column 2 - number of packets lost

lossdiff = diff(loss_ts(:,2));
loss_start = find(lossdiff == 1) + 1;
loss_end = find(lossdiff == -1) + 1;

if loss_ts(1,2) == 1
   loss_start = [1; loss_start];
end
if loss_ts(end,2) == 1
   loss_end = [loss_end; size(lossdiff,1)+2];
end

consecloss_ts = zeros(size(loss_start,1),2);
consecloss_ts(:,1) = loss_ts(loss_start,1);
consecloss_ts(:,2) = loss_end - loss_start;

clear loss_start loss_end lossdiff;

% time series of consecutive arrivals
% column 1 - dep time of the first lost packet in the arrival event 
%            in seconds relative to the first packet
% column 2 - number of packets

arrdiff = diff(~loss_ts(:,2));
arr_start = find(arrdiff == 1) + 1;
arr_end = find(arrdiff == -1) + 1;

if ~loss_ts(1,2) == 1
   arr_start = [1; arr_start];
end
if ~loss_ts(end,2) == 1
   arr_end = [arr_end; size(arrdiff,1)+2];
end

consecarr_ts = zeros(size(arr_start,1),2);
consecarr_ts(:,1) = loss_ts(arr_start,1);
consecarr_ts(:,2) = arr_end - arr_start;

clear arr_start arr_end arrdiff;

% time series of delay
% column 1 - dep time of the packet in seconds relative to the first packet
% column 2 - delay in seconds (-1 indicates lost packet)

delay_ts = loss_ts;
delay_ts(:,2) = -1;
tmp = find(base(:,2) ~= -1);
tmp = find(base(tmp,4) ~= -1);
delay_ts(tmp,2) = base(tmp,4) - base(tmp,2) + (base(tmp,3) - base(tmp,1)).*(2^-32);

clear tmp;

if plotflag == 1
   subplot(3,2,1);
   hold off
   plot(loss_ts(find(loss_ts(:,2) == 0),1),loss_ts(find(loss_ts(:,2) == 0),2),'*','MarkerSize',2);
   hold on
   plot(loss_ts(find(loss_ts(:,2) == 1),1),loss_ts(find(loss_ts(:,2) == 1),2),'r*','MarkerSize',2);
   hold off
   title('loss\_ts');
   ylabel('loss flag');
   subplot(3,2,2);
   plot(rt_loss_ts(:,1),'*','MarkerSize',2);
   title('rt\_loss\_ts in discrete time');
   ylabel('# of lost packets');
   subplot(3,2,3);
   plot(consecloss_ts(:,1),consecloss_ts(:,2),'r*','MarkerSize',2);
   title('consecloss\_ts');
   ylabel('# of lost packets');
   subplot(3,2,4);
   title('consecarr\_ts');
   plot(consecarr_ts(:,1),consecarr_ts(:,2),'g*','MarkerSize',2);
   title('consecarr\_ts');
   ylabel('# of packets');
   subplot(3,2,5);
   plot(delay_ts((find(delay_ts(:,2)>-1)),1),delay_ts((find(delay_ts(:,2)>-1)),2),'*','MarkerSize',2); 
   title('delay\_ts');
   ylabel('delay [sec]');
   xlabel('dep time [sec]');
   subplot(3,2,6);
   plot(loss_ts(:,1),[0;diff(loss_ts(:,1))],'*','MarkerSize',2);
   title('inter departure times');
   ylabel('inter departure time [sec]');
   xlabel('dep time [sec]');
end