Smooth histograms

smooth-histogram.py

#!/usr/bin/env python

# Generates a curve for drawing a smooth histogram for a set of real values
# provided on the standard input.

# Copyright (C) 2016  Morten Hustveit
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

import argparse
import math
import sys

parser = argparse.ArgumentParser(description='Binary name classifier')
parser.add_argument('--min', metavar='min', type=float, nargs='?',
                    help='minimum value', default=0.0)
parser.add_argument('--max', metavar='max', type=float, nargs='?',
                    help='maximum value', default=0.0)
parser.add_argument('--smoothness', metavar='smoothness', type=int, nargs='?',
                    help='smoothness of graph', default=10)
parser.add_argument('--subdivisions', metavar='subdivisions', type=int, nargs='?',
                    help='max number of subdivisions on line', default=200)
args = parser.parse_args()

assert args.smoothness > 1

values = map(lambda x: float(x.strip()), filter(lambda x: len(x) > 0, sys.stdin.read().split('\n')))

min_value = min(args.min, min(values))
max_value = max(args.max, max(values))
step_size = (max_value - min_value) / (args.subdivisions - 1)

shape = []
for i in range(0, args.smoothness):
  f = i / (args.smoothness * 0.5)
  shape.append(math.exp(-(f*f)))

scale = 1.0 / (shape[0] + 2.0 * sum(shape[1:]))
scale /= len(values)
shape = map(lambda x: x * scale, shape)

buckets = [0] * (args.subdivisions + args.smoothness * 2)

area = 0

for value in values:
  bucket_idx = int(round((value - min_value) / step_size) + args.smoothness)
  buckets[bucket_idx] += shape[0]
  area += shape[0]
  for i in range(1, args.smoothness):
    buckets[bucket_idx + i] += shape[i]
    buckets[bucket_idx - i] += shape[i]
    area += 2.0 * shape[i]

for value, idx in zip(buckets, range(0, len(buckets))):
  sys.stdout.write("%.9g\t%.9g\n" % ((idx - args.smoothness) * step_size + min_value, value))