<filename>src/summariseParticleOutput.py<gh_stars>0
from __future__ import division
from __future__ import print_function
# Postprocess the .h5 file output of crimsonParticleTracking.py to summarise the data recorded on particles.
# Will read the particle_config.json file from the simulation folder (which should be the working dir for this
# script) in order to provide per-bin stats. Expects the relative paths from the working dir to the bins to match those
# given in the particle_config.json

from past.utils import old_div
import h5py
import particles
import numpy
import pandas
import matplotlib.pyplot
from matplotlib.backends.backend_pdf import PdfPages
import re
from multi_key_dict import multi_key_dict
import copy


def pairUpSystoleAndDiastoleBins(config_manager):
    systole_regex = re.compile('systole', re.IGNORECASE)
    diastole_regex = re.compile('diastole', re.IGNORECASE)

    systole_diastole_bin_pairs = {}
    paired_bin_joint_names = multi_key_dict()

    for bin_index, bin_name in enumerate(config_manager.particleBinNames()):
        if systole_regex.match(bin_name):
            for bin_index2, bin_name2 in enumerate(config_manager.particleBinNames()):
                bin_name_without_systole_substring = systole_regex.sub('', bin_name)  # removes all occurences of "systole" from the string
                bin_name_without_diastole_substring = diastole_regex.sub('', bin_name2)

                if bin_name_without_systole_substring == bin_name_without_diastole_substring:
                    print("Pairing bins {} and {}, because I think they're systole and diastole of the same spatial location. A combined analysis on these will additionally be performed.".format(bin_name, bin_name2))
                    systole_diastole_bin_pairs[bin_index] = bin_index2
                    systole_diastole_bin_pairs[bin_index2] = bin_index

                    paired_bin_joint_names[bin_index, bin_index2] = systole_regex.sub('systole+diastole', bin_name)  # same as bin_name_without_diastole_substring

    return systole_diastole_bin_pairs, paired_bin_joint_names


