%{ 
CaliPro implementation for lymph node infection model used 
for the paper: 
    Distinguishing multiple roles for T-cell and Macrophage
    involvement in determining lymph node fates during 
    Mycobacterium tuberculosis infection

    Paper authors: Kathryn C. Krupinsky, Christian T. Michael, 
                   Pariksheet Nanda, Joshua Mattila, and Denise Kirschner

Code written by: Pariksheet Nanda 
Code edited by: Paul Wolberg

The implementation of a CaliPro "class" object described in this file
is an extension of the CaliPro method published by Louis Joslyn et al. as:
    "CaliPro: A Calibration Protocol That Utilizes Parameter Density Estimation 
    to Explore Parameter Space and Calibrate Complex Biological Models"
( https://doi.org/10.1007/s12195-020-00650-z )

Originally implemented in MATLAB for a specific model solving
framework, this code has been edited to highlight where a 
user's calls to an external ODE platform will be made - 
you can search for these in this code by searching for "USER-PLAT"

Please note that this file has been modified and will not work "out of the box". 
Your parameter names, experimental data files, and calls to your solver framework
may differ from those referred to in the examples, 
based on your specific implementation.
%}
classdef CaliPro < handle
    % Calculates new parameters ranges for model constraints.
    %
    % Uses a Lymph node model run result and a set of constraints to adjust
    % parameter ranges using either methods of alternate density
    % subtraction (ADS).  There is only one public method to write the
    % parameter file, but interactively changing the properties provides
    % useful insights into how changing the set of constraints propagates
    % through to the pass percentage and adjusted parameter ranges.
    %
    % Dimension symbols:
    %
    % H = hosts.
    % L = lymph nodes.
    % S = samples = H x L.
    % SP = samples in pass set.
    % SF = samples in fail set.
    % C = constraints.
    % PA = all parameters.
    % PV = varied parameters.
    % PP = varied parameter pool chosen for calibration (e.g. LN.hs1).
    % P = calibration parameters expanded from parameter pool (e.g. 
    %     LN1.hs1, LN2.hs1, .., LN5.hs1).
    %
    % Example:
    %
    %   % Use a reproducible seed to run lymph node model:
    %   rng(182819758, 'twister');
    %
    %   % Run lymph node model:
    %   USER-PLAT_SOLVER_CALL_HERE('lymph_granuloma_mult_ln_settings.m', 'Model_LHS.mat');
    %
    %   % Initialize CaliPro with the lymph node model run:
    %   cal = CaliPro(result);
    %
    %   % Check how well the run meets all constraints:
    %   cal.pass_percent
    %   table(cal.constraints_names', cal.pass_percents, ...
    %         'VariableNames', {'Constraint' 'Pass %'})
    %
    %   % See new parameter ranges:
    %   cal.params_calibrate
    %   cal.ranges
    %   cal.ranges_diff_percent
    %
    %   % Write out the new settings file with the modified parameters:
    %   cal.write_settings_file('lymph_granuloma_mult_ln_settings.m', ...
    %       'calipro_iter_1_settings.m')
    %   % See an output settings file:
    %   system 'cat calipro_iter_1_settings.m';
    %
    %   % Check how well the run meets only the experimental constraints:
    %   cal.constraints(4:end) = 0;
    %   cal.pass_percent
    %
    %   % Cleanup the generated file:
    %   delete calipro_iter_1_settings.m
    %
    % References:
    %
    % [1] Joslyn, Louis R., Denise E. Kirschner, and Jennifer J.
    %     Linderman. 2021. "CaliPro: A Calibration Protocol That Utilizes
    %     Parameter Density Estimation to Explore Parameter Space and
    %     Calibrate Complex Biological Models." Cellular and Molecular
    %     Bioengineering 14 (1): 31–47.
    %     https://doi.org/10.1007/s12195-020-00650-z.
    %
    % [2] Joslyn, Louis R., Jennifer J. Linderman, and Denise E. Kirschner.
    %     2022. "A Virtual Host Model of Mycobacterium Tuberculosis
    %     Infection Identifies Early Immune Events as Predictive of
    %     Infection Outcomes." Journal of Theoretical Biology 539 (April):
    %     111042. https://doi.org/10.1016/j.jtbi.2022.111042.
    %
    % Copyright (c) 2022-2023 Denise Kirschner Laboratory at the University
    % of Michigan Medical School.

    % Public editable properties.
    properties (SetObservable, AbortSet)
        % Logical vector of constraints.  See constraint names for a
        % description of these indices.
        constraints(4, 1) logical = [1; 1; 1; 1]
        % Character vector of P parameters to be calibrated.  Must be a
        % subset of params_varied.
        params_calibrate
    end % properties, public, editable.

    % Public read-only constants.
    properties (Constant)
        % Names of constraints.
        constraints_names = [...
            "CFUs", ...
            "Macs", ...
            "CD4+ T-cells", ...
            "CD8+ T-cells" ...
            ]
    end % properties, public, read-only.

    % Public read-only properties.
    properties (SetAccess = protected)
        % Result structure from a completed run of lymph node model.
        result
        % Product of pass_percents for constraints applied.
        pass_percent double = NaN
        % Pass percent of granulomas for each constraint, including
        % contraints that are not applied.
        pass_percents(4, 1) double = NaN
        % H x L logical matrix of hosts x lymph nodes that meet the set
        % specified by the C x 1 constraints vector.
        pass_granulomas
        % P x 2 adjusted parameter ranges by ADS.
        ranges
        % P x 2 percent difference between ranges before and after
        % calibration.
        ranges_diff_percent
        % PV table of all varied parameters that are potentially 
        % calibratable with their name, group, and pool.
        params_varied
    end % properties, public, read-only.

    % Hidden properties.
    properties (SetAccess = protected, GetAccess = protected)
        % SP x PV matrix of parameters in the pass granulomas set.
        % SP + SF = S (samples) = H x L.
        params_pass_
        % SF x PV matrix of parameters in the fail granulomas set.
        % SP + SF = S (samples) = H x L.
        params_fail_
        % H x PV floating point matrix of all varied parameters.
        params_samples_
        % C x H x L logical matrix of constraint per granuloma.
        pass_constraint_matrices_
        % 1 x P integer vector containing PV indices of params_varied.
        idx_calibrate_
        % 1 x P integer vector containing PA indices of
        % odeSettings.parameters.
        idx_params_varied_
    end % properties, private.

    % Public methods.
    methods
        function obj = CaliPro(result, varargin)
            % Initialize the CaliPro object with the lymph node model result, any
            % additional arguments, and run all calculations to get new
            % parameter ranges.  The event listeners run any affected
            % calculations from updating the public properties.
            %
            % Input:
            %
            % result: Result structure from a completed run of lymph node model.
            %
            % Optional Input:
            %
            % contraints: C x 1 logical vector of constraints.  See
            %             constraint names for a description of these
            %             indices.
            %
            % Optional Input Intended for Internal Use with Unit Tests:
            %
            % params_calibrate: 1 x P character vector of parameters to be 
            %                   calibrated.  Must be a subset of 
            %                   params_varied.
            %
            % force_pass_granulomas: H x L logical matrix of hosts x lymph 
            %                        nodes.
            %
            % Output:
            %
            % obj: CaliPro object that automatically recalculates output
            %      parameters based on changing any of its properties.

            % Assign the lymph node model run.
            obj.result = result;
            % Assign any extra arguments.
            pre_assign_params_calibrate = {};
            defer_assign_pass_granulomas = logical([]);
            if nargin > 1
                while numel(varargin) > 0
                    switch varargin{1}
                        case 'constraints'
                            obj.constraints = varargin{2};
                        case 'params_calibrate'
                            pre_assign_params_calibrate = ...
                                varargin{2};
                        case 'force_pass_granulomas'
                            % Only meant for unit tests!
                            defer_assign_pass_granulomas = ...
                                varargin{2};
                        otherwise
                            msg = ['CaliPro constructor argument %s ' ...
                                'not understood'];
                            err = 'CaliPro:InvalidConstructorArgument';
                            MException(err, msg, varargin{1});
                    end
                    if numel(varargin) > 2
                        varargin = varargin(3:end);
                    else
                        varargin = {};
                    end
                end
            end
            % Initialize the remaining fields.
            obj.initialize_fail_constraint_matrices_();
            obj.calc_pass_granulomas_();
            if numel(defer_assign_pass_granulomas)
                obj.pass_granulomas = ...
                    defer_assign_pass_granulomas;
            end
            obj.initialize_params_varied_();
            % Update parameters if requested.
            if numel(pre_assign_params_calibrate)
                obj.params_calibrate = ...
                    pre_assign_params_calibrate;
            end
            obj.initialize_params_calibrate_();
            obj.initialize_params_samples_();
            obj.initialize_split_samples_();
            obj.calc_ranges_();
            % Register public property setter listeners.
            addlistener(obj, 'constraints', 'PostSet', ...
                @CaliPro.handle_constraints_change_);
            addlistener(obj, 'params_calibrate', 'PostSet', ...
                @CaliPro.handle_params_calibrate_change_);
        end

        function write_settings_file(obj, path_settings_ref, ...
                path_settings_new)
            % Write parameters file with adjusted ranges.
            %
            % Input:
            %
            % file_params: The output file .m to write with updated
            %              parameters, similar in content to
            %              lymph_granuloma_settings.m.

            % Read in the entire reference file.
            str_settings = readlines(path_settings_ref);
            % File location of parameters.
            str_idx = @(str, expr) find(...
                ~ cellfun(@isempty, regexp(str, expr)), 1);
            idx_params_start = str_idx(...
                str_settings, '^[ ]*ode[.]parameters[ ]*=[ ]*');
            idx_params_end = str_idx(...
                str_settings(idx_params_start:end), '^[ ]*}[ ]*;') ...
                + idx_params_start - 1;
            % Collect varied parameter names, values, and file locations.
            tbl_file = table(...
                'VariableNames', ...
                {'param', 'lb', 'ub', 'idx_f'}, ...
                'VariableTypes', ...
                {'string', 'double', 'double', 'int16'}, ...
                'Size', [height(obj.params_varied) 4]);
            % Using two successive single-quotes escapes the single-quote.
            fprintf("USER-PLAT - Your specific framework syntax to find configure parameters goes here.\n");
            %{
            expr = ...
                ['[{]\s*'...
                '''(?<param>\S+)''' ...
                '\s*[,]\s*' ...
                '''(?<distr>\S*)''' ...
                '\s*[,]\s*' ...
                '(?<lb>\S+)' ...
                '\s*[,]\s*' ...
                '(?<ub>\S+)' ...
                '\s*[}]'];
            %}
            idx = idx_params_start;
            tbl = 1;
            while idx < idx_params_end
                % Loop through the file until we encounter a varied 
                % parameter.
                match = [];
                while isempty(match) && idx < idx_params_end
                    match = regexp(str_settings(idx), expr, 'names');
                    idx = idx + 1;
                    % Ignore unvaried parameters
                    if ~ isempty(match) && match.distr == ""
                        match = [];
                    end
                end
                if ~isempty(match) && idx <= idx_params_end
                    tbl_file(tbl, :) = ...
                        {match.param, match.lb, match.ub, idx - 1};
                    tbl = tbl + 1;
                end
            end
            fprintf("USER-PLAT - Your specific framework syntax to write a new parameter file goes here.\n");
            %{
            assert(height(tbl_file) == height(obj.params_varied), ...
                'CaliPro:write_settings_file', ...
                ['Number of varied parameters in template file ' ...
                '%s is %d but the number of varied parameters ' ...
                'in the lymph node model results file is %d; ' ...
                'the numbers should be the same.'], ...
                path_settings_ref, ...
                height(tbl_file), ...
                height(obj.params_varied));
            %}
            %{
            % Overwrite the file lines with the CaliPro adjusted ranges.
            expr = [
                '[,]\s*' ...
                '(?<lb>\S+)' ...
                '\s*[,]\s*' ...
                '(?<ub>\S+)' ...
                '\s*[}]'];
            for idx_p = 1:numel(obj.params_calibrate)
                replacement = [ ...
                    ', ' ...
                    char(vpa(obj.ranges(idx_p, 1))) ...
                    ', ' ...
                    char(vpa(obj.ranges(idx_p, 2))) ...
                    ' }'];
                mask = obj.params_varied.idx_p == idx_p;
                idx_fs = tbl_file.idx_f(mask)';
                for idx_f = idx_fs
                    str_settings(idx_f) = ...
                        regexprep(str_settings(idx_f), expr, replacement);
                end
            end
            %}
            % Write the new file.
            write_lines(str_settings, path_settings_new);
        end
    end % methods, public.

    % Public static method to expose local functions for unit tests as
    % described in:
    % https://www.mathworks.com/matlabcentral/answers/361609-testing-private-functions-in-classes#answer_327241
    methods (Static)
        function function_handles = localfunctions()
            function_handles = localfunctions();
        end
    end % methods, public, static.

    methods (Static, Access = protected)
        % Event handlers.

        function handle_constraints_change_(~, event)
            obj = event.AffectedObject;
            obj.calc_pass_();
            obj.calc_ranges_();
        end

        function handle_params_calibrate_change_(~, event)
            obj = event.AffectedObject;
            % The validation has to be shoehorned in here instead of in the
            % properties block at the top of the class, because validation
            % there requires using constant arguments whereas params needs
            % to be compared against the non-constant argument result.
            obj.validate_params_calibrate_();
            % Now run the usual recalculations affected by changing the
            % params.
            obj.calc_idx_varied_();
            obj.calc_ranges_();
        end

        function handle_method_change_(~, event)
            obj = event.AffectedObject;
            obj.calc_ranges_();
        end

        % Static methods which don't need access to the class object.

        function validate_params_calibrate_(params_calibrate, ...
                params_varied)
            n = numel(params_calibrate);
            found = zeros(1, n);
            for param = 1:n
                found(param) = ismember(...
                    params_calibrate{param}, params_varied);
            end
            if ~all(found)
                err = ...
                    'validate_params_calibrate_:ParamsNotFound';
                msg = ...
                    'Could not find these parameters among params_varied:';
                msg = [msg, repmat(' %s', 1, sum(~found))];
                throw(MException(err, msg, params_calibrate{~found}));
            end
        end
    end % methods, private, static.

    % Hidden methods.
    methods (Access = protected)
        function initialize_fail_constraint_matrices_(obj)
            % Calculate matrices for each constraint.
            %
            % This trades off memory by having a matrix for each constraint
            % to avoid recalculations when playing with different
            % constraints.  The matrices evaluate model result for
            % granulomas meeting a set of experimental and intuitive
            % constraints.
            %
            % Modifies the class object member:
            %
            % pass_constraint_matrices_: C x H x L logical matrix of
            %                            constraint per granuloma.  C is
            %                            the number of constraints, H is
            %                            number of hosts, and L is the
            %                            number of lymph nodes.

            numConstraints = numel(obj.constraints);
            numHosts = obj.result.odeSettings.NR;
            lnPrefix = {'Ln1', 'Ln2', 'Ln3', 'Ln4', 'Ln5'};
            lnOffsets = cellfun(@(s) find(~ cellfun(@isempty, ...
                strfind(obj.result.odeSettings.icNames, s)), 1), ...
                lnPrefix) - 1;
            numLns = numel(lnOffsets);
            obj.pass_constraint_matrices_ = ...
                ones(numConstraints, numHosts, numLns);

            % Load and process experimental data limits.  All experimental
            % data is from the "gauntlet" DropBox data folder.  The .mat
            % files in this directory use the same name as the .csv files
            % from which the data has been memoized.
            fprintf("USER-PLAT - Your specific framework syntax to retrieve experimental data goes here.\n");

            % Granuloma processing loops.
            for host = 1:numHosts
                outputs = ...
                    obj.result.modelOutput(:, host, :);
                % The output index mapping is in the settings file in
                % the ode.initialConditions cell array.
                BlE4C = outputs(2, :, :);
                BlE8C = outputs(6, :, :);

                for ln = 1:numLns
                    % The output index mapping is in the settings file in
                    % the ode.initialConditions cell array.
                    APC = outputs(1 + lnOffsets(ln), :, :);
                    N4C = outputs(2 + lnOffsets(ln), :, :);
                    P4C = outputs(3 + lnOffsets(ln), :, :);
                    E4C = outputs(4 + lnOffsets(ln), :, :);
                    G4C = outputs(5 + lnOffsets(ln), :, :);
                    CM4C = outputs(6 + lnOffsets(ln), :, :);
                    EM4C = outputs(7 + lnOffsets(ln), :, :);
                    N8C = outputs(8 + lnOffsets(ln), :, :);
                    P8C = outputs(9 + lnOffsets(ln), :, :);
                    E8C = outputs(10 + lnOffsets(ln), :, :);
                    G8C = outputs(11 + lnOffsets(ln), :, :);
                    CM8C = outputs(12 + lnOffsets(ln), :, :);
                    EM8C = outputs(13 + lnOffsets(ln), :, :);
                    N4NC = outputs(14 + lnOffsets(ln), :, :);
                    CM4NC = outputs(15 + lnOffsets(ln), :, :);
                    N8NC = outputs(16 + lnOffsets(ln), :, :);
                    CM8NC = outputs(17 + lnOffsets(ln), :, :);
                    MR = outputs(18 + lnOffsets(ln), :, :);
                    MI = outputs(19 + lnOffsets(ln), :, :);
                    MA = outputs(20 + lnOffsets(ln), :, :);
                    BI = outputs(21 + lnOffsets(ln), :, :);
                    BE = outputs(22 + lnOffsets(ln), :, :);

                    % Constraint 1: CFU upper bound after day 90 as
                    % described in the CaliPro [1] and HostSim [2] papers.
                    % CFU = BI + BE.
                    % ----------------------------------UPDATE CONSTRAINT
                    %
                    % Note that the CaliPro paper limits CFU to 36,000 for
                    % the Multigran model, but the HostSim paper methods
                    % *incorrectly* also state a limit of 36,000 which is
                    % clearly more than a magnitude too low than what is
                    % reported in Figure 2B.  So we assume our usual
                    % granuloma limit of 10^6.

                    bugsDaysIndex = intersect(UninfectedBugs.days, ...
                        1:numel(BI(1,:)));
                    if any(BI(:, :, bugsDaysIndex) + ...
                            BE(:, :, bugsDaysIndex)  < ...
                            UninfectedBugs.lowerBound(1:numel(bugsDaysIndex)) | ...
                            BI(:, :, bugsDaysIndex) + ...
                            BE(:, :, bugsDaysIndex) > ...
                            UninfectedBugs.upperBound(1:numel(bugsDaysIndex)))
                        obj.pass_constraint_matrices_(1, host, ln) = 0;
                    end

                    % Constraint 2: Macrophage limits.
        		    macDaysIndex = intersect(UninfectedMacs.days, ...
                        1:numel(MR(1,:)));
                    if any(MR(:, :, macDaysIndex) + ...
                            MI(:, :, macDaysIndex) + ...
                            MA(:, :, macDaysIndex) > abs(1)) && ...
                            any(MR(:, :, macDaysIndex) + ...
                            MI(:, :, macDaysIndex) + ...
                            MA(:, :, macDaysIndex) > ...
                            UninfectedMacs.lowerBound(1:numel(macDaysIndex)) | ...
                            MR(:, :, macDaysIndex) + ...
                            MI(:, :, macDaysIndex) + ...
                            MA(:, :, macDaysIndex) < ...
                            UninfectedMacs.upperBound(1:numel(macDaysIndex)))
                        obj.pass_constraint_matrices_(2, host, ln) = 0;
                    end

                    % Constraint 3: Experimental CD4 T-cell limits.
                    %
                    % The set of T-cells being added may be wrong and Sayed
                    % needs to double check because right now nearly all
                    % runs fail this constraint!
        		    tCellDaysIndex = intersect(...
                        UninfectedCD4Tcells.days, ...
                        1:numel(N4C(1,:)));
                    if any( squeeze(N4C(:, :, tCellDaysIndex) + ...
                            P4C(:, :, tCellDaysIndex) + ...
                            E4C(:, :, tCellDaysIndex) + ...
                            G4C(:, :, tCellDaysIndex) + ...
                            CM4C(:, :, tCellDaysIndex) + ...
                            EM4C(:, :, tCellDaysIndex) + ...
                            N4C(:, :, tCellDaysIndex) + ...
                            CM4C(:, :, tCellDaysIndex) + ...
                            N4NC(:, :, tCellDaysIndex) + ...
                            CM4NC(:, :, tCellDaysIndex))' ...
                            < UninfectedCD4Tcells.lowerBound(...
                            1:numel(tCellDaysIndex)) | ...
                            squeeze(N4C(:, :, tCellDaysIndex) + ...
                            P4C(:, :, tCellDaysIndex) + ...
                            E4C(:, :, tCellDaysIndex) + ...
                            G4C(:, :, tCellDaysIndex) + ...
                            CM4C(:, :, tCellDaysIndex) + ...
                            EM4C(:, :, tCellDaysIndex) + ...
                            N4C(:, :, tCellDaysIndex) + ...
                            CM4C(:, :, tCellDaysIndex) + ...
                            N4NC(:, :, tCellDaysIndex) + ...
                            CM4NC(:, :, tCellDaysIndex))' ...
                            > UninfectedCD4Tcells.upperBound(...
                            1:numel(tCellDaysIndex)))
                        obj.pass_constraint_matrices_(3, host, ln) = 0;
                    end
                    % Constraint 4: Experimental CD8 T-cell limits.
                    %
                    % The set of T-cells being added may be wrong and Sayed
                    % needs to double check because right now nearly all
                    % runs fail this constraint!
        		    tCellDaysIndex = intersect(...
                        UninfectedCD8Tcells.days, ...
                        1:numel(N8C(1,:)));
                    if any( squeeze(N8C(:, :, tCellDaysIndex) + ...
                            P8C(:, :, tCellDaysIndex) + ...
                            E8C(:, :, tCellDaysIndex) + ...
                            G8C(:, :, tCellDaysIndex) + ...
                            CM8C(:, :, tCellDaysIndex) + ...
                            EM8C(:, :, tCellDaysIndex) + ...
                            N8C(:, :, tCellDaysIndex) + ...
                            CM8C(:, :, tCellDaysIndex) + ...
                            N8NC(:, :, tCellDaysIndex) + ...
                            CM8NC(:, :, tCellDaysIndex))' ...
                            < UninfectedCD8Tcells.lowerBound(...
                            1:numel(tCellDaysIndex)) | ...
                            squeeze(N8C(:, :, tCellDaysIndex) + ...
                            P8C(:, :, tCellDaysIndex) + ...
                            E8C(:, :, tCellDaysIndex) + ...
                            G8C(:, :, tCellDaysIndex) + ...
                            CM8C(:, :, tCellDaysIndex) + ...
                            EM8C(:, :, tCellDaysIndex) + ...
                            N8C(:, :, tCellDaysIndex) + ...
                            CM8C(:, :, tCellDaysIndex) + ...
                            N8NC(:, :, tCellDaysIndex) + ...
                            CM8NC(:, :, tCellDaysIndex))' ...
                            > UninfectedCD8Tcells.upperBound(...
                            1:numel(tCellDaysIndex)))
                        obj.pass_constraint_matrices_(4, host, ln) = 0;
                    end
                end
            end
        end

        function calc_pass_granulomas_(obj)
            % Multiply the memoized constraint matrices.
            %
            % Modifies the class object member:
            %
            % pass_granulomas: H x L logical matrix of hosts x initial
            %                  granulomas that meet the set specified by
            %                  the constraints_ array.

            obj.pass_granulomas = logical(squeeze(prod(...
                obj.pass_constraint_matrices_(obj.constraints, :, :))));
            obj.calc_pass_percent_();
        end

        function calc_pass_percent_(obj)
            % Summarize passing granulomas given the constraints.
            %
            % Modifies class object members:
            %
            % pass_percent: Pass percent of granulomas for the selected
            %               constraints.
            %
            % pass_percents: Pass percent of granulomas for each
            %                constraint.

            n = numel(obj.pass_granulomas);
            pass = sum(obj.pass_granulomas, 'all');
            obj.pass_percent = pass / n * 100;
            pass_per_constraint = sum(...
                obj.pass_constraint_matrices_, [2, 3]);
            obj.pass_percents = pass_per_constraint / n * 100;
        end

        function initialize_params_varied_(obj)
            % Extract the names of the varied parameters that can
            % be calibrated from the result object.
            %
            % Modifies the class object member:
            %
            % params_varied: P x 5 table of varied parameter names where P
            %                is the number of parameters whose ranges are
            %                not being calibrated.

            % Read all parameters (e.g. hs1), including those that are not
            % varied, from the ODE file, because the ODE file also assigns
            % the parameter group (e.g. LN1.hs1).  From the group, create a
            % paramter "pool" which strips the group number (.e.g LN.hs1)
            % to compare against the parameters to be calibrated.
            file_ode = ...
                which(func2str(obj.result.odeSettings.odeModelHandle));
            % Read in the entire ODE file.
            tbl = read_ode_grouped_params(file_ode);
            % Filter to only varied parameters.
            params = obj.result.odeSettings.parameters';
            assert(numel(params) == height(tbl), ...
                'CaliPro:NumberOfParametersMismatch', ...
                ['Current ODE function %s lists %d parameters, but ' ...
                'the result contains %d parameters.'], ...
                file_ode, height(tbl), numel(params));
            params = vertcat(params{:});
            tbl.varied = ~ cellfun(@isempty, params(:, 2));
            tbl = removevars(tbl(tbl.varied == 1, :), 'varied');
            % Save the table;
            obj.params_varied = tbl;
        end

        function initialize_params_calibrate_(obj)
            % Read calibration parameters from a lymph node model result.
            %
            % Modifies class object members:
            %
            % params_calibrate: 1 x P character array of parameters chosen
            %                   for calibration during the lymph node model run
            %                   subset from the PV x 5 params_varied
            %                   table.
            %
            % params_varied.idx_p : PV x 1 integer vector containing
            %                       indices of params_calibrate.

            if isempty(obj.params_calibrate) && ...
                ~ismember("parameters_calibrate", ...
                    fieldnames(obj.result.odeSettings))
                err = ...
                    'initialize_params_calibrate_:ParamsCalibrateNotFound';
                msg = ...
                    ['Could not find parameters_calibrate in ' ...
                    'results.  This usually means that the results ' ...
                    'file was generated before the ' ...
                    'parameters_calibrate variable was added, in ' ...
                    'which case you should pass the optional argument ' ...
                    '''params_calibrate'' when initializing the ' ...
                    'CaliPro object.'];
                throw(MException(err, msg));
            end
            if isempty(obj.params_calibrate)
                obj.params_calibrate = ...
                    obj.result.odeSettings.params_calibrate;
            end
            % Check all params in params_calibrate are in params_varied.
            params_varied_pool = ...
                obj.params_varied.pool + "." + obj.params_varied.param;
            [~, idx] = ...
                ismember(params_varied_pool, obj.params_calibrate);
            obj.params_varied.idx_p = idx;
            idx = unique(idx);
            idx = idx(idx ~= 0);
            if numel(idx) ~= numel(obj.params_calibrate)
                err = ...
                    'initialize_params_calibrate_:ParamsDoNotMatch';
                msg = ['All of the requested parameters to be ' ...
                    'calibrated must be present in the list of ' ...
                    'varied parameters!'];
                throw(MException(err, msg)) 
            end
        end

        function initialize_params_samples_(obj)
            % Extract convenience matrix of all parameter LHS samples from
            % the result object.
            %
            % Modifies the class object member:
            %
            % params_samples_: H x PV floating point matrix of all
            %                  varied parameters.

            obj.params_samples_ = obj.result.lhsMatrix;
        end

        function initialize_split_samples_(obj)
            % Split lung granuloma calibration parameter LHS samples into
            % pass and fail set matrices.
            %
            % Modifies class object members:
            %
            % params_pass_: SP x P double matrix of parameter LHS samples
            %               of granulomas pooled across hosts that met the
            %               pass criteria.
            %
            % params_fail_: SF x P double matrix of parameter LHS samples
            %               of granulomas pooled across hosts that failed
            %               the pass criteria.
            
            % Split params_samples array into passed and failed sets.
            num_params = height(obj.params_varied);
            num_passed = sum(obj.pass_granulomas, 'all');
            num_failed = numel(obj.pass_granulomas) - num_passed;
            obj.params_pass_ = zeros(num_passed, num_params);
            obj.params_fail_ = zeros(num_failed, num_params);
            num_hosts = size(obj.params_samples_, 1);
            num_grans = 5;
            i_passed = 1;
            i_failed = 1;
            for host = 1:num_hosts
                for gran = 1:num_grans
                    if obj.pass_granulomas(host, gran)
                        obj.params_pass_(i_passed, :) = ...
                            obj.params_samples_(...
                            host, :);
                        i_passed = i_passed + 1;
                    else
                        obj.params_fail_(i_failed, :) = ...
                            obj.params_samples_(...
                            host, :);
                        i_failed = i_failed + 1;
                    end
                end
            end
        end

        function calc_ranges_(obj)
            % Adjust lung parameter ranges using ADS.
            %
            % Note that no tolerance threshold is applied when updating
            % these ranges because the overall pass percentage is the
            % convergence threshold. Therefore if a range is nearly
            % identical to the original range, it is still tweaked by the
            % ADS algorithm because no assumption is made about the
            % sensitivity of the range.
            %
            % The only time that a range is not tweaked, is if the ADS
            % algorithm fails to calculate a new range due to too few
            % passing parameters.  In this case the original range is used
            % and a warning is thrown.
            %
            % Modifies class object members:
            %
            % ranges: P x 2 double matrix of parameter min and max values.
            %
            % ranges_diff_percent: P x 2 percent difference between ranges
            %                      before and after calibration.

            tbl_params = obj.result.odeSettings.parameters';
            tbl_params = cell2table(vertcat(tbl_params{:}), ...
                'VariableNames', {'param', 'distr', 'lb', 'ub'});
            tbl_params = tbl_params( ...
                ~ cellfun(@isempty, tbl_params.distr), :);
            assert(height(obj.params_varied) == height(tbl_params));
            param_ranges = table2array(tbl_params(:, 3:4));
            num_params = numel(obj.params_calibrate);
            obj.ranges = zeros(num_params, 2);
            ranges_orig = zeros(num_params, 2);
            reverted = ones(num_params, 1) * NaN;
            for param = 1:num_params
                idx = find(obj.params_varied.idx_p == param);
                range = unique(param_ranges(idx, :), 'rows');
                assert(height(range) == 1);
                ranges_orig(param, :) = range;
                pass = obj.params_pass_(:, idx);
                fail = obj.params_fail_(:, idx);
                % Handle the exception of too few or ill-fitting pass sets.
                try
                    obj.ranges(param, :) = find_ads(pass(:), fail(:));
                catch exception
                    err_id = split(exception.identifier, ':');
                    if strcmp(...
                            err_id(end), ...
                            'CannotFitNewRangeFromPassFailSets')
                        obj.ranges(param, :) = ranges_orig(param, :);
                        reverted(param, 1) = param;
                    else
                        rethrow(exception);
                    end
                end
            end
            extents = diff(ranges_orig, 1, 2);
            obj.ranges_diff_percent = ...
                100 * (obj.ranges - ranges_orig) ./ extents;
            reverted = reverted(~isnan(reverted));
            if numel(reverted)
                warning(...
                    'calc_ranges_:RevertedParameterRange', ...
                    ['Reverted %d parameter ranges ', ...
                    'because updating them failed:', ...
                    repmat(' %s', 1, numel(reverted))], ...
                    numel(reverted), ...
                    obj.params_calibrate{reverted});
            end
        end
    end % methods, private.
end % CaliPro class.

function limits = experimentalLimits(fileMat)
% Load and process experimental data limits.
%
% Input:
%
% fileMat: Experimental data is from the "gauntlet" DropBox data folder.
%          The .mat files in this directory use the same name as the .csv
%          files from which the data has been memoized.
%
% Output:
%
% limits: Structure with 3 members of parallel lengths: 'days',
%         'lowerBound' and 'upperBound'.  The 'days' array are the
%         timpoints at which the experimental data was collected,
%         'lowerBound; and 'upperBound' are the relaxed-by-1-magnitude min
%         and max counts at that timepoint.

table = load(fileMat);
name = fieldnames(table);
table = table.(name{1});
repeatedDays = table(:, 1);
repeatedValues = table(:, 2);
assert(issorted(repeatedDays));
limits.days = uint16(unique(repeatedDays));
groups = findgroups(repeatedDays);
% Relax both these limits by 1 order of magnitude per the HostSim paper
% [2].
limits.lowerBound = splitapply(@min, repeatedValues, groups)' / 10;
limits.upperBound = splitapply(@max, repeatedValues, groups)' * 10;
assert(numel(limits.days) == numel(limits.lowerBound));
assert(numel(limits.days) == numel(limits.upperBound));
end

function write_lines(lines, filename)
% The opposite of readlines().
%
% writelines() is only in R2022a but I'm using the older R2021b.  See:
% <a href="matlab:web('https://www.mathworks.com/help/matlab/ref/writelines.html')")writelines</a>.
%
% Input:
%
% lines: String of character vectors, e.g. produced by readlines().
%
% filename: New file to write the lines input to.
%
% Output:
%
% None; writes lines to filename.
    fid = fopen(filename, 'w');
    for i = 1:numel(lines)
        if (i == numel(lines) && isempty(lines{end}))
            % Don't add extra newlines to the end of a file.
            fprintf(fid, '%s', lines{i});
        else
            fprintf(fid, '%s\n', lines{i});
        end
    end
    fclose(fid);
end

function [hdi, xdensityPass, xPass, xdensityFail, xFail, xdensitySub] = ...
    find_ads(passedtotalParamsVec, failedtotalParamsVec)
% For a given set of parameter values, this file will find the density
% distribution of those parameter values, and calculate the ADS according
% to the algorithm set forth by Louis Joslyn in 2020.
%
% Input:
%
% passedtotalParamsVec: Number of passed runs by 1 (the parameter that is
%                       being calculated to find the HDR).
%
% failedtotalParamsVec: Number of failed runs by 1 (the parameter that is
%                       being calculated to find the HDR, the same as the
%                       passedtotalParamsVec parameter).
%
% coverage: Value that determines the region that we are searching for.
%           Value from 0-1, equal to value calibrationInput.portionParams.
%
% Output:
%
% hdi: 1 by 2. The min and max parameter values that corespond to the ADS
%      identified region of parameter space where more pass runs exist than
%      fail runs.  These values will be used as the min and max for
%      sampling the parameter range on the next iteration of CaliPro.
%
% xdensityPass: 1 by 100. 1 by 100 set of density values that describe the
%               underlying pass dataset distribution.
%
% xPass: The xvalues (parameter values) that correspond to the xdensityPass
%        values.  1 by 100.
%
% xdensityFail: 1 by 100. 1 by 100 set of density values that describe the
%               underlying fail dataset distribution.
%
% xFail: The xvalues (parameter values) that correspond to the xdensityFail
%        values.  1 by 100.

if numel(passedtotalParamsVec) == 0
    throw(MException(...
        'find_ads:CannotFitNewRangeFromPassFailSets', ...
        'passedtotalParamsVec is insufficient to propose a new range.'));
end

% Find density estimation of the passed params.
[xdensityPass,xPass] = ks_density_twice(passedtotalParamsVec);

% Find density of the failed params.
%
% If there are no failedtotalParamsVec:
if (isempty(failedtotalParamsVec))
    % Make the density empty.
    xdensityFail = [];
    % Make the xvalues that go with the density empty too.
    xFail = [];
else
    [xdensityFail,xFail] = ks_density_twice(failedtotalParamsVec);
end
% Only keep values from density that are included in the bounds of original
% parameter ranges:
xdensityPass(...
    xPass < min(passedtotalParamsVec) | ...
    xPass > max(passedtotalParamsVec)) = 0;
xdensityFail(...
    xFail < min(failedtotalParamsVec) | ...
    xFail > max(failedtotalParamsVec)) = 0;
% Subtract the failed from the passed:
if isempty(xdensityFail)
    throw(MException(...
        'find_ads:CannotFitNewRangeFromPassFailSets', ...
        'empty xdensityFail is insufficient to propose a new range.'));
else
    xdensitySub = xdensityPass - xdensityFail;
end
% Find all points along range where passed runs are greater than failed
% runs, so that we know the ranges where the runs (on the majority)
% performed well:
hdiRange = xPass(xdensitySub > 0);
% Put the min and max into the highest density interval output.
hdi = [min(hdiRange), max(hdiRange)];
if min(passedtotalParamsVec) == max(passedtotalParamsVec)
      hdi = [max(passedtotalParamsVec), max(passedtotalParamsVec)];
end
% If ADS fails to find a new range, throw an exception instead of returning
% an empty range.
if size(hdi, 2) == 0
    throw(MException(...
        'find_ads:CannotFitNewRangeFromPassFailSets', ...
        'passedtotalParamsVec is insufficient to propose a new range.'));
end

end

function [xDen, xVal] = ks_density_twice(paramVec)
% ks_density_twice.m is a function that calls the ksdensity function twice
% with a slight change in bandwidth the second time.  We found, through
% empirical research of comparing distribution histograms against the
% built-in matlab ks density estimator, that  we needed to run the function
% twice, but with the bandwidth prespecified to be 1/2 of what the embedded
% selection process of the bandwidth is for the unedited version.
%
% Input:
%
% paramVec: A vector of parameter values that correspond to either the pass
%           or fail set of runs. (number of varied parameters by number of
%           runs that satisfied the criteria)
%
% Output:
%
% xDen: 1 by 100 set of density values that describe the underlying data
%       distribution.
%
% xVal: The xvalues (parameter values) that correspond to the xDen values.
%       1 by 100.

% Use the built-in matlab function, ksdensity to create a kernel based
% density estimation of the underlying data distribution.
[~, ~, bandwidth1_P] = ksdensity(paramVec);
% Run again, but with the bandwidth prespecified to be 1/2 of what the
% embedded selection process of the bandwidth is for the unedited version.
% Empirical research done with Dr. Marissa Renardy shows that, for whatever
% reason, the bandwidth selected by this algorithm is too smooth.  Dividing
% that number by 1/2 makes the density approximation a much more realistic
% representation of the underlying data distributon.
[xDen, xVal] = ksdensity(paramVec, 'bandwidth', bandwidth1_P / 2);

end
