function [PatientData,PRxEventTimes,PTM]=CHARISGUIReportPRxEvents(patientID,thresholdPressure,minEventTime,minPreEventTime,numBars,rawICP, rawTime, PRx,PacketAvg_Time)

% The input patientID must be a string
% Example: '101314101514' will be output as ReportPRxEvent 101314101514.mat 

% numBars is in 5 second segments. 1 minute = 12 numBars


%Run PRx functions with inputs
[PatientData] = icpEventFinder(patientID,rawICP,rawTime,thresholdPressure,minEventTime,minPreEventTime);
eventTimes=PatientData.eventTimes';
[ PRxEventTimes , PTM] = eventIcp2eventPrx(eventTimes,rawTime,PacketAvg_Time,PRx);

end

function [PatientData] = icpEventFinder(patientID,rawICP,rawTime,thresholdPressure,minEventTime,minPreEventTime)

%make sure icp data is in mmHg before input
%rawTime must be in seconds
%thresholdPressure must be in mmHg
%minEventTime and minPreEventTime must be in minutes
%eventPreTime and eventTimeLengths are in minutes
%outTime is in seconds

if(length(rawICP)==length(rawTime))
    %average every minute, assumes 50 Hz
    [meanM,varM,x,~] = statEveryK(rawICP,50*60);
    
    %adjust the mean
    lowVarMeanX = find(varM<50);
    lowVarMean = meanM(lowVarMeanX);
    adjustedMean = interp1(lowVarMeanX,lowVarMean,1:length(meanM));
    
    %find possible event start locations
    eventThresholdPressure = thresholdPressure;
    highX = find(adjustedMean>=eventThresholdPressure);
    DhighX = diff(highX);
    highXstarters = [0,find(DhighX > 1)]+1;
    highXstartersCount = zeros(1,length(highXstarters));
    for starter = 1:length(highXstarters)-1
        highXstartersCount(starter) = length(find(DhighX(highXstarters(starter):highXstarters(starter+1)-1)==1));
    end
    highXstartersCount(end) = length(find(DhighX(highXstarters(end):end)==1));
    highXlocations = highX(highXstarters);
    
    %Choose the possible events that last for n minutes or more
    nMinutes = minEventTime;
    eventLocations = highXlocations(find(highXstartersCount>=nMinutes));
    eventTimeLengths = highXstartersCount(find(highXstartersCount>=nMinutes));
    
    %Choose the events that proceed a segment of regular icp
    regTimeThresh = minPreEventTime;
    regLowThresh = 0;
    eventRegTimeLengths = zeros(1,length(eventLocations));
    tempSeg = adjustedMean(1:eventLocations(1));
    eventRegTimeLengths(1) = backBandCounter(tempSeg(1:end-1),regLowThresh,eventThresholdPressure);
    for event = 2:length(eventLocations)
        tempSeg = adjustedMean(eventLocations(event-1):eventLocations(event));
        eventRegTimeLengths(event) = backBandCounter(tempSeg(1:end-1),regLowThresh,eventThresholdPressure);
    end
    eventsWithReg = find(eventRegTimeLengths>=regTimeThresh);
    eventLocations = eventLocations(eventsWithReg);
    eventTimeLengths = eventTimeLengths(eventsWithReg);
    
    outIndex = x(eventLocations);
    outTime = rawTime(x(eventLocations));
    eventPreTime = eventRegTimeLengths(eventsWithReg);
    
    %Go through each event
    goodEvents = ones(length(eventLocations),1);
    set(0,'DefaultFigureWindowStyle','docked');
    figure;
    for event = 1:length(eventLocations)
        clc;
        domain = x(eventLocations(event)):x(eventLocations(event)+eventTimeLengths(event));
        plot(rawTime(domain),rawICP(domain));
        hold on;
        plot(rawTime(x),adjustedMean,'g');
        axis([rawTime(x(eventLocations(event)-eventPreTime(event))) rawTime(x(eventLocations(event)+eventTimeLengths(event))) 0 50])
        xlabel('Time in seconds');
        ylabel('mmHg');
        title(sprintf('Event %d of %d.\n',event,length(eventLocations)));
        zoom xon;
        pan xon;
        hold off;
       
choice=questdlg('Is this a non-artifact Event?','Find Events','Yes','No','Yes');
switch choice
    case 'Yes'
        answer=1;
    case 'No'
        answer=2;
    case 'Cancel'
        