def generate_summaries(config_manager):

    data_file_base_name = "{}-particles-{}".format(config_manager.fname(), config_manager.cycle_start())
    data_file_name_h5 = "{}.h5".format(data_file_base_name)

    systole_diastole_bin_pair_indices, paired_bin_names = pairUpSystoleAndDiastoleBins(config_manager)
    number_of_systole_diastole_bin_pairs = old_div(len(systole_diastole_bin_pair_indices),2)

    if 2*number_of_systole_diastole_bin_pairs != len(systole_diastole_bin_pair_indices):
        raise RuntimeError("Odd number of bin index pairs computed. Code must be buggy...")

    figure_output_pdf = PdfPages(data_file_base_name + ".pdf")
    with h5py.File(data_file_name_h5, 'r') as simulationOutputFile:

        number_of_output_dataframe_rows = 3*config_manager.numberOfParticleBins() + 2*number_of_systole_diastole_bin_pairs  # space for combining bin pair data too

        output_data_column_names = ['Mean over Particles', 'Minimum over Particles', 'Maximum over Particles', 'Number of Particles Involved']
        output_data_row_names = ['none'] * number_of_output_dataframe_rows

        particle_summary_data = numpy.zeros((number_of_output_dataframe_rows, len(output_data_column_names)))

        simulation_end_data_field_name = None
        for bin_index, bin_name in enumerate(config_manager.particleBinNames()):

            # on the first pass through the enclosing loop, get the index of the final data field
            # that we need to access on all subsequent calls:
            if simulation_end_data_field_name is None:
                for field_name in list(simulationOutputFile[bin_name]['Final'].keys()):
                    pass
                simulation_end_data_field_name = field_name  # i.e. set this to the last field name from the previous for loop

            current_bin_data = numpy.array(simulationOutputFile[bin_name]['Final'][simulation_end_data_field_name])[:, 1]
            current_bin_residence_times = numpy.array(simulationOutputFile[bin_name + '_particle_residence_time']['Final'][simulation_end_data_field_name])[:, 1]

            nonzero_residency_time_mask = numpy.where(current_bin_residence_times > 0.001)
            current_bin_nonzero_residence_times = current_bin_residence_times[nonzero_residency_time_mask]
            current_bin_data_for_particles_with_nonzero_residence_time = current_bin_data[nonzero_residency_time_mask]

            data_per_unit_time_for_current_bin = old_div(current_bin_data_for_particles_with_nonzero_residence_time, current_bin_nonzero_residence_times)

            # append histogram for current dataset to the pdf:
            number_of_bins = 256
            x_axis_length = 500
            matplotlib.pyplot.hist(data_per_unit_time_for_current_bin, bins=numpy.linspace(0, x_axis_length, num=number_of_bins), density=True)
            matplotlib.pyplot.axis([0, x_axis_length, 0, 0.05])
            matplotlib.pyplot.xlabel("Time-Averaged {} Acquisition Rate".format(config_manager.getDataNameForBin(bin_name)))
            matplotlib.pyplot.ylabel("Probability Density over Particles")
            matplotlib.pyplot.title("Time-Averaged {} Acquisition Rate for bin {}".format(config_manager.getDataNameForBin(bin_name), bin_name))
            matplotlib.pyplot.plot()
            figure_output_pdf.savefig()
            matplotlib.pyplot.cla()

            x_axis_length = 1.5
            matplotlib.pyplot.hist(current_bin_nonzero_residence_times, bins=numpy.linspace(0, x_axis_length, num=number_of_bins),
                                   density=True, facecolor='g')
            matplotlib.pyplot.axis([0, x_axis_length, 0, 10])
            matplotlib.pyplot.xlabel(
                "Residence Time in {}".format(bin_name))
            matplotlib.pyplot.ylabel("Probability Density over Particles")
            matplotlib.pyplot.title(
                "Particle Residence Times in bin {}".format(bin_name))
            matplotlib.pyplot.plot()
            figure_output_pdf.savefig()
            matplotlib.pyplot.cla()


            print(bin_name, "mean per unit time:", numpy.mean(data_per_unit_time_for_current_bin), "min:", numpy.min(data_per_unit_time_for_current_bin), "max:", numpy.max(data_per_unit_time_for_current_bin))

            output_data_row_names[bin_index] = "Time-Averaged {} Acquisition Rate for bin {} (excluding zeros)".format(config_manager.getDataNameForBin(bin_name), bin_name)
            particle_summary_data[bin_index, 0] = numpy.mean(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 1] = numpy.min(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 2] = numpy.max(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 3] = len(data_per_unit_time_for_current_bin)

            shifted_index = bin_index + config_manager.numberOfParticleBins()
            output_data_row_names[shifted_index] = "residence times " + bin_name + " (excluding zeros)"
            particle_summary_data[shifted_index, 0] = numpy.mean(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 1] = numpy.min(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 2] = numpy.max(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 3] = len(current_bin_nonzero_residence_times)

            twice_shifted_index = bin_index + 2*config_manager.numberOfParticleBins()
            output_data_row_names[twice_shifted_index] = "total " + bin_name + " (excluding zeros)"
            particle_summary_data[twice_shifted_index, 0] = numpy.mean(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 1] = numpy.min(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 2] = numpy.max(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 3] = len(current_bin_data_for_particles_with_nonzero_residence_time)


        # Do the combined systole/diastole bin info:
        systole_diastole_bin_pair_indices_mutating = copy.deepcopy(systole_diastole_bin_pair_indices)  # we'll edit this data structure as we go - defensively, let's not break the original...
        output_row_placement_index_shift = 3*config_manager.numberOfParticleBins()
        for bin_index, bin_name in enumerate(config_manager.particleBinNames()):
            once_shifted_index = bin_index + 1*config_manager.numberOfParticleBins()
            twice_shifted_index = bin_index + 2*config_manager.numberOfParticleBins()

            if bin_index in systole_diastole_bin_pair_indices_mutating:  # checks amongst key entries in the dictionary
                # Recover the original total from the previously computed mean by multiplying by the number of particles:
                number_of_particles_involved_in_bin = particle_summary_data[twice_shifted_index, 3]
                total_bin_data = particle_summary_data[twice_shifted_index, 0] * number_of_particles_involved_in_bin
                
                # do the same for the other bin (the systole bin for diastole bins; the diastole bin for systole bins):
                other_bin_index = systole_diastole_bin_pair_indices_mutating[bin_index]
                other_bin_twice_shifted_index = other_bin_index + 2*config_manager.numberOfParticleBins()
                number_of_particles_involved_in_other_bin = particle_summary_data[other_bin_twice_shifted_index, 3]
                total_other_bin_data = particle_summary_data[other_bin_twice_shifted_index, 0] * number_of_particles_involved_in_other_bin

                total_particles_in_bin = number_of_particles_involved_in_bin + number_of_particles_involved_in_other_bin

                # Total plap output for this bin pair:
                combined_bin_mean_data = old_div((total_bin_data + total_other_bin_data), total_particles_in_bin)

                particle_summary_data[output_row_placement_index_shift, 0] = combined_bin_mean_data
                particle_summary_data[output_row_placement_index_shift, 1] = min(particle_summary_data[twice_shifted_index, 1], particle_summary_data[other_bin_twice_shifted_index, 1])
                particle_summary_data[output_row_placement_index_shift, 2] = max(particle_summary_data[twice_shifted_index, 2], particle_summary_data[other_bin_twice_shifted_index, 2])
                particle_summary_data[output_row_placement_index_shift, 3] = total_particles_in_bin

                output_data_row_names[output_row_placement_index_shift] = "total {} (excluding zeros)".format(paired_bin_names[bin_index])


                output_row_for_combined_bin_residence_time = output_row_placement_index_shift + number_of_systole_diastole_bin_pairs

                total_bin_residence_time = particle_summary_data[once_shifted_index, 0] * number_of_particles_involved_in_bin
                other_bin_once_shifted_index = other_bin_index + 1*config_manager.numberOfParticleBins()
                total_other_bin_residence_time = particle_summary_data[other_bin_once_shifted_index, 0] * number_of_particles_involved_in_other_bin

                combined_bin_mean_residence_time = old_div((total_bin_residence_time + total_other_bin_residence_time), total_particles_in_bin)
                particle_summary_data[output_row_for_combined_bin_residence_time, 0] = combined_bin_mean_residence_time
                particle_summary_data[output_row_for_combined_bin_residence_time, 1] = min(particle_summary_data[once_shifted_index, 1], particle_summary_data[other_bin_once_shifted_index, 1])
                particle_summary_data[output_row_for_combined_bin_residence_time, 2] = max(particle_summary_data[once_shifted_index, 2], particle_summary_data[other_bin_once_shifted_index, 2])
                particle_summary_data[output_row_for_combined_bin_residence_time, 3] = total_particles_in_bin

                output_data_row_names[output_row_for_combined_bin_residence_time] = "residence times {} (excluding zeros)".format(paired_bin_names[bin_index])

                # Remove this bin pair so we don't do it twice:
                systole_diastole_bin_pair_indices_mutating.pop(bin_index)
                systole_diastole_bin_pair_indices_mutating.pop(other_bin_index)

                output_row_placement_index_shift += 1


        figure_output_pdf.close()

        output_data_frame = pandas.DataFrame(particle_summary_data, columns=output_data_column_names, index=output_data_row_names)

        print(output_data_frame)

        data_file_name_xlsx = "{}.xlsx".format(data_file_base_name)
        with pandas.ExcelWriter(data_file_name_xlsx) as writer:
            output_data_frame.to_excel(writer)


