Calculates overall statistical uncertainty in a TOPAS distribution using the batch method

topas_batch_uncertainty

#!/usr/bin/env python

###############################################################
##  Author:   David Hall
##  Email:    david(dot)hall(dot)physics(at)gmail(dot)com
##  Date:     March 1st, 2016
##  Program:  topas_batch_uncertainty
##
##  This program reads a set of TOPAS result files and uses 
##  the batch method to compute the relative statistical
##  uncertainty in the scored distribution.
##
##  Generally the high-dose region is used to compute an
##  overall uncertainty, and this can be controlled by setting
##  the threshold from the command-line.
###############################################################

from __future__ import print_function
import sys
import logging
import numpy as np
from argparse import ArgumentParser
from topas2numpy import BinnedResult

def main():
    desc = 'Calculates overall statistical uncertainty in a TOPAS distribution using the batch method'
    parser = ArgumentParser(description=desc)
    parser.add_argument('infile', nargs='+', help='TOPAS result files')
    parser.add_argument('-t', '--threshold', type=float, default=50.0,
                        help='Percentage of maximum value above which voxels contribute [default=50]')
    parser.add_argument('--plot', action='store_true')
    args = parser.parse_args()


    data = [BinnedResult(path) for path in args.infile]

    # input validation
    if len(data) <= 1:
        logging.critical('Cannot employ batch method with only one simulation')
        return 1
    try:
        for d in data[1:]:
            assert d.quantity == data[0].quantity
            assert d.unit == data[0].unit
            assert d.dimensions == data[0].dimensions
    except AssertionError:
        logging.critical('Input files are incomparable (e.g. different dimensions)')
        return 1

    if all('Mean' in d.statistics for d in data):
        stat_name = 'Mean'
    elif all('Sum' in d.statistics for d in data):
        stat_name = 'Sum'
    else:
        logging.critical('Not all files contain the same first moment statistic')
        return 1
      
      
    # compute aggregated mean and standard deviation from sub-samples
    N_samples = len(data)

    mu_tot = np.zeros_like(data[0].data[stat_name])
    for d in data:
        mu_tot = np.add(mu_tot, d.data[stat_name])
    mu_tot /= N_samples

    sigma_tot = np.zeros_like(data[0].data[stat_name])
    for d in data:
	sigma_tot = np.add(sigma_tot, np.square(d.data[stat_name] - mu_tot))
    sigma_tot = np.sqrt(sigma_tot / (N_samples * (N_samples-1)))


    # compute average statistical uncertainty
    # Rogers and Mohan 2000 Questions for comparison of clinical Monte Carlo codes. The use of computers in radiation therapy, 
    # 13th int Conf ICCR. ed T Bortfgld and W Schlegel (Heidelberg Springer) 120-122
    abs_threshold = args.threshold * np.amax(mu_tot) / 100.
    mask = mu_tot >= abs_threshold

    rel_unc = sigma_tot[mask] / mu_tot[mask]
    average_sigma = np.sqrt(np.mean(np.square(rel_unc)))


    print('Threshold of {0}% yields {1}/{2} voxels'.format(args.threshold, mu_tot[mask].size, mu_tot.size))
    print('Average statistical uncertainty = {0:.2%}'.format(average_sigma))
    
    
    if args.plot:
        # find centroid of high-dose region
        sh = mu_tot.shape
        x, y, z = np.meshgrid(np.arange(sh[0]), np.arange(sh[1]), np.arange(sh[2]), indexing='ij')
        x0 = round(np.mean(x[mask]))
        y0 = round(np.mean(y[mask]))
        z0 = round(np.mean(z[mask]))

        # plot 2D slice of mean and standard deviation
        import matplotlib.pyplot as plt
        fig = plt.imshow(mu_tot[:,y0,:])
        fig.axes.get_xaxis().set_ticks([])
        fig.axes.get_yaxis().set_ticks([])
        plt.colorbar()
	plt.savefig('mean.eps', bbox_inches='tight')

        plt.clf()
        fig = plt.imshow(sigma_tot[:,y0,:])
        fig.axes.get_xaxis().set_ticks([])
        fig.axes.get_yaxis().set_ticks([])
        plt.colorbar()
        plt.savefig('std.eps', bbox_inches='tight')


if __name__ == '__main__':
    logging.basicConfig(format='%(levelname)s: %(message)s')
    sys.exit(main())