function [VrEst, dVuEst, VEEst] = EstResidVols_v1p2(Vv,Vi,dT,VrN,Vr0,EstMeth)
% Estimate residual volumes in serial cystometry from file or arrays.
%
%
% %% SYNTAX 1: obtain data from Excel file %%
% [EstVr, EstdVu, VEEst] = EstResidVols_v1p2(ExcelFile,worksheet,EstMeth)
%
% -- INPUTS --
% ExcelFile:        File name (or complete path) to the Excel file
%                   containing the cystometry data [string]. Use [] as this
%                   argument to initiate manual file selection. Note that 
%                   if the file cannot be loaded the function will prompt
%                   the user to select the file via GUI.
% worksheet:        Name of the worksheet in ExcelFile to use [string]
% EstMeth:          Name of the method to use for estimating the post-void
%                   residual volumes, options described below [string]
%
% = Format of Excel Data File =
%       The first row of the data file must be a list of headers for each
%       column. Each row must be data from the cystometry, ordered
%       chronologically. The data file must contain the following headers
%       with the corrisponding data in that column in any order:
%       "Vol Voided":       column of measured voided volumes
%       "Vol Infused":      column of experimenter infused volume
%       "Intervoid Time":   column of time between voids
%       "Last Resid":       one final measured residual volume
%       "Initial Resid":    one initial starting volume of the bladder
%
%       These columns can be in any order (although rows must be
%       chronological), any text can appear in the heading name after the
%       required text, and it is permissible to include other columns not
%       listed here (they will be ignored). See provided data file for an
%       example.
%
% %% SYNTAX 2: directly input arrays %%
% [EstVr, EstdVu, VEEst] = EstResidVols_v1p2(Vv,Vi,dT,VrN,Vr0,EstMeth)
%
% -- INPUTS --
% Vv:               Measured voided volumes [1D array]
% dT:               Time between void (or measurement) events [1D array]
% Vi:               Infused bladder volume (if any)
% VrN:              Final measured residual volume [scalar]
% Vr0:              Starting bladder volume [scalar]
% EstMeth:          Name of the method to use for estimating the post-void
%                   residual volumes, options described below [string]
%
%
% -- OUTPUTS --
% VrEst:            Estimated residual volumes
% dVuEst:           Estimated ureter flow rates
% VEEst:            Estimated voiding efficiency
%
%
%
% Estimation Method Descriptions: Algorithm to use for estimation
%       1) 'ConstantFlow'
%               Solves the conservation of mass equation for the unique
%               solution that results in the final Vr being equal to the
%               measured value, VrN, assuming the ureter flow rate is equal
%               between all voiding events.
%               This is the method described in: Danziger, ZC, Grill, WM.
%               Estimating postvoid residual volume without measuring
%               residual bladder volume during serial cystometrograms.
%               Am J Physiol Renal Physiol. 2016;311 (2):F459-68.
%               doi: 10.1152/ajprenal.00516.2015
%               
%       2) 'MinFlowVariance'
%               Estimates (unequal) ureter flow rates by minimizing the
%               prediction error of the final measured Vr, subject to the
%               constraints that 1) every ureter flow is positive, 2) all
%               estimated Vr are non-negative, 3) variability of ureter
%               flow rates between voids should be as low as possible.
%               This is the method described in: Angoli, D, Geramipour, A,
%               Danziger, ZC. Validation of an efficient and continuous
%               urodynamic monitoring system for awake, unrestrained,
%               chronic rodent studies. Am J Physiol Renal Physiol. 2020;
%               318 (1):F86-F95. doi: 10.1152/ajprenal.00349.2019
%
%
%
% EXAMPLES
% Estimate residuals from the 1st sheet in the example data file.
%
% fileName = 'Data1-AcuteWithReplacement.xlsx';  % file must be on your MATLAB path
% [VrEst, dVuEst] = EstResidVols(fileName,1,'MinFlowVariance')
%
% 
% Estimate residuals using syntax 2, directly inputing arrays. For the
% volume infused to the bladder (Vi), intervoid interval (dT), measured
% voided volumes (Vv), the final measured residual volume of 0.15, and the
% bladder starting volume for this batch of 0.
%
% Vi = [0.648 0.289 0.507 0.542 0.424 0.419 0.569 0.617 0.547];
% dT = [623.54 410.9 632.7 661.8 547.9 511.45 671.05 719.3 672.9];
% Vv = [0.545 0.39 0.5 0.53 0.42 0.43 0.61 0.555 0.475];
% EstVr = EstResidVols_v1p2(Vv,Vi,dT,0.15,0,'MinFlowVariance')
%
%
% NOTE:
%   The user is responsible for ensuring consistent units! All times and
%   volumes must have consistent units.
%
%
%
% %%% ZC Danziger Nov 2018 %%%
%
%
% EDITS ZC Danziger Dec 2021: a) Updated documentation. b) Updated input
% handling to minimize user error. c) Changed computation method of
% constant flow rate algorithm to the analytical solution for speed. d)
% Changed non-constant flowrate initial optimization guess to conform to
% constraints to improve optimizaiton stability. e) Added check on data
% fidelity to warn user if there is too little measured voided and residual
% volumes to make physical sense.
%




%% Handle Inputs

% determine input mode
inputMode = isnumeric(Vv) && ~isempty(Vv);

