#!/usr/bin/perl
#
#  histo3.pl   perl script for histogram binning 
#  (c) Hans-Joachim Drescher    hajo@free.fr  
#  licence: GPL Version 2 
#
#
#
#
$opt_c=1;
$opt_d=0;
$opt_w=0;
use Getopt::Long;
use POSIX qw(ceil floor);
GetOptions(qw(dx:f dy:f f logx logy cond:s c:i d:i w:i h norm:i xmin:f xmax:f ymin:f ymax:f nd:f x0:f y0:f s:f ));
$opt_c--;
$opt_d--;
$opt_w--;
if ($opt_h){
print '

   histo3.pl  -c col_x -d col_y -w col_weight -dy dy -dx dx -logx -logy -nd i -norm inorm -x0 X0 -y0 Y0
              -f 
	      
         -c col_x       column of x values [ 1 ]
         -d col_y       column of y values [ 0 means none ]
         -w col_weight  column of weight [weight is 1 if not specified]
         -dx dx         bin width in x (is factor in log-case)  [ .1 ]
         -dy dy         bin width in y (is factor in log-case)  [ .1 ]

         -nd  i         computes dx as number of bins per decade in log case
         -log           log-binning
         -x0            bin-middle of any bin on x-axis [1]
         -y0            bin-middle of any bin on y-axis [1]
         -f             print empty bins 

         -norm ijk     normalization [0]
            k=0   no normalization
            k=1   normalize by number of entries
            k=4   normalize by number of entries in bin
            j=1   by bin-width
            i     not yet implemented 
         -s f     scale by 1/factor
';
exit;
}


$x0=1;
$x0 = $opt_x0 if defined($opt_x0);
$dx=.1;
$dx=$opt_dx if defined($opt_dx);
$dx=10.0**(1.0/$opt_nd) if defined($opt_nd);
if( $opt_logx) {
    $lx0=log($x0);
    $ldx=log($dx);
}

$y0=1;
$y0 = $opt_y0 if defined($opt_y0);
$dy=.1;
$dy=$opt_dy if defined($opt_dy);
$dy=10.0**(1.0/$opt_nd) if defined($opt_nd) && defined($opt_dy);
if( $opt_logy) {
    $ly0=log($y0);
    $ldy=log($dy);
}

$n=0;
while (<>){
    next if /^\s*\#/;
    next if /^\s*$/;
    $n++;
    if ( $n % 1000 == 0 ) { print STDERR "." ;}
    (@A)=split;

    if ( defined($opt_cond)) {
	($u=$opt_cond) =~ s/\$([0-9]+)/\$A\[-1+\1\]/g;
	# print "here",eval($u),"--->$u<---$opt_cond\n";
	if( ! eval($u) ) {next;}
    }
    
    if( $dx <= 0 ) {
	$a=$A[$opt_c];
    } else {
	$a=xbin(ixbin($A[$opt_c]));
    }

    if( $opt_d > 0 ) {
	if( $dy <= 0 ) {
	    $b=$A[$opt_d];   
	} else {
	    $b=ybin(iybin($A[$opt_d]));
	}
    } else {
	$b=' ';
    }
    
    next if ( defined($opt_xmin) && $a < $opt_xmin ) ;
    next if ( defined($opt_xmax) && $a > $opt_xmax ) ;
    next if ( defined($opt_xmin) && $b < $opt_ymin ) ;
    next if ( defined($opt_ymax) && $b > $opt_ymax ) ;
    
    if ( $opt_w > 0 ) {
	$bin{$a}{$b} += $A[$opt_w];     # with weight 
    } else {
	$bin{$a}{$b} += 1;	        # no weight (one)
    }
    $nbin{$a}{$b} += 1;
}

# bin width for normalization #######################
if( $opt_logy ){
    $deltay=$y*(sqrt($dy)-1.0/sqrt($dy));
} else {
    $deltay=$dy;
}
if( $opt_logx ){
    $deltax=$x*(sqrt($dx)-1.0/sqrt($dx));
} else {
    $deltax=$dx;
}


