function [co, to, told, fea] = estimateCO(fname,estID,filt_order)
% ESTIMATECO  Cardiac output estimator.
%   [CO, TO, TOLD, FEA] = ESTIMATECO(FNAME, ESTID, FILT_ORDER) is the 
%   main function for estimating cardiac output.  
%
%   In:   FNAME      <1xn>  --- file where ABP and features are located
%         ESTID      <1x1>  --- the CO estimator to use
%         FILT_ORDER <1x1>  --- order of running avg LPF to use on output 
%                                      (enter 0 or 1 for no LPF)
%
%   Out:  CO         <kx1>  --- estimated CO (uncalibrated)
%         TO         <kx1>  --- time [minutes] (not evenly sampled!)
%         TOLD       <mx1>  --- time [minutes], not sqi filtered
%         FEA        <kx11> --- feature matrix
% 
%   Usage:
%   This function is only a wrapper.  Actual CO estimation computation is 
%   done in the following required functions:
%   ESTID 1: est01_MAP      - Mean pressure
%         2: est02_WK       - Windkessel 1st order LTI RC circuit model
%         3: est03_SA       - Systolic area distributed model
%         4: est04_SAwarner - Warner systolic area with time correction
%         5: est05_Lilj     - Liljestrand PP/(Psys+Pdias) estimator
%         6: est06_Herd     - Herd estimator
%         7: est07_SAwessCI - Wesseling systolic area with impedance correction
%         8: est08_Pulsion  - Pulsion non-linear compliance model
%         9: est09_LidCO    - LidCO root-mean-square model
%        10: est10_RCdecay  - RC exponential decay fit
%        11: est11_mf       - Wesseling non-linear, time-varying 3-element model
%
%   Written by James Sun (xinsun@mit.edu) on Nov 19, 2005.
%   - v2.0 - 1/18/06 - added output for detecting percentage of good beats
%   - v3.0 - 1/20/06 - added "feature" output, CO no longer normalized

% load features from fname
load(fname,'F','onset','sqi','t0','m2t');
load agegender

co = [];
to = [];
told=[];
fea = [];
for seg=1:size(t0,1)  % loop through each ABP waveform segment
%seg
    onset1 = onset{seg};

    % skip segment if too little data
    if length(onset1)<50 || length(onset1)<(3*filt_order)
        continue
    end

    % read features
    F1      = F{seg};
    Psys    = F1(:,2);
    Pdias   = F1(:,4);
    PP      = F1(:,5);
    MAP     = F1(:,6);
    Period  = F1(:,7);
    HR      = 60*125./Period; 
    tSArr   = F1(:,9);
    SArr    = F1(:,10);
    tSAcosm = F1(:,11);
    SAcosm  = F1(:,12);

    % select systolic area estimate
    SA      = SArr;
    tSA     = tSArr;

    % special items for estID 8,9,11
    if estID==11
        caseid = m2t.caseid;
        ind = find(agegender(:,1)==caseid);
        % quit if no age, gender info found
        if isempty(ind)
            continue
        end
        age = agegender(ind,3);
        gender = agegender(ind,4);
    end
    
    if estID==8 || estID==9 || estID==11
        str = sprintf('load %s abp%d',fname,seg); eval(str);
        str = sprintf('abp=abp%d; clear abp%d',seg,seg); eval(str);
    end

    
    % estimate CO!
    switch estID
        case  1,   x = est01_MAP(MAP);
        case  2,   x = est02_WK(PP,HR);
        case  3,   x = est03_SA(SA,HR);
        case  4,   x = est04_SAwarner(SA,HR,onset1,tSA);
        case  5,   x = est05_Lilj(PP,HR,Psys,Pdias);
        case  6,   x = est06_Herd(MAP,Pdias,HR);
        case  7,   x = est07_SAwessCI(MAP,SA,HR);
        case  8,   x = est08_Pulsion(abp, onset1, MAP, HR, tSA);
        case  9,   x = est09_LidCO(abp,onset1,HR);
        case 10,   x = est10_RCdecay(Period,Psys,Pdias,MAP);
        case 11,   x = est11_mf(abp,onset1,MAP,HR,age,gender);
        case 12,   x = est12_coalees(PP,HR,onset1,tSA,Psys,Pdias);
        case 13,   x = est13_trivial(PP,HR,MAP);
        otherwise, x = nan;
    end
    
    % synchronize time of of all segments
    offset = (t0(seg,1)-t0(1,1))*24*60;
    t      = onset1(1:end-1)/(125*60)+offset; %[min]

    % remove all datapoints with bad SQI
    sqi1   = sqi{seg}(:,1);
    ind    = find(sqi1);
    t_old  = t;
    x(ind) = [];
    t(ind) = [];

    % apply zero-phase moving avg LPF
    if filt_order<2, x_filt = x;
    else             x_filt = filtfilt(ones(filt_order,1)/filt_order,1,x);
    end

    % features
    ff = F{seg};
    ff(ind,:) = [];
    
    % append segment data to pool
    co = [co; x_filt];
    to = [to; t];
    told = [told; t_old];
    fea = [fea; ff];
end

% normalize estimated CO to 1st point
%co = co/co(1);