if inputMode
    % Case 2 where user is supplying an array:
    
    % input checks
    svv = size(Vv);
    if svv(1)~=1, error('Voided Volumes must be of size 1 x NumOfVoids'); end
    if ~all(svv==size(dT)), error('Voids and Time array must have same number of elements'); end
    if ~all(svv==size(Vi)), error('Voids and Infusion array must have same number of elements'); end
    if any(svv<0) || any(dT<0) || any(Vi<0), error('All voided volumes and times must be positive'); end
    if ~all(size(VrN)==[1 1]) || ~all(size(Vr0)==[1 1]) || VrN<0 || Vr0<0
        error('Initial bladder volume and final residual measurement must be non-negative scalars')
    end
    
else
    % Case 1 where user is supplying an Excel file:
    fn = Vv;
    sn = Vi;
    EstMeth = dT;
    try
        % load file
        [~,~,rawData] = xlsread(fn,sn);
    catch
        % prompt user for file - loading failed
        [fn,pn] = uigetfile({'*.xls;*.xlsx;*.cvs'});
        [~,~,rawData] = xlsread([pn '\' fn],sn);
    end
    
    % locate needed data
    
    % infer number of voids in this data table
    iVv  = strncmpi(rawData(1,:),'Vol Void',7);
    NVs = [];
    for k=2:size(rawData,1)
        % Scan down all rows looking for voiding data. When we find a blank
        % or reach the end of the column, assume that is the last void.
        if numel(rawData{k,iVv})==1 && isnan(rawData{k,iVv})
            NVs = k-1;
            break
        elseif k==size(rawData,1)
            NVs = k;
        end
    end
    
    % resume populating voiding data arrays
    Vv  = [rawData{2:NVs,strncmpi(rawData(1,:),'Vol Void',7)}];
    Vi  = [rawData{2:NVs,strncmpi(rawData(1,:),'Vol Inf',7)}];
    dT  = [rawData{2:NVs,strncmpi(rawData(1,:),'Intervoid',7)}];
    VrN = [rawData{2,strncmpi(rawData(1,:),'Last Resid',7)}];
    Vr0 = [rawData{2,strncmpi(rawData(1,:),'Initial Resid',7)}];
end


%% Perform Estimation

% check on data fidelity
if Vr0 + sum(Vi) > sum(Vv) + VrN
    warning(['Data Unreliable: ' ...
             'Infused volume is greater than measured voided and residual volumes; ' ...
             'it will not be possible to satisfy non-negative flow rate constraints.'])
end


switch EstMeth
    case 'ConstantFlow'
        % constant flow rate constraint formulation
        % Danziger & Grill 2016
        
        % compute ureter infusion without constraints based on final PVR (VrN) 
        dVuEst = (sum(Vv)-sum(Vi)-Vr0+VrN)/sum(dT);
        % replicate constant dVuEst to be consistant outputs across estimation methods 
        dVuEst = repmat(dVuEst,[1 length(Vv)]);
        
    case 'MinFlowVariance'
        % non-constant dVu with CV cost function and nonnegativity
        % Angoli, Geramipour, & Danziger 2019
        
        % compute an initial guess:
        % begin with constant dVu estimate
        dVuConst = (sum(Vv)-sum(Vi)-Vr0+VrN)/sum(dT);
        % increment the guess until no Vr estimates are negative
        cc = 0;
        while any( (dVuConst*cumsum(dT) + cumsum(Vi) - cumsum(Vv) + Vr0)<0 ) && cc<1e4
            % add increments of 1% of typical voids (divide by time for rate)
            dVuConst = dVuConst + 0.01*mean(Vv)/mean(dT); 
            cc = cc+1; 
        end
        % create initial guess vector, one element per voiding event
        u0 = zeros(size(dT)) + dVuConst;
        
        % define cost function
        fcost = @(u,dT,Vv,Vr0,VrN) ...              
            (1/mean(u.*dT)^2)*(sum(u.*dT)+sum(Vi)-sum(Vv)+Vr0 - VrN)^2 + std(u)/mean(u);
        
        % define constraints (A*dVu < b) and all dVu > 0 (lbm)
        lbm = zeros(size(u0));                    	
        Am = -tril(repmat(dT,[length(dT) 1]));      
        bm = Vr0+cumsum(Vi)-cumsum(Vv);            	
        
        % silence output
        options = optimoptions('fmincon','Display','none');
        
        % use the same constraints as CV
        exflg = -1; cc = 0;
        while exflg<0 && cc<10
            % flow estiamtes
            [dVuEst,~,exflg] = fmincon(@(u) fcost(u,dT,Vv,Vr0,VrN), u0,Am,bm,[],[],lbm,[],[],options);
            % reseed initial conditions and try again if we failed to find a solution 
            u0 = zeros(size(dT)) + dVuConst + rand(size(dT))*dVuConst;
            
            % break the loop if we are not able to minimize the cost
            % function without violating constraints within 10 attempts
            cc=cc+1;
            if cc==10
                warning('Unable to find ureter flow rate estimates that satisfy constraints!')
            end
        end
        
     
    otherwise
        error('Unecognized Estimation Method')
        
end


% -- Generate Supporting Outputs --
% compute all PVR estimates
VrEst = cumsum(dVuEst.*dT) + cumsum(Vi) - cumsum(Vv) + Vr0;
% compute all Voiding Efficiency estimates
VEEst = Vv./(Vv+VrEst);