end
         if(answer == 2)
            goodEvents(event)=0;
        end
    end
    
    outIndex = outIndex(find(goodEvents==1));
    outTime = outTime(find(goodEvents==1));
    eventPreTime = eventPreTime(find(goodEvents==1));
    eventTimeLengths = eventTimeLengths(find(goodEvents==1));
    close
    
    %create eventResult PatientDataect
    crThreshold = sprintf('%d mmHg Threshold',thresholdPressure);
    crPreTime = sprintf('%d minutes pre-time',minPreEventTime);
    crPostTime = sprintf('%d minutes post-time',minEventTime);
    criteria = struct('threshold',crThreshold,'preTime',crPreTime,'postTime',crPostTime);
    atPreTimes = eventPreTime;
    atPostTimes = eventTimeLengths;
    eventAttributes = struct('postTimes',atPostTimes,'preTimes',atPreTimes);
    PatientData = CHARISicpEventResult(patientID,criteria,outTime,eventAttributes);
    
  
else
    fprintf('The ICP and time data are of different lengths.\n');
end
end

function [meanM,varM,indexStatStart,indexStatMid] = statEveryK(matrix,k)

[R,C] = size(matrix);
numSeg = floor(R./k);
if(numSeg<1)
    numSeg = 1;
end
meanM = zeros(numSeg+1,C);
varM = zeros(numSeg+1,C);
x = 1:numSeg+1;
for c = 1:C
    for s = 1:numSeg
        fprintf('   Segment %d of %d \n',s,numSeg);
        x(s) = (s.*k)-k+1;
        tempData = matrix((s.*k)-k+1:s.*k,c);
        meanM(s,c) = mean(tempData(find(~isnan(tempData))));
        varM(s,c) = var(tempData(find(~isnan(tempData))));
    end
    if((numSeg)+1<=R)
        x(numSeg+1) = numSeg.*k;
        tempData = matrix((numSeg.*k)+1:end,c);
        meanM(numSeg+1,c) = mean(tempData(find(~isnan(tempData))));
        varM(numSeg+1,c) = var(tempData(find(~isnan(tempData))));
    end
end
indexStatStart = x;
indexStatMid = x+(round(k/2));
end

function [ count ] = backBandCounter(vector,low,high)
%Counts the number of values that fall within the band from last to first.
%Count stops when a value falls outside the band.
flippedVector = flip(vector);
count = 0;
for k = 1:length(vector);
    if(flippedVector(k)>= low && flippedVector(k)<=high)
        count = count + 1;
    else
        break;
    end
end
end

function [location]=closestPoint(Vector,GivenPoint)

% Check unique Vector 
if length(unique(Vector))==length(Vector);
    

distance=nan(length(Vector),1);

% Get distances
for i=1:length(Vector)
    
    distance(i,1)=abs(Vector(i)-GivenPoint);
end
[~,location]=min(distance);

else
    fprintf('ERROR: Values are not Unique')
    location=[];
    
end
end

function [ PRxEventTimes , PTM] = eventIcp2eventPrx(eventTimes,rawTime,PRxTime,PRx)

PTM = getPRxIcpTimeMatrix(PRx,PRxTime);
PRxEventTimes = nan(length(eventTimes),3);
for c = 1:3;
    vector = PTM(:,c);
    for event = 1:length(eventTimes);
        time = eventTimes(event);
        location = closestPoint(vector,time);
        if(location == length(find(~isnan(vector))))
            PRxEventTimes(event,c) = nan;
        else
        PRxEventTimes(event,c) = location;
        end
    end
end
end

function [PTM] = getPRxIcpTimeMatrix(PRx,time)

[R,C] = size(PRx);
PTM = nan(R,C);
time = correctPRxTime(time);
min5Start = length(time)-length(PRx)+1;
min10Start = min5Start+60;
min20Start = min5Start+180;
PTM(:,1) = time(min5Start:end);
PTM(1:end-60,2) = time(min10Start:end);
PTM(1:end-180,3) = time(min20Start:end);
end

function [outTime] = correctPRxTime(inTime)


dTime = diff(inTime);
delta = mode(dTime);
X = find(dTime == delta);
V = inTime(X);
Xq = 1:length(inTime);
outTime = interp1(X,V,Xq)';
if(isnan(outTime(end)))
    outTime(end) = outTime(end-1)+5;
end

end