%% Overview
% This script walks through some examples using our published data to
% illustrate how to use the AMEASURE algorithm to estimate post-void
% residuals (PVR) from cystometry data.
%
% Walkthroughs are included for estimating PVR from data stored in Excel
% files or MATLAB variables. Instructions to properly format your Excel
% file for the EstResidVols_v1p2 function, see help EstResidVols_v1p2.
%
% This file contains 3 examples:
% 1) Using AMEASURE algorithms to estimate PVR during an acute anesthetized
% preparation where infusion of the bladder is used to increase voiding
% frequency. In this example, there is only a single batch of voids on
% which to perform the estimation.
% 2) Estimating and plotting PVR in an animal left in a metabolic cage for
% three weeks. Bladder volume was measured approximately every 12 hours,
% creating many batches, each one requring its own PVR estimation.
% 3) A simple simulation of a batch of voids to futher demonstrate the
% estimation method.
%
%
% 
% Note on PVR measurement noise: Sucking bladder volume out through a
% catheter to measure PVR is probably less accurate than you think. To
% get a sense of the noise associated with this type of measurement, you
% can look back at your own data in your cystometry studies. For each void,
% subtract the PVR and voided volume from the volume you infused during the
% bladder fill period, that is Vi-(Vv+Vr). You will find that often this
% quantity is negative, ... but that is impossible!
% The Vv+Vr term must always be at least as large as Vi (and exactly equal
% to it if you have diverted the ureters). This is an example where you can
% see clearly that physical PVR measurement is subject to substantial
% variation. It is discussed in the 2019 paper and incorporated into the
% simulations.
%
%
% ZC Danziger Oct 2019
%



%% Estimating Residual Volumes in Anesthetized Serial Cystometrograms
%
% Data example 1 (Data1-AcuteWithReplacement.xlsx) is an Excel file
% containing acute cystometry on an anesthetized rat. After each void, the
% PVR was removed for measurement, then placed back into the bladder and
% the bladder infusion continued.

% == Run the estimation algorithm ==
Data1 = 'Data1-AcuteWithReplacement.xlsx';
sheetNum = 1;
% Directly compute PVR estimates from the properly formatted Excel data
% file
[EstVr, EstdVu, VEEst] = EstResidVols_v1p2(Data1,sheetNum,'MinFlowVariance');
% display the estimated residuals
disp('Estimated PVRs'); disp(EstVr); disp(' -------- ')

% == Compare the PVR estimates and the measured values: ==

% load data from Excel and place into variables
measuredValues = xlsread(Data1,sheetNum);
measuredPVR = measuredValues(:,7);
measuredVoid = measuredValues(:,3);
measuredVE = measuredVoid ./ (measuredVoid + measuredPVR);

% Compute differences between measurements and estimates of voiding
% efficiency. Typically, measured PVR will be unavailable since the goal of
% the algorithm is to estimate PVR; however, these PVRs were measured as
% part of the algorithm validation protocol.
disp('   Msrd VE   Estmd VE   |error|')
disp([measuredVE, VEEst', abs(measuredVE-VEEst')])
disp(' -------- ')

% Reproduce figure 1A from Danziger, Z. C. and W. M. Grill (2016).
% "Estimating postvoid residual volume without measuring residual bladder
% volume during serial cystometrograms." Am J Physiol Renal Physiol
% 311(2): F459-468.
figure; hold on
plot(VEEst,'-o','linewidth',2,'markersize',10,'MarkerFaceColor','w')
plot(measuredVE,'-o','linewidth',2,'markersize',10,'MarkerFaceColor','w')
ylim([0.5 1])
legend('Estimated','Measured','location','best')
xlabel('Void Number'); ylabel('Voiding Efficiency')
grid on



%% Estimate PVR from a Longitudinal Data Set

% Data2-Longitudinal contains a sereis of void batches over time. Each
% batch is placed in its own sheet in the Excel workbook. At the start of
% each batch, the bladder was (usually) emptied, and the batch is
% considered to be over when we removed the bladder volume for measurement.
% The AMEASURE algorithm estimates the PVR of all voids in between those
% time points. This is reapeated for each batch.

% Data storage parameters:
Data2 = 'Data2-Longitudinal.xlsx';
numSheets = 43;