if __name__ == "__main__":
    config_manager_in = particles.Configuration('particle_config.json')
    generate_summaries(config_manager_in)
<filename>src/summariseParticleOutput.py<gh_stars>0
from __future__ import division
from __future__ import print_function
# Postprocess the .h5 file output of crimsonParticleTracking.py to summarise the data recorded on particles.
# Will read the particle_config.json file from the simulation folder (which should be the working dir for this
# script) in order to provide per-bin stats. Expects the relative paths from the working dir to the bins to match those
# given in the particle_config.json

from past.utils import old_div
import h5py
import particles
import numpy
import pandas
import matplotlib.pyplot
from matplotlib.backends.backend_pdf import PdfPages
import re
from multi_key_dict import multi_key_dict
import copy


def pairUpSystoleAndDiastoleBins(config_manager):
    systole_regex = re.compile('systole', re.IGNORECASE)
    diastole_regex = re.compile('diastole', re.IGNORECASE)

    systole_diastole_bin_pairs = {}
    paired_bin_joint_names = multi_key_dict()

    for bin_index, bin_name in enumerate(config_manager.particleBinNames()):
        if systole_regex.match(bin_name):
            for bin_index2, bin_name2 in enumerate(config_manager.particleBinNames()):
                bin_name_without_systole_substring = systole_regex.sub('', bin_name)  # removes all occurences of "systole" from the string
                bin_name_without_diastole_substring = diastole_regex.sub('', bin_name2)

                if bin_name_without_systole_substring == bin_name_without_diastole_substring:
                    print("Pairing bins {} and {}, because I think they're systole and diastole of the same spatial location. A combined analysis on these will additionally be performed.".format(bin_name, bin_name2))
                    systole_diastole_bin_pairs[bin_index] = bin_index2
                    systole_diastole_bin_pairs[bin_index2] = bin_index

                    paired_bin_joint_names[bin_index, bin_index2] = systole_regex.sub('systole+diastole', bin_name)  # same as bin_name_without_diastole_substring

    return systole_diastole_bin_pairs, paired_bin_joint_names