$s=0.0;$s2=0.0;$w=0.0;$ny=0.0;

# min and max for filling up empty bins 
$ymin= 1e30;
$ymax=-1e30;
$xmin= 1e30;
$xmax=-1e30;
foreach $x ( sort {$a<=>$b} keys(%bin)) {
    $xmin=$x if $x<$xmin;
    $xmax=$x if $x>$xmax;
    if ( $opt_d > 0 ) {
	foreach $y ( sort {$a<=>$b} keys(%{$bin{$x}})) {
	    $ymin=$y if $y<$ymin;
	    $ymax=$y if $y>$ymax;
	}    
    }
}

foreach $x ( sort {$a<=>$b} keys(%bin)) {
    $bin{$x}{$ymin}=0 unless defined($bin{$x}{$ymin}) ;
    $bin{$x}{$ymax}=0 unless defined($bin{$x}{$ymax}) ;
}

foreach $x ( sort {$a<=>$b} keys(%bin)) {
    # print empty bins ################
    if( defined($xlast) && $opt_f ) {
	$i1=ixbin($xlast);
	$i2=ixbin($x);
	for ($ii=$i1+1;$ii<$i2;$ii++){
	    $j1=iybin($ymin);
	    $j2=iybin($ymax);
	    for ($jj=$j1;$jj<=$j2;$jj++){
		print xbin($ii)," ",ybin($jj)," 0 0  <--\n";
	    }
	    print "\n";
	}
    }
    
    foreach $y ( sort {$a<=>$b} keys(%{$bin{$x}})) {
	# normalization ##################################
	$z=$bin{$x}{$y};
	$z /= $deltax*$deltay   if $opt_norm/10 % 10 == 1;
	$z /= $n                if $opt_norm    % 10 == 1;
	
	$z /= $nbin{$x}{$y}     if  $nbin{$x}{$y}  != 0 && $opt_norm    % 10 == 4;
	$z /= $opt_s            if defined($opt_s);

	# again empty bins #########################
	if( defined($ylast) && $opt_f ) {
	    $j1=iybin($ylast);
	    $j2=iybin($y);
	    for ($jj=$j1+1;$jj<$j2;$jj++){
		print $x," ",ybin($jj)," 0 0  <--\n";
	    }
	}
	print "$x $y   ",$z," ",$nbin{$x}{$y},"\n";
	$s+=$x*$z;
	$s2+=$x*$x*$z;
	$w+=$z;
	$ylast=$y;
    }
    $ylast=$undefined;
    $xlast=$x;
    print "\n";
}

if ( $n != 0 ) {
    printf STDERR "%s %e %s %e %s %e\n","total:",$n,"mean:",$s/$w,"  ",sqrt($s2/$w-($s/$w)**2),"  ",sqrt($s2/$w-($s/$w)**2)/sqrt($n),"\n";
}


#
# returns number of bin as a function of x
#
sub ixbin {
    if( $opt_logx ){
	return int(log($_[0]/$xmin)/$ldx+0.5);
    } else {
	return int(($_[0]-$xmin)/$dx+0.5);	
    }
}

#
# returns number of bin as a function of y
#
sub iybin {
    if( $opt_logy ){
	return int(log($_[0]/$ymin)/$ldy+0.5);
    } else {
	return int(($_[0]-$ymin)/$dy+0.5);
	
    }
}

#
# return x as function of ix (number of bin)
#
sub xbin {
    if( $opt_logx ) {
	return $xmin*$dx**$_[0];
    } else {
	return $xmin+$dx*$_[0];
    }
}
# return y 
sub ybin {
    if( $opt_logy ) {
	return $ymin*$dy**$_[0];
    } else {
	return $ymin+$dy*$_[0];
    }
}