#!/usr/bin/perl
#
use POSIX;
use List::Util;

sub print_help {
    print STDERR <<EOF;
usage:	histogram filename [-i interval] [-w]
	generate the histogram of column 1 of a numeric data file
	default interval based on a number of equi-spaced divisions equaling
	        min[3+n^(1/3),1000], where n=number of available data
	-i interval : set the desiderd interval for the histogram
	-w          : read the next input column as weights
output: 3 columns with x N[x] f[x]
	where x	   = position of the center of the interval
	      N[x] = number of occurences in the interval centered at x
	      f[x] = normalized number of occurrences: f[x] integrates to 1
			      			by Nick Manini, 25-02-2014
EOF
return;
}

$BIG=9.99e99;
$intervalproposed=0;
$weightFlag = 0;

while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){#filename unless starting with -
        if(-e $p){
            push(@file_di_input,$p);
        }else{
            print STDERR "File $p does not exist!!!\n";
        }
        next;
    }
    if($p eq "-h"){# [-h] print help
        print_help;
        exit;
    }
    if($p eq "-w"){# [-w] takes into account a coulumn of weights
        $weightFlag = 1;
        next;
    }
    if($p eq "-i"){# [-i int] forces to build a coulumn at interval int
        $intervalproposed = shift(@ARGV);
        next;
    }
    print_help;
    print STDERR "invalid option $p!\n ";
    exit;
}

if($#file_di_input < 0){
   $filename="/tmp/histogram_tmp".$$;
   `cat > $filename`;
   push(@file_di_input,$filename);
   $rmflag=1;
}

$mainfile="/tmp/histogram".$$;
`cat "@file_di_input" > $mainfile`;# put everything in a single file

# preliminary computation of ndata, min and max:
$min=$BIG;
$max=-$BIG;
$ndata=0;
open(INPUT,$mainfile) || die "cannot open file $mainfile\n";
while(<INPUT>){
    s/^\s*//;	# remove leading spaces
    my @fi =  split(/\s+/);
    if($#fi >= $weightFlag && !/^#/){	# ignore empty and commented lines
	$ndata++;
	if ($fi[0]<$min){$min=$fi[0];}
	if ($fi[0]>$max){$max=$fi[0];}
    }
}
close(INPUT);

if($ndata<=0){die "histogram: no data available to process\n"}

if($intervalproposed==0){
    $ninterval=POSIX::floor(3.+pow($ndata,0.333333333333333));
    if($ninterval>1000){$ninterval=1000;}
    $interval=($max-$min)/$ninterval;
}else{
    $interval=$intervalproposed;
    $ninterval=POSIX::floor(0.+($max-$min)/$interval)
}
$max+=0.75*$interval;
$min-=0.75*$interval;

$delta=$max-$min;
for($i=300;$i>-300;$i--)       # the range of double                   
{
    $delscale=10**$i;
    if($delscale<$interval){last}
}

$count=0;
@facts=(2.,10.,20.,100.);
while($interval/$delscale>1.5){
    $count++;
    $delscale*=$facts[$count]/$facts[$count-1];
}
if($min>0){$scalemin=$delscale*(int($min/$delscale+0.5)-1)}
else      {$scalemin=$delscale*(int($min/$delscale-0.5)-1)}

if($intervalproposed==0){
    $interval=$delscale;
}

$min=$scalemin;
$ninterval=POSIX::floor(($max-$min)/$interval);

for($i=0;$i<=$ninterval;$i++){
    push(@noccurrences,0);
}

open(INPUT,$mainfile) || die "cannot open file $mainfile\n";
if($weightFlag){
    $totweight=0;
    while(<INPUT>){
	s/^\s*//;	# remove leading spaces
	my @fi =  split(/\s+/);
	if($#fi >= 1 && !/^#/){	# ignore too short and commented lines
	    $noccurrences[POSIX::floor(($fi[0]-$min)/$interval)]+=$fi[1];
	    $totweight+=$fi[1];
	}
    }
}else{
    $totweight=$ndata;
    while(<INPUT>){
	s/^\s*//;	# remove leading spaces
	my @fi =  split(/\s+/);
	if($#fi >= 0 && !/^#/){	# ignore empty and commented lines
	    $noccurrences[POSIX::floor(($fi[0]-$min)/$interval)]++;
	}
    }
}
close(INPUT);

if(abs($totweight)>1.e-30){# so that histogram integrates to 1
    $normalization=1.0/($interval*$totweight);
}else{ # in case the weights total 0, one has better not divide!!
    $normalization=1.0/$interval;
}

$ntot=List::Util::sum @noccurrences;
print "# histogram of ",$ntot," data";
if($weightFlag){
    print " using weights, of total value ",$totweight,"\n";
}else{
    print "\n";
}

for($i=0;$i<=$ninterval;$i++){
    print $min+($i+0.5)*$interval,"\t",$noccurrences[$i],"\t",$normalization*$noccurrences[$i],"\n";
}

`rm $mainfile`;
if($rmflag){`rm $filename`;}