Post-process output from sig_RecessionAnalysis

[RY4GIT]

TOSSH/TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m

Lines 162 to 166 in a7012da

	RecessionParameters(i,1) = median((RecessionParametersTemp(:,1)),'omitnan');
	RecessionParameters(i,2) = median(RecessionParametersTemp(:,2),'omitnan');
	RecessionParametersT0Temp = 1./(RecessionParametersTemp(:,1).*median(Q_mat{i}(Q_mat{i}>0),'omitnan').^(RecessionParametersTemp(:,2)-1));
	ReasonableT0 = and(RecessionParametersTemp(:,2)>0.5,RecessionParametersTemp(:,2)<5);
	RecessionParameters(i,3) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan');

Issue description:

The current code ./TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m post-processes recession parameters (a and b) from individual fit from sig_RecessionAnalysis to derive the site median a and b, as well as the timescale of recession decay T0, that can be calculated from a and b. However, this post-processing is not included in calc_ALL.m.

This causes some confusions, as:

• RecessionParameters is a 1-by-2 matrix if you use calc_ALL.m and don't postprocess, and 1-by-3 matrix if you use calc_McMillan_Groundwater.m and post-process the parameters
• RecessionParameters is a 1-by-2 or 1-by-3 matrix, whose element get converted as "RecessionParmaeters_1", "RecessionParmaeters_2", and "RecessionParmaeters_3" in output CSV file by default in Matlab. Therefore, in the output CSV file, it is not labeled which item represent what parameter, either a, b, or T0.
• Additionally, T0 is a scaled timescale of recession by flow magnitude, but the process is not described in comments.

I'm not looking for an immediate fix right now—left the comment for future improvement.

Some improvements to be done:

• In calculation functions (e.g., calc_ALL.m, calc_McMillan_Groundwater.m), instead of storing the output as RecessionParameters (n-by-2 or n-by-3 matrix ), prepare output variables as literal name of each parameters (n-by-1 matrix for each: RecessionParameters_a, RecessionParameters_b, and RecessionParameters_T0) and take the an element of RecessionParameters output like:

    [RecessionParametersTemp,~,~,RecessionParameters_error_str(i)] = ...
        sig_RecessionAnalysis(Q_mat{i},t_mat{i}, ...
        'recession_length', recession_length, ...
        'n_start', n_start, ...
        'eps', eps, ...
        'fit_individual',true);

    % Get the estimated parameters for the power function
    RecessionParameters_a(i) = median((RecessionParametersTemp(:,1)),'omitnan');
    RecessionParameters_b(i) = median(RecessionParametersTemp(:,2),'omitnan'); % This is the b parameter: RecessionParameters_b
    
    % Calculate the timescale of the recession from the parmaeters a and b
    RecessionParametersT0Temp = 1./(RecessionParametersTemp(:,1).*median(Q_mat{i}(Q_mat{i}>0),'omitnan').^(RecessionParametersTemp(:,2)-1)); 
    ReasonableT0 = and(RecessionParametersTemp(:,2)>0.5,RecessionParametersTemp(:,2)<5);
    RecessionParameters_T0(i) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan'); % This is the T0 parameter: RecessionParameters_T0

• The process of scaling by median Q can be either (1) included in ./TOSSH_code/calculation_functions/calc_All.m as well, or (2) included as an option in sig_RecessionAnalysis.
• Also additional commenting and citation to this conversion would be helpful. The post-processing converts the estimated a parameter; it scales a according to flow magnitude (in other words, a gets normalized by median flow) as described in McMillan et al., (2014). b is (supposed to be) insensitive to flow magnitude, so no scaling needed.

[RY4GIT]

Temporary solution has been worked out in this commit on the fork

RY4GIT@5fd7060

It saves the output variables as literal name of each parameters (RecessionParameters_a, RecessionParameters_b, and RecessionParameters_T0).

Also I have sorted out the conversion equation with a help from Hilary in my handwirtten note, but haven't coded them as comment yet

[RY4GIT]

• calc_ALL
• calc_McMillan_Groundwater
• Description in the sig_RecessionAnalysis

[RY4GIT]

Addressed in #14

[SebastianGnann]

This is great, thank you!!

I have two questions. First, would it make sense to write a short signature function that calculates T0? Or provide it as an output with the sig_RecessionAnalysis function? That might perhaps be more consistent than the post-processing.

And second, if we extract the median of the recession parameters a and b separately, this means they come from different recession segments. To get the parameters from the same recession, I used this code snippet:
[~, ind] = min(abs(Recession_Parameters(:,2) - median(Recession_Parameters(:,2),'omitnan')));
Recession_Parameters = [Recession_Parameters(ind,1) Recession_Parameters(ind,2)];
Not sure what the best solution is, but that just came to my mind.

[RY4GIT]

Thank you Sebastian for checking it through!

First, would it make sense to write a short signature function that calculates T0? Or provide it as an output with the sig_RecessionAnalysis function?

I agree either of the options would be more helpful and consistent!

I wasn't entirely sure how the recession signature works with different options so I just went ahead with post processing for temporary solution.

if we extract the median of the recession parameters a and b separately, this means they come from different recession segments. To get the parameters from the same recession, I used this code snippet:
[~, ind] = min(abs(Recession_Parameters(:,2) - median(Recession_Parameters(:,2),'omitnan')));
Recession_Parameters = [Recession_Parameters(ind,1) Recession_Parameters(ind,2)];

Good idea to use it from the same recession, based on median b! It's more like you gonna pick a representative recession out of all recessions, rather than median of the parameters a and b, though. Should we phrase it like a and b from median curve or recession?

I'll be out of office next 2 weeks (I'm actually writing this from an airport), so if you wanna work on the above changes please go ahead! If not I'm gonna pick up later