function y = evco(fname,estID,winsize,clbnum)
% EVCO  Evaluate CO estimation
%   Y = EVCO(FNAME, ESTID, WINSIZE, CLBNUM) calibrates and evaluates CO
%   estimation performance, using thermodilution CO as gold-standard.
%
%   In:   FNAME     (string)   file name -- e.g. fname='~/170.mat';
%         ESTID     (integer)  choose CO estimation algorithm
%         WINSIZE   (integer)  ABP window size (in seconds) to use
%         CLBNUM    (integer)  choose calibration method:
%                              = 1  single k MMSE
%                              = 2  multi  k causal MMSE
%                              = 3  single k 1st point
%                              = 4  multi  k previous point
%                              = 5  not calib, but normalize near 1
%                              = other --- no calib
%
%   Out:  Y
%         Col 1:  enumerate CO number
%             2:  time                      [min]
%             3:  TCO                       [L/min]
%             4:  estimated CO              [L/min] - possibly uncalibrated
%             5:  sigma of estimated CO     [L/min] - possibly uncalibrated
%             6:  sqi percent of good beats [mmHg]
%             7:  Systolic BP       [mmHg]
%             8:  Diastolic BP      [mmHg]
%             9:  Pulse pressure    [mmHg]
%            10:  Mean pressure     [mmHg]
%            11:  heart rate        [bpm]
%            12:  noisiness (sqi thresh at -3)
%            13:  Area under systole   0.3*sqrt(RR)  method
%            14:  Area under systole   1st min-slope method
% 
%   Example:
%   y = evco('170',5,60,1);
%
%   Written by James Sun (xinsun@mit.edu) on Jan 18, 2006.
%   - v2.0 - 1/30/06 - add features (cols 7 thru 14) to output

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get TCO, estimated CO data
load(fname,'t0','m2t');

% obtain TCO, sync TCO time with ABP time
t0_co = 24*60*(m2t.t0 - t0(1,1)); % minutes offset
r     = [(t0_co+m2t.CO(:,1))  m2t.CO(:,2)];

% run bb CO estimate
[est_val,t_est,t_raw,F] = estimateCO(fname,estID,0);
x     = [t_est*60 double(est_val)];
t_raw = t_raw*60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% windowing --- obtain one-to-one estimates with tco
mu    = zeros(size(r,1),1);
sigma = zeros(size(r,1),1);
sqi   = zeros(size(r,1),1);
fea   = zeros(size(r,1),8);
for i=1:size(r,1)
    T  = r(i,1)*60;

    ind = find(x(:,1)<T & x(:,1)>T-winsize);
    ind_raw = find(t_raw<T & t_raw>T-winsize);
    if isempty(ind)
        mu(i)    = NaN;
        sigma(i) = NaN;
        sqi(i)   = 0;
        fea(i,:) = NaN;
    else
        mu(i)    = mean(x(ind,2));
        sigma(i) = std(x(ind,2));
        sqi(i)   = round(length(ind)/length(ind_raw) * 100);
        fea(i,:) = mean(F(ind,[2 4 5 6 7 8 10 12]));
    end
end

% heart rate
fea(:,5) = 60*125./fea(:,5); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calibration
ind_nan = ~isnan(mu);

k = calib(r(ind_nan,2),mu(ind_nan),clbnum);
y = [(1:size(r,1))' r k.*mu k.*sigma sqi fea];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% format output
%y(isnan(mu),:)=[];