def generate_summaries(config_manager):

    data_file_base_name = "{}-particles-{}".format(config_manager.fname(), config_manager.cycle_start())
    data_file_name_h5 = "{}.h5".format(data_file_base_name)

    systole_diastole_bin_pair_indices, paired_bin_names = pairUpSystoleAndDiastoleBins(config_manager)
    number_of_systole_diastole_bin_pairs = old_div(len(systole_diastole_bin_pair_indices),2)

    if 2*number_of_systole_diastole_bin_pairs != len(systole_diastole_bin_pair_indices):
        raise RuntimeError("Odd number of bin index pairs computed. Code must be buggy...")

    figure_output_pdf = PdfPages(data_file_base_name + ".pdf")
    with h5py.File(data_file_name_h5, 'r') as simulationOutputFile:

        number_of_output_dataframe_rows = 3*config_manager.numberOfParticleBins() + 2*number_of_systole_diastole_bin_pairs  # space for combining bin pair data too

        output_data_column_names = ['Mean over Particles', 'Minimum over Particles', 'Maximum over Particles', 'Number of Particles Involved']
        output_data_row_names = ['none'] * number_of_output_dataframe_rows

        particle_summary_data = numpy.zeros((number_of_output_dataframe_rows, len(output_data_column_names)))

        simulation_end_data_field_name = None
        for bin_index, bin_name in enumerate(config_manager.particleBinNames()):

            # on the first pass through the enclosing loop, get the index of the final data field
            # that we need to access on all subsequent calls:
            if simulation_end_data_field_name is None:
                for field_name in list(simulationOutputFile[bin_name]['Final'].keys()):
                    pass
                simulation_end_data_field_name = field_name  # i.e. set this to the last field name from the previous for loop

            current_bin_data = numpy.array(simulationOutputFile[bin_name]['Final'][simulation_end_data_field_name])[:, 1]
            current_bin_residence_times = numpy.array(simulationOutputFile[bin_name + '_particle_residence_time']['Final'][simulation_end_data_field_name])[:, 1]

            nonzero_residency_time_mask = numpy.where(current_bin_residence_times > 0.001)
            current_bin_nonzero_residence_times = current_bin_residence_times[nonzero_residency_time_mask]
            current_bin_data_for_particles_with_nonzero_residence_time = current_bin_data[nonzero_residency_time_mask]

            data_per_unit_time_for_current_bin = old_div(current_bin_data_for_particles_with_nonzero_residence_time, current_bin_nonzero_residence_times)

            # append histogram for current dataset to the pdf:
            number_of_bins = 256
            x_axis_length = 500
            matplotlib.pyplot.hist(data_per_unit_time_for_current_bin, bins=numpy.linspace(0, x_axis_length, num=number_of_bins), density=True)
            matplotlib.pyplot.axis([0, x_axis_length, 0, 0.05])
            matplotlib.pyplot.xlabel("Time-Averaged {} Acquisition Rate".format(config_manager.getDataNameForBin(bin_name)))
            matplotlib.pyplot.ylabel("Probability Density over Particles")
            matplotlib.pyplot.title("Time-Averaged {} Acquisition Rate for bin {}".format(config_manager.getDataNameForBin(bin_name), bin_name))
            matplotlib.pyplot.plot()
            figure_output_pdf.savefig()
            matplotlib.pyplot.cla()

            x_axis_length = 1.5
            matplotlib.pyplot.hist(current_bin_nonzero_residence_times, bins=numpy.linspace(0, x_axis_length, num=number_of_bins),
                                   density=True, facecolor='g')
            matplotlib.pyplot.axis([0, x_axis_length, 0, 10])
            matplotlib.pyplot.xlabel(
                "Residence Time in {}".format(bin_name))
            matplotlib.pyplot.ylabel("Probability Density over Particles")
            matplotlib.pyplot.title(
                "Particle Residence Times in bin {}".format(bin_name))
            matplotlib.pyplot.plot()
            figure_output_pdf.savefig()
            matplotlib.pyplot.cla()


            print(bin_name, "mean per unit time:", numpy.mean(data_per_unit_time_for_current_bin), "min:", numpy.min(data_per_unit_time_for_current_bin), "max:", numpy.max(data_per_unit_time_for_current_bin))

            output_data_row_names[bin_index] = "Time-Averaged {} Acquisition Rate for bin {} (excluding zeros)".format(config_manager.getDataNameForBin(bin_name), bin_name)
            particle_summary_data[bin_index, 0] = numpy.mean(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 1] = numpy.min(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 2] = numpy.max(data_per_unit_time_for_current_bin)
            particle_summary_data[bin_index, 3] = len(data_per_unit_time_for_current_bin)

            shifted_index = bin_index + config_manager.numberOfParticleBins()
            output_data_row_names[shifted_index] = "residence times " + bin_name + " (excluding zeros)"
            particle_summary_data[shifted_index, 0] = numpy.mean(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 1] = numpy.min(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 2] = numpy.max(current_bin_nonzero_residence_times)
            particle_summary_data[shifted_index, 3] = len(current_bin_nonzero_residence_times)

            twice_shifted_index = bin_index + 2*config_manager.numberOfParticleBins()
            output_data_row_names[twice_shifted_index] = "total " + bin_name + " (excluding zeros)"
            particle_summary_data[twice_shifted_index, 0] = numpy.mean(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 1] = numpy.min(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 2] = numpy.max(current_bin_data_for_particles_with_nonzero_residence_time)
            particle_summary_data[twice_shifted_index, 3] = len(current_bin_data_for_particles_with_nonzero_residence_time)


        # Do the combined systole/diastole bin info:
        systole_diastole_bin_pair_indices_mutating = copy.deepcopy(systole_diastole_bin_pair_indices)  # we'll edit this data structure as we go - defensively, let's not break the original...
        output_row_placement_index_shift = 3*config_manager.numberOfParticleBins()
        for bin_index, bin_name in enumerate(config_manager.particleBinNames()):
            once_shifted_index = bin_index + 1*config_manager.numberOfParticleBins()
            twice_shifted_index = bin_index + 2*config_manager.numberOfParticleBins()

            if bin_index in systole_diastole_bin_pair_indices_mutating:  # checks amongst key entries in the dictionary
                # Recover the original total from the previously computed mean by multiplying by the number of particles:
                number_of_particles_involved_in_bin = particle_summary_data[twice_shifted_index, 3]
                total_bin_data = particle_summary_data[twice_shifted_index, 0] * number_of_particles_involved_in_bin
                
                # do the same for the other bin (the systole bin for diastole bins; the diastole bin for systole bins):
                other_bin_index = systole_diastole_bin_pair_indices_mutating[bin_index]
                other_bin_twice_shifted_index = other_bin_index + 2*config_manager.numberOfParticleBins()
                number_of_particles_involved_in_other_bin = particle_summary_data[other_bin_twice_shifted_index, 3]
                total_other_bin_data = particle_summary_data[other_bin_twice_shifted_index, 0] * number_of_particles_involved_in_other_bin

                total_particles_in_bin = number_of_particles_involved_in_bin + number_of_particles_involved_in_other_bin

                # Total plap output for this bin pair:
                combined_bin_mean_data = old_div((total_bin_data + total_other_bin_data), total_particles_in_bin)

                particle_summary_data[output_row_placement_index_shift, 0] = combined_bin_mean_data
                particle_summary_data[output_row_placement_index_shift, 1] = min(particle_summary_data[twice_shifted_index, 1], particle_summary_data[other_bin_twice_shifted_index, 1])
                particle_summary_data[output_row_placement_index_shift, 2] = max(particle_summary_data[twice_shifted_index, 2], particle_summary_data[other_bin_twice_shifted_index, 2])
                particle_summary_data[output_row_placement_index_shift, 3] = total_particles_in_bin

                output_data_row_names[output_row_placement_index_shift] = "total {} (excluding zeros)".format(paired_bin_names[bin_index])


                output_row_for_combined_bin_residence_time = output_row_placement_index_shift + number_of_systole_diastole_bin_pairs

                total_bin_residence_time = particle_summary_data[once_shifted_index, 0] * number_of_particles_involved_in_bin
                other_bin_once_shifted_index = other_bin_index + 1*config_manager.numberOfParticleBins()
                total_other_bin_residence_time = particle_summary_data[other_bin_once_shifted_index, 0] * number_of_particles_involved_in_other_bin

                combined_bin_mean_residence_time = old_div((total_bin_residence_time + total_other_bin_residence_time), total_particles_in_bin)
                particle_summary_data[output_row_for_combined_bin_residence_time, 0] = combined_bin_mean_residence_time
                particle_summary_data[output_row_for_combined_bin_residence_time, 1] = min(particle_summary_data[once_shifted_index, 1], particle_summary_data[other_bin_once_shifted_index, 1])
                particle_summary_data[output_row_for_combined_bin_residence_time, 2] = max(particle_summary_data[once_shifted_index, 2], particle_summary_data[other_bin_once_shifted_index, 2])
                particle_summary_data[output_row_for_combined_bin_residence_time, 3] = total_particles_in_bin

                output_data_row_names[output_row_for_combined_bin_residence_time] = "residence times {} (excluding zeros)".format(paired_bin_names[bin_index])

                # Remove this bin pair so we don't do it twice:
                systole_diastole_bin_pair_indices_mutating.pop(bin_index)
                systole_diastole_bin_pair_indices_mutating.pop(other_bin_index)

                output_row_placement_index_shift += 1


        figure_output_pdf.close()

        output_data_frame = pandas.DataFrame(particle_summary_data, columns=output_data_column_names, index=output_data_row_names)

        print(output_data_frame)

        data_file_name_xlsx = "{}.xlsx".format(data_file_base_name)
        with pandas.ExcelWriter(data_file_name_xlsx) as writer:
            output_data_frame.to_excel(writer)


if __name__ == "__main__":
    config_manager_in = particles.Configuration('particle_config.json')
    generate_summaries(config_manager_in)
# Postprocess the .h5 file output of crimsonParticleTracking.py to summarise the data recorded on particles. # Will read the particle_config.json file from the simulation folder (which should be the working dir for this # script) in order to provide per-bin stats. Expects the relative paths from the working dir to the bins to match those # given in the particle_config.json # removes all occurences of "systole" from the string # same as bin_name_without_diastole_substring # space for combining bin pair data too # on the first pass through the enclosing loop, get the index of the final data field # that we need to access on all subsequent calls: # i.e. set this to the last field name from the previous for loop # append histogram for current dataset to the pdf: # Do the combined systole/diastole bin info: # we'll edit this data structure as we go - defensively, let's not break the original... # checks amongst key entries in the dictionary # Recover the original total from the previously computed mean by multiplying by the number of particles: # do the same for the other bin (the systole bin for diastole bins; the diastole bin for systole bins): # Total plap output for this bin pair: # Remove this bin pair so we don't do it twice: