hist

#!/usr/bin/env python
import numpy
import fileinput
import argparse
import math
from sys import exit, stderr

parser = argparse.ArgumentParser()
parser.add_argument('-d', '--density', dest='density', 
                    default=False, action='store_true')
parser.add_argument('--min', dest='min', default=None, type=float)
parser.add_argument('--max', dest='max', default=None, type=float)
parser.add_argument('-n', '--number-of-bins', dest='bins', type=int, 
                    help='number of bins', default=None)
parser.add_argument('-w', '--bin-width', type=float, help='bin width',
                    default=0.0)
parser.add_argument('-g', '--gnuplot', help='gnuplot output', 
                    action='store_true', default=False)
parser.add_argument('-t', '--gnuplot-terminal', help='gnuplot plotting terminal', type=str, 
                    default='')
parser.add_argument('-ly', '--logscale-y', help='set logscale for y', 
                    action='store_true', default=False)
parser.add_argument('-lx', '--logscale-x', help='set logscale for x', 
                    action='store_true', default=False)
parser.add_argument('-l', '--lines', help='use lines', 
                    action='store_true', default=False)
parser.add_argument('-u', '--uncertainty', action='store_true',
                    default=False, help='calculate uncertainty')
parser.add_argument('files', metavar='FILE', nargs='*')
args = parser.parse_args()

values = []
for line in fileinput.input(args.files):
    try:
        x = float(line)
    except ValueError:
        filename = fileinput.filename()
        lineno = fileinput.lineno()
        msg = 'floating point conversion error: %s line %i: %s' % (filename, lineno, line)
        stderr.write(msg)
        exit(1)
    values.append(x)


values = numpy.array(values)

if args.min:
    minval = args.min
else:
    minval = min(values)

if args.max:
    maxval = args.max
else:
    maxval = max(values)

if args.bin_width > 0:
    nbins = (maxval-minval)/args.bin_width
else:
    if args.bins is None:
        # So called "Rice Rule".
        nbins = numpy.ceil(2.0*len(values)**(1.0/3.0))
    else:
        nbins = args.bins
if args.logscale_x:
    bins = numpy.logspace(math.log(minval,10), math.log(maxval,10), nbins+1)
else:
    bins = numpy.linspace(minval, maxval, nbins+1)

hist, bins = numpy.histogram(values, bins=bins, density=args.density)

center = (bins[:-1]+bins[1:])/2.0

def display(args):
    for i in xrange(len(hist)):
        if not args.gnuplot and args.logscale_y:
            if hist[i] == 0:
                continue
            output = "%e %e" % (center[i], numpy.log10(hist[i]))
        else:
            if args.density:
                output = "%e %e" % (center[i], hist[i])
            else:
                output = "%e %i" % (center[i], hist[i])
        if args.uncertainty:
            output += " %.6f" % numpy.sqrt(hist[i])
        print output

if args.gnuplot:
    if len(args.gnuplot_terminal) != 0:
        print 'set term %s' % args.gnuplot_terminal
    if args.logscale_y:
        print 'set logscale y'
    if args.logscale_x:
        print 'set logscale x'
    if args.logscale_x and args.logscale_y:
        print 'set logscale xy'
    if not args.uncertainty:
        if args.lines:
            print "plot '-' u 1:2 w lines notitle"
        else:
            print "plot '-' u 1:2 w boxes notitle"
        display(args)
    else:
        if args.lines:
            print "plot '-' u 1:2 w lines notitle, '-' u 1:2:3 w errorbars ls 1 notitle"
        else:
            print "plot '-' u 1:2 w boxes notitle, '-' u 1:2:3 w errorbars ls 1 notitle"
        display(args)
        print 'e'
        display(args)
else:
    display(args)