function [output,errors] = ProcessMetabolics(paths,varargin)
    % PROCESSMETABOLICS   Process COSMED data from an excel file
    %   Written by Grant Elliott, June 2011
    %
    %   [Analysis, errors] = ProcessMetabolics(filenames) retrieves data from excel spreadsheets 
    %   produced by COSMED with paths specified by the strings stored in filenames and 
    %   calculates metabolic demands with error bars following user input.
    %   
    %   May not work if Excel is not installed.
    %
    %   filenames may be a string or cell array of strings containing paths to Excel spreadsheets.
    %
    %   Analysis is a structure array of results corresponding to those files which were processed without
    %   error, a list given by filenames{~errors}
    %
    %   error is a logical array, true wherever an error prevented property processing a datafile.
    %
    %   [Analysis, errors] = ProcessMetabolics(filenames, paramName, paramValue,...) takes additional optional
    %   parameters, which may be cell array of the same length as filenames
    %
    %   Accepts the following optional parameters:
    %       Bounds                   - 1x2 Array specifying time range to analyze (Defaults to prompt)
    %       PlotResults              - Leaves a figure of results on screen (Defaults to true)
    %       SubjectMass              - Subject mass in kg (Defaults to cell B7 +/- 1kg)
    %       SubjectHeight            - Subject height in cm (Defaults to cell B6 +/- 1cm)
    %       SubjectAge               - Subject age in years (Defaults to cell B5)
    %       SubjectSex               - Subject sex (M/F) (Defaults to cell B4)
    %       SubjectID                - Subject ID number (Defaults to cell B1)
    %       SubjectFirstName         - Subject first name (Defaults to cell B2)
    %       SubjectLastName          - Subject last name (Defaults to cell B3)
    %       Date                     - Date vector of trial time (Defaults to cells E2 and E3)
    %       Pressure                 - Barometric pressure in mmHG during trial (Defaults to cell H1)
    %       Temperature              - Temperature in C during trial (Defaults to cell H2)
    %       Humidity                 - Relative humidity during trial (Defaults to cell H3/100)
    %       VO2Weight                - O2 exchange contributon to power in W/L (Defaults to Brockway's 16580)
    %       VCO2Weight               - CO2 exchange contributon to power in W/L (Defaults to Brockway's 4510)
    %       StartLockout             - Time from start of data before data may not be used in analysis in s (Defaults to 20)
    %       StopLockout              - Time from end of data after which data may not be used in analysis in s (Defaults to 10)
    %       BreathHoldLockout        - Time after a breath hold during which data may not be used in analysis in s (Defaults to 10)
    %       NotesField               - Annotations Field name for notes in cell B8 (Defaults to not used)
    %       User1Field               - Annotations Field name for User1 data in cell H19 (Defaults to not used)
    %       User2Field               - Annotations Field name for User2 data in cell H20 (Defaults to not used)
    %       User3Field               - Annotations Field name for User3 data in cell H21 (Defaults to not used)
    %       Annotations              - Cell array of the form {'FieldName', value,...} to be stored in the final output
    %       AnnotationsFromSheet     - Cell array of the form {'FieldName', Rows, Cols,...} to be stored in the final output
    %                                  where Rows and Cols are indices or arrays of indices into the spreadsheet
    %
    %   Analysis contains the following fields:
    %       FileName                 - String containing the filename of the original excel spreadsheet
    %       Annotations              - Structure containing fields specified in Annotations, AnnotationsFromSheet, NotesField, User#Field
    %       Subject                  - Structure containing subject data
    %       TrialConditions          - Structure describing conditions during data acquisition
    %       AnalysisConditions       - Structure describing options used for the analysis
    %       Data                     - Structure containing raw data from the spreadsheet
    %       Window                   - Structure containing the user selected window
    %       Results                  - Structure containing intermediate and final results of metabolic calculations
    %
    %   Analysis.Subject contains the following fields
    %       Mass                     - Measurement of subject's mass in kg
    %       Height                   - Measurement of subject's height in cm
    %       Age                      - Number representing subject age in years
    %       Sex                      - String containing subject's sex (M/F)
    %       FirstName                - String containing subject's first name
    %       LastName                 - String containing subject's last name
    %
    %   Analysis.TrialConditions contains the following fields
    %       Date                     - Date vector representing trial time
    %       Pressure                 - Number representing atmospheric pressure in mmHg
    %       Temperature              - Number representing temperature in degrees Celcius
    %       Humidity                 - Number representing relative humidity (0 to 1)
    %
    %   Analysis.AnalysisConditions contains the following fields
    %       LockoutStart             - Number representing seconds from start before data may be used for analysis
    %       LockoutStop              - Number representing seconds from end after which data may not be used for analysis
    %       LockoutHold              - Number representing seconds from end of a breath hold before data may be used for analysis
    %       VO2Weight                - Number representing O2 exchange contribution to energy in J/L
    %       VCO2Weight               - Number representing CO2 exchange contribution to energy in J/L
    %
    %   Analysis.Data contains the following fields
    %       VO2                      - Array of O2 exchange readings in L/s
    %       VCO2                     - Array of CO2 exchange readings in L/s
    %       TimeStamped              - Array of original COSMED time stamps in s
    %       TimeResampled            - Array of resampled high precision time stamps in s
    %       BreathTotalDuration      - Array of breath durations in s
    %       BreathInDuration         - Array of inhalation times in s
    %       BreathOutDuration        - Array of exhalation times in s
    %       BreatchCount             - Number of breaths
    %       GlitchCount              - Number of sensor glitches
    %       HoldCount                - Number of breath holds
    %       GlitchIndices            - Array of indexes to sensor glitches
    %       HoldIndices              - Array of indexes to end of breath holds
    %
    %   Analysis.Window contains the following fields
    %       Bounds                   - Array containing lower and upper bounds of the window in s
    %       GlitchCount              - Number of sensor glitches excluded from the window
    %       HoldCount                - Number of breath holds excluded from the window
    %       WindowIndices            - Array of indexes to breaths used in the window
    %       BreathCount              - Number of breaths used in the window
    %       Duration                 - Number of seconds of breathing used in the window
    %
    %   Analysis.Results contains the following fields
    %       VO2                      - Measurement representing mean O2 exchange rate in L/s
    %       VCO2                     - Measurement representing mean CO2 exchange rate in L/s
    %       R                        - Measurement repreesnting mean ratio of VO2 and O2 exchange
    %       Power                    - Measurement representing total metabolic power in W
    %       NormalizedPower          - Measurement representing total metabolic power per mass in W/kg
    %
    %   See also PlotMetabolics, Measurement.
    
    [output, errors] = ProcessPaths(paths,{'xls','xlsx'},@(path)ProcessMetabolicsCallback(path,varargin{:}));
end

function [output,abort]=ProcessMetabolicsCallback(filename,varargin)
    abort=false;
    
    try
        %not using numdat and txtdat because of indexing issue if the stock excel file has been modified
        %This should work in basic mode on machines without excel installed, but in some cases it doesn't.  That's your problem...
        warning off MATLAB:xlsread:ActiveX;
        [~,~,rawdat] = xlsread(filename);
    catch c
        error('Metabolics:PlotMetabolics:BadFileName','Failed to open file as Excel data');
    end

    p=inputParser;
    p.addParamValue('Bounds',NaN,@(x) isnumeric(x) && all(x>0) && numel(x)==2);
    p.addParamValue('SubjectMass',Measurement(cell2mat(rawdat(7,2)),1));
    p.addParamValue('SubjectHeight',Measurement(cell2mat(rawdat(6,2)),1));
    p.addParamValue('SubjectAge',cell2mat(rawdat(5,2)), @(x) isnumeric(x) && numel(x)==1 && x>=0);
    p.addParamValue('SubjectFirstName',cell2mat(rawdat(3,2)), @(x) ischar(x));
    p.addParamValue('SubjectLastName',cell2mat(rawdat(2,2)), @(x) ischar(x));
    p.addParamValue('SubjectID',cell2mat(rawdat(1,2)), @(x) isnumeric(x) && numel(x)==1);
    p.addParamValue('SubjectSex',cell2mat(rawdat(4,2)), @(x) ischar(x) && numel(x)==1 && (upper(x)=='M' || upper(x)=='F'));
    p.addParamValue('Date',datevec([cell2mat(rawdat(2,5)) ' ' cell2mat(rawdat(3,5))]), @(x) isnumeric(x) && all(size(x)==[1,6]));
    p.addParamValue('Pressure',cell2mat(rawdat(1,8)), @(x) isnumeric(x) && numel(x)==1 && x>=0);
    p.addParamValue('Temperature',cell2mat(rawdat(2,8)), @(x) isnumeric(x) && numel(x)==1);
    p.addParamValue('Humidity',cell2mat(rawdat(3,8))/100, @(x) isnumeric(x) && numel(x)==1 && x>=0 && x<=1);
    p.addParamValue('VO2Weight',16580, @(x) isnumeric(x) && numel(x)==1);
    p.addParamValue('VCO2Weight',4510, @(x) isnumeric(x) && numel(x)==1);
    p.addParamValue('StartLockout',20, @(x) isnumeric(x) && numel(x)==1 && x>=0);
    p.addParamValue('StopLockout',10, @(x) isnumeric(x) && numel(x)==1 && x>=0);
    p.addParamValue('BreathHoldLockout',10, @(x) isnumeric(x) && numel(x)==1 && x>=0);
    p.addParamValue('PlotResults',true, @(x) islogical(x) && numel(x)==1);
    p.addParamValue('NotesField','',@(x) ischar(x) && ~isempty(x));
    p.addParamValue('User1Field','',@(x) ischar(x) && ~isempty(x));
    p.addParamValue('User2Field','',@(x) ischar(x) && ~isempty(x));
    p.addParamValue('User3Field','',@(x) ischar(x) && ~isempty(x));
    p.addParamValue('Annotations',{}, @(x) iscell(x) && mod(numel(x),2)==0 && all(ischar([x{1:2:end}])));
    p.addParamValue('AnnotationsFromSheet',{}, @(x) iscell(x) && mod(numel(x),3)==0 && all(ischar([x{1:3:end}])) && (length([x{2:3:end}])==1 || length([x{3:3:end}])==1 || length([x{2:3:end}])==length([x{3:3:end}])) && all(isnumeric([x{2:3:end}])) && all(isnumeric([x{3:3:end}])))
    p.parse(varargin{:});
    settings = p.Results;

    settings.SubjectSex = upper(settings.SubjectSex);
    settings.SubjectFirstName = upper(settings.SubjectFirstName);
    settings.SubjectLastName = upper(settings.SubjectLastName);
    settings.SubjectMass = Measurement(settings.SubjectMass);
    settings.SubjectHeight = Measurement(settings.SubjectHeight);
    
    VO2 = cell2mat(rawdat(4:end,14))'/(1000*60); %L/s
    VCO2 = cell2mat(rawdat(4:end,15))'/(1000*60); %L/s
    inhalationDurations = cell2mat(rawdat(4:end,35))'; %s
    exhalationDurations = cell2mat(rawdat(4:end,36))'; %s
    breathDurations = inhalationDurations+exhalationDurations;
    timeVec = datevec(cell2mat(rawdat(4:end,10)),'HH:MM:SS');
    stampTime = (timeVec(:,6)+60*timeVec(:,5)+60*60*timeVec(:,4))'; %s
    time = [0,cumsum(inhalationDurations(2:end)+exhalationDurations(1:end-1))];
    
    %Append AnnotationsFromSheet
    tempAnnotations = cell(1,2/3*numel(settings.AnnotationsFromSheet));
    for i=1:numel(settings.AnnotationsFromSheet)/3
        xcoord = settings.AnnotationsFromSheet{(i-1)*3+2};
        ycoord = settings.AnnotationsFromSheet{(i-1)*3+3};
        val = rawdat(xcoord,ycoord);
        if all(cellfun(@(x)isnumeric(x), val))
            val=cell2mat(val);
        end
        tempAnnotations{(i-1)*2+1}=settings.AnnotationsFromSheet{(i-1)*3+1};
        tempAnnotations{(i-1)*2+2}=val;
    end
    settings.Annotations=[settings.Annotations tempAnnotations];
    %Append Notes field
    if ~isempty(settings.NotesField)
        note = cell2mat(rawdat(8,2));
        if isnan(note)
            note='';
        end
        settings.Annotations=[settings.Annotations, {settings.NotesField, note}];
    end
    %Append User fields
    for i=1:3
        field = settings.(sprintf('User%dField',i));
        if ~isempty(field)
            note = cell2mat(rawdat(18+i,8));
            if isnan(note)
                note=0;
            end
            settings.Annotations=[settings.Annotations, {field, note}];
        end
    end
    
    count = length(stampTime);
    holds=[];
    sectionStart=1;
    while(sectionStart<count)
        time(sectionStart:end) = time(sectionStart:end) - (time(sectionStart)-stampTime(sectionStart));
        gap = find(stampTime-time > 4.9,1);
        if (isempty(gap))
            gap = count+1;
        else
            holds=[holds,gap]; %#ok<AGROW>
        end
        time(sectionStart:gap-1) = time(sectionStart:gap-1) - mean(time(sectionStart:gap-1)-stampTime(sectionStart:gap-1));
        sectionStart = gap;
    end
    glitches=find(or(VO2==0,VCO2==0));

    output.FileName = filename;
    output.Annotations = cell2struct(settings.Annotations(2:2:numel(settings.Annotations)),settings.Annotations(1:2:numel(settings.Annotations)),2);
    output.Subject.ID = settings.SubjectID;
    output.Subject.LastName = settings.SubjectLastName;
    output.Subject.FirstName = settings.SubjectFirstName;
    output.Subject.Sex = settings.SubjectSex;
    output.Subject.Age = settings.SubjectAge;
    output.Subject.Mass = settings.SubjectMass;
    output.Subject.Height = settings.SubjectHeight;
    output.TrialConditions.Date = settings.Date;
    output.TrialConditions.Pressure = settings.Pressure;
    output.TrialConditions.Temperature = settings.Temperature;
    output.TrialConditions.Humidity = settings.Humidity;
    output.AnalysisConditions.LockoutStart = settings.StartLockout;
    output.AnalysisConditions.LockoutStop = settings.StopLockout;
    output.AnalysisConditions.LockoutHold = settings.BreathHoldLockout;
    output.AnalysisConditions.VO2Weight = settings.VO2Weight;
    output.AnalysisConditions.VCO2Weight = settings.VCO2Weight;
    output.Data.VO2 = VO2;
    output.Data.VCO2 = VCO2;
    output.Data.TimeStamped = stampTime;
    output.Data.TimeResampled = time;
    output.Data.BreathTotalDuration = breathDurations;
    output.Data.BreathInDuration = inhalationDurations;
    output.Data.BreathOutDuration = exhalationDurations;
    output.Data.BreathCount = count;
    output.Data.GlitchCount = length(glitches);
    output.Data.HoldCount = length(holds);
    output.Data.GlitchIndices = glitches;
    output.Data.HoldIndices = holds;

    windowMask=ones(size(time));
    windowMask(time<time(1)+settings.StartLockout)=0;                                %Do not use records from the beginning of the trial
    windowMask(time>time(end)-settings.StopLockout)=0;                               %Do not use records from the end of the trial
    windowMask(glitches)=0;                                                          %Do not use samples that glitched
    windowMask(holds-1)=0;                                                           %Do not use samples immediately prior to breath holds
    windowMask(holds)=0;                                                             %Do not use samples immediately after breath holds
    for i=holds
        windowMask(and(time>=time(i),time<time(i)+settings.BreathHoldLockout))=0;    %Do not use samples following breath holds
    end
    
    if isnan(settings.Bounds)
        PlotMetabolics(output);
        lowerIndex=[];
        upperIndex=[];
        while(isempty(lowerIndex) || isempty(upperIndex) || lowerIndex>upperIndex)
            fprintf('Select lower and upper bounds\n');
            try
                %bounds = ginput(2);
                grid on;
                bounds = input('Type input time range: [min(time) max(time)]:\n');
            catch c %#ok<NASGU>
                abort = true;
                return;
            end
            %bounds = sort(bounds(:,1))';
            bounds = transpose(sort(bounds));
            lowerIndex=find(time>bounds(1),1);
            upperIndex=find(time<bounds(2),1,'last');
        end
    else
        sort(settings.Bounds);
        lowerIndex=find(time>settings.Bounds(1),1);
        upperIndex=find(time<settings.Bounds(2),1,'last');
    end
    bounds(1)=time(lowerIndex)-inhalationDurations(lowerIndex);
    bounds(2)=time(upperIndex)+exhalationDurations(upperIndex);
    fprintf('Using Bound: %0.3fs to %0.3fs\n',bounds(1),bounds(2));
    %If a nulled region overlaps either end of the window, redefine the effective bounds
    %(If the window didn't actually include a gap, there's no reason to warn the user)
    if (windowMask(lowerIndex)==0 || windowMask(upperIndex)==0)
        lowerIndex = lowerIndex + find(windowMask(lowerIndex:end),1)-1;
        upperIndex = find(windowMask(1:upperIndex),1,'last');
        if (isempty(lowerIndex) || isempty(upperIndex) || lowerIndex>upperIndex)
            fprintf('Failed to find usable data within bounds.  Not calculating metabolic power');
            return;
        end
        bounds(1)=time(lowerIndex)-inhalationDurations(lowerIndex);
        bounds(2)=time(upperIndex)+exhalationDurations(upperIndex);
        fprintf('Revising Bound to Legal Region: %0.3fs to %0.3fs\n',bounds(1),bounds(2));
    end
    windowMask(or(time<bounds(1),time>bounds(2)))=0;
    windowIndices = find(windowMask);
    windowCount = length(windowIndices);
    windowTime = sum(breathDurations(windowIndices));
    windowGlitches = sum(and(glitches>=lowerIndex,glitches<=upperIndex));
    windowHolds = sum(and(holds>=lowerIndex,holds<=upperIndex));
    if (windowGlitches==1)
        fprintf('Warning: 1 Glitch Excluded From Window\n')
    elseif (windowGlitches~=0)
        fprintf('Warning: %d Glitches Excluded From Window\n',windowGlitches);
    end
    if (windowHolds==1)
        fprintf('Warning: 1 Breath Hold Excluded From Window\n')
    elseif (windowHolds~=0)
        fprintf('Warning: %d Breath Holds Excluded From Window\n',windowHolds);
    end

    VO2Pts = VO2(windowIndices);
    VCO2Pts = VCO2(windowIndices);
    powerPts = VO2Pts*settings.VO2Weight + VCO2Pts*settings.VCO2Weight;
    RPts = VCO2Pts./VO2Pts;
    weightPts = breathDurations(windowIndices);
    [~,VO2Est] = Measurement.meanUncorrelatedWeighted(VO2Pts, weightPts);
    [~,VCO2Est] = Measurement.meanUncorrelatedWeighted(VCO2Pts, weightPts);
    [~,powerEst] = Measurement.meanUncorrelatedWeighted(powerPts, weightPts);
    [~,REst] = Measurement.meanUncorrelatedWeighted(RPts, weightPts);
    normalizedPowerEst = powerEst/output.Subject.Mass;

    pm=char(177);
    fprintf('Calculated metabolic demand using %d breaths totaling %0.3fs\n',windowCount,windowTime);
    fprintf('Mean VO2:  %0.2f %c %0.2f mL/s\n',VO2Est.Value*1000,pm,VO2Est.Uncertainty*1000);
    fprintf('Mean VCO2: %0.2f %c %0.2f mL/s\n',VCO2Est.Value*1000,pm,VCO2Est.Uncertainty*1000);
    fprintf('Mean ReR: %0.2f %c %0.2f\n',REst.Value,pm,REst.Uncertainty);
    fprintf('Power: %0.1f %c %0.1f W\n',powerEst.Value,pm,powerEst.Uncertainty);
    fprintf('Normalized Power: %0.2f %c %0.2f W/kg\n',normalizedPowerEst.Value,pm,normalizedPowerEst.Uncertainty);
    
    output.Window.Bounds = bounds;
    output.Window.GlitchCount = windowGlitches;
    output.Window.HoldCount = windowHolds;
    output.Window.WindowIndices = windowIndices;
    output.Window.BreathCount = windowCount;
    output.Window.Duration = windowTime;
    output.Results.VO2 = VO2Est;
    output.Results.VCO2 = VCO2Est;
    output.Results.R = REst;
    output.Results.Power = powerEst;
    output.Results.NormalizedPower = normalizedPowerEst;
    
    if (settings.PlotResults)
        PlotMetabolics(output,'axes',gca);
    else
        close(gcf);
    end
end