% Extract voiding data and place into array:
VoidData = [];
for k=1:numSheets
    % estimate PVR
    EstVr = EstResidVols_v1p2(Data2,k,'MinFlowVariance');
    % load in measured data from sheet k
    excelData = xlsread(Data2,k);
    % concatenate measured voiding data with the estimated PVR
    VoidData = [VoidData; [excelData(:,1:6), EstVr']];
    fprintf('Sheet %d of %d loaded.\n',k,numSheets)
end


% Prepare to plot results:
% Separate the data by whether the void occurred in the light cycle (0) or
% the dark cycle (1), stored in column 6.
% Often, when we remove the bladder volume for measurement at the end of a
% batch it is not taken right after the animal actually voids because we do
% not have a scientist continually watching the animal. As a result, this
% final measurement is recorded as a voiding event with a "voided volume" of
% zero with a "residual" of the measured amount (terms in quotes because
% this is not an actual voiding event).
% It is important to make this measurement to account for any extra
% urine accumulation in the bladder since the final actual void in the batch
% to ensure the estimations account for the extra filling, but this extra
% "voiding event" was not produced by the animal and should not be included
% in an analysis of LUT behavior. Therefore, this "false void" event
% is used in our PVR estimation of all other voids, but is not itself
% included in the physiological analysis - all voids with a zero voided
% volume are excluded.
ixLight = VoidData(:,6) == 0 & VoidData(:,2) ~= 0;
ixDark  = VoidData(:,6) == 1 & VoidData(:,2) ~= 0;

% Plot Longitudinal Data:
% figure resize callback
ResClbk = @(u,evt) arrayfun(@(u) datetick(u,'x','mmm-dd'),findobj(get(u,'children'),'type','axes'));       
figure('Name','Longitudinal Voids','SizeChangedFcn',ResClbk,'NumberTitle','off','Position',[270 550 1500 420]); 

subplot(2,1,1); hold on                             % subplot for voided volume
plot(VoidData(ixLight,5),VoidData(ixLight,2),'.')   % plot light cycle voids vs time of day
plot(VoidData(ixDark,5), VoidData(ixDark,2), '.')   % plot dark  cycle voids vs time of day
datetick('x','mmm-dd')                              % change xlabels from datenumbers to interpretable strings
ylabel('Voided Volume (ml)'); title('Longitudinal Voiding Parameters')
ylim([0 1.5])
legend('Light Cycle','Dark Cycle')

subplot(2,1,2); hold on                             % subplot for PVR
plot(VoidData(ixLight,5),VoidData(ixLight,7),'.')   
plot(VoidData(ixDark,5), VoidData(ixDark,7), '.')   
datetick('x','mmm-dd')   
xlabel('Day'); ylabel({'Estimated Post-Void'; 'Residual Volume (ml)'})
ylim([0 1.5])






%% Construct a Simple Simulation of Serial Voiding


% set up simulation voiding parameters
% (notation follows equation 1 in Angoli, Geramipour, & Danziger 2019)
NumVoids = 12;          % number of voids in the batch
dVu = 0.37;             % average volume flow rate through ureters to bladder (ml/h)
dVuStd = 0.075;         % standard deveiation of dVu (ml/h)
t = 45/60;              % average inter-void interval (h)
tStd = 10/60;           % standard deviation of t (h)
Vr = 0.06;              % average PVR in the batch (ml)
VrStd = 0.01;           % standard deviation of Vr (ml)
% set up measurment parameters
VrMsrStd = 0.05;        % noise when actually measuring residual bladder volume (ml)
VvMsrStd = 0.005;       % noise when actually measuring voided volume (ml)

% simulate a batch of voids
n = 1;
% Track the ground-truth residual volumes (VrTrue), the PVR obtained by
% directly removing and measuring bladder volume (VrMsur), and the PVR
% estimated by the algorithm (VrEstm) using only the final measured PVR in
% the batch of voids.
varNames = {'VrTrue','VrMsur','VrEstm','IVI','VvTrue','VvMsur','dVu'};
tmpFill = mat2cell( nan(NumVoids,length(varNames)), NumVoids, ones(1,length(varNames)) );
SimVoids = table(tmpFill{:},'VariableNames',varNames);
while n <= NumVoids
    % if this is the first void, Vr starts at zero, otherwise Vr is
    % whatever was left over from the previous void
    if n==1
        VrPrev = 0;
    else
        VrPrev = SimVoids{n-1,'VrTrue'};
    end
    
    % = generate ground-truth voiding data =
    SimVoids{n,'VrTrue'}    = Vr  + randn*VrStd;	% residual for this void
    SimVoids{n,'IVI'}       = t   + randn*tStd;     % intervoid interval (time bladder has to fill)
    SimVoids{n,'dVu'}       = dVu + randn*dVuStd;	% bladder fill rate from ureters
    % compute what the voided volume must be given all other parameters (see manuscript equation 1):
    SimVoids{n,'VvTrue'}    = SimVoids{n,'dVu'}*SimVoids{n,'IVI'} + VrPrev - SimVoids{n,'VrTrue'};
    
    % = generate simulated voiding data measurements =
    SimVoids{n,'VrMsur'} = SimVoids{n,'VrTrue'} + randn*VrMsrStd;
    SimVoids{n,'VvMsur'} = SimVoids{n,'VvTrue'} + randn*VvMsrStd;
    
    % before moving to the next void, check to make sure no parameters violate conservation laws
    if SimVoids{n,'VrTrue'}>0 && SimVoids{n,'IVI'}>0.1*t && SimVoids{n,'dVu'}>0.1*dVu ...
            && SimVoids{n,'VvMsur'}>0 && SimVoids{n,'VrMsur'}>0
        n = n+1;
    end
end

% = generate estimated PVRs =
% Once the batch is over, algorithmically estimate PVR of voids. Also note,
% we are using the measurement values for estimation (which have noise),
% not the ground-truth values.
% Here we use the array input syntax instead of the Excel input syntax.
[VrEst, dVuEst, VEEst] = EstResidVols_v1p2(SimVoids{:,'VvMsur'}',zeros(1,NumVoids),SimVoids{:,'IVI'}',SimVoids{NumVoids,'VrMsur'}',0,'MinFlowVariance');
% Store data in a simulation table
SimVoids{:,'VrEstm'} = VrEst';

disp(SimVoids)

% Plot the voiding efficiency results:
figure
plot(cumsum(SimVoids{:,'IVI'}),SimVoids{:,'VvTrue'}./(SimVoids{:,'VrTrue'}+SimVoids{:,'VvTrue'}),'k',...
     cumsum(SimVoids{:,'IVI'}),SimVoids{:,'VvMsur'}./(SimVoids{:,'VrMsur'}+SimVoids{:,'VvMsur'}),...
     cumsum(SimVoids{:,'IVI'}),SimVoids{:,'VvMsur'}./(SimVoids{:,'VrEstm'}+SimVoids{:,'VvMsur'}),'linewidth',3)
xlabel('Time (hours)')
ylabel('VE of Each Void (ml)')
legend('Ground-Truth','Measured VE','Estimated VE','location','SE')
set(gca,'box','off','ylim',[0 1]); grid on










