%n - number of points to be predicted
%plotoption allows to plot every step of the way with boundary
%if not plotoption - just plotting data and outliers
%index of missing points

function [outliers,outdates,sigmas] = spikedetsystem_mwar(fid,data,initpoints,pdlist,n,rank,D,missindex,DCTthresh)

%fid=fopen(filename,'w+');                   
ld = length(data);                          %length of data
%Preprocessing
index_zeros = find(data == 0);              %list of indexes of data points = 0
if ~isempty(index_zeros)
    data = interp_outliers(data, index_zeros);  %interp 0s in the dataset
end;
all_outliers_list = [];                     %list of indexes of all outliers detected
sigmas = [];                                %list of how many variances the outliers were from the predicted value 

i = initpoints;%ceil(ld*.5);

while i<=ld-n                  %run the detector starting with approximately half the points
  if ~isempty(missindex) & i == missindex(mi)     
        mi = min(length(missindex(:,1)),mi + 1);
        i = i+missindex(mi,ah)-missindex(mi,1);
  else
      boundary = normalstdev(pdlist,'n');
      model = newmodel(data,i,DCTthresh);%DCTthresholdSelection(data,i-n,n,'n');
      P = PredictionwithWavelets(model,n,rank); 
      pdlist = [pdlist;data(i+1:i+n) P(i+1:i+n)];
      sigmas = [sigmas;(i+1:i+n) (3/boundary)*(data(i+1:i+n)-P(i+1:i+n))];
      ind_outliers = find(data(i+1:i+n)>(P(i+1:i+n)+ones(n,1)*boundary)); %index of all outliers for the given n,model
      ind_outliers = ind_outliers + i*ones(size(ind_outliers));
    if ~isempty(ind_outliers)                 %if there are outliers      
      all_outliers_list = [all_outliers_list;ind_outliers];
    end;
    i = i+n;
  end;
end;
%fclose(fid);
outdates = D(all_outliers_list,1:8);
outliers = all_outliers_list;



%used to interpolate over 0s
function out = interp_outliers(data, index)
out = data;
temp = ones(length(index),1);
index+temp
index-temp
out(index) = (data(index+temp)+data(index-temp))./2;
