#!/usr/bin/perl

# PURPOSE: Histogram calculation script

# There is NO guarantee that this works in every case

# You can do whatever you want with this script, but
# try to keep these comment lines to give the credit to
# the author

# Version: 3.0
# AUTHOR: P. Nurmi
# DATE: 18.10.2006
# COMMENTS AND FURTHER DEVELOPMENT: contact pasinurmi@netscape.net

# ============================================================================
# START PROGRAM => DEFINITIONS
# ============================================================================

use Getopt::Long;
use POSIX qw(log10);

# Basic Input operations
my $log;
my $sub;
my $min;
my $max;
my $dsample;
my $proces;
my $cumul;
my $times;
my $user;
my $visual;
my $yvals;
GetOptions ( 'logaritmic' => \$log, 'subtract' => \$sub,
	     'dsample' => \$dsample, 'process' => \$proces, "min=f" => \$min,
	     "xmax=f" => \$max, "cumulative" => \$cumul , "times=f" => \$times,
             "user:s" => \$user, 'visualize' => \$visual, 'yvalues' => \$yvals);
my $step = shift;
my $column = shift || 1;my $inputfile= shift;

# Test command-line parameters
unless ( defined( $step )  && 
	 defined( $column )  && defined( $inputfile )) {
  print STDERR "$0: $1: $2: usage: $0 $1 $2

   *******************************************
   *            HISTOGRAM PROGRAM            *
   *******************************************

[ -l (for log10)]
[ -s (to subtract/add/multply/divide 2 columns before calculation)] 
[ -d (to take sample of col 1 based on col 2 and its upper and lower limits)] 
[ -p (NOTHING YET)] 
[ -c (to calculate cumulative disribution)] 
[ -y (take y-values from another column)] 
[ -m= (minimum value included in the calculation)] 
[ -x= (maximum value included in the calculation)] 
[ -t= (multiply y values in the histogram by some arbitrary value)]
[ -u= (user specified additional designation for the files)]  
[ -v (prints the data to the screen and make simple plot of the histogram)]  
  -m=<min> -x=<max> -t=<times> -u=<desig> <step> <column> <inputfile>\n
   At least <step> <column> <inputfile> must be given 
********************************
 OUTPUT:
********************************
 File 1: <inputfile><desig>.hist
       1     , 2,   3   ,    4    ,      5     ,             6
: bin(center), N, N/step, N/all(N), N(no times),[ybin(center)/bin(center)]
       OR  if -l 
       1     ,    2  ,         3        ,    4    ,       5       , 6 ,   7 
: bin(center), N/step, N/(10^upp-10^low), N/all(N), 10^bin(center), N ,[N(y)/N]

********************************
 File 2 (if option -c is given): <inputfile><desig>.cumhist
       1     ,    2      ,            3          ,               4              ,    5  (if -l)
: bin(center), sum(N)<bin, sum(N)*bin(center)<bin, sum(N)*bin(center)/all(N)<bin,10^bin(center) 

********************************
 OTHER REMARKS:
********************************

 N is multiplied by value if -t is give \n";

  exit 1;
}

$column--;
$column2=0;
$column3=0;
$column4=0;
$column5=0;

$limits[0] = $min;
$limits[1] = $max;

if ( not $times ) {
  $times = 1.0;
}

if ( $yvals ) {
  print "Give the column for y-values \n";
  chop ($column5 = <STDIN>);
  $column5--;
  if ( $proces ) {
     print "You have chosen to get yvalues from divison of two columns\n";
     print "Give second column  (value taken is col2/col1)  \n";
     chop ($column3 = <STDIN>);
     $column3--;
#  print "Give fourth column\n";
#  chop ($column4 = <STDIN>);
#  $column4--;
  }

}


if ( $sub ) {
  print "Give second column\n";
  chop ($column2 = <STDIN>);
  $column2--;
}



if ( $dsample ) {
  print "Give sample limiting column\n";
  chop ($column3 = <STDIN>);
  $column3--;
  print "Give lower limit for sample\n";
  chop ($samplelow = <STDIN>);
  print "Give upper limit for sample\n";
  chop ($sampleup = <STDIN>);
}

if ( $sub ) {
  print "Write - to subtract, + to add, x to multiply or / to divide\n";
  chop ($whattodo = <STDIN>);
  print "Column ", $column+1 ," and ",$column2+1," are processed before histogram calculation in this order\n";
}

# ============================================================================
# BASIC INPUT OPERATIONS ARE DONE => OPEN FILES
# ============================================================================

open(INLIST, $inputfile)  || die "\n $0 Cannot open $! \n";
open(OUTLIST1, ">$inputfile".$user.".hist")  || die "\n $0 Cannot open $! \n";
$i=0;
$outlow=0;
$outup=0;
$histogram = undef;
$nofsteps = ($limits[1]-$limits[0])/$step;

for ($j=0;$j<$nofsteps;$j++)
  {
    $histogram[$j] = 0;
    $xvalues[$j] = $limits[0]+$j*$step+($step/2)
  }
print "Number of steps is $nofsteps\n";

# ============================================================================
# HISTOGRAM CALCULATION BEGINS
# ============================================================================

LINE: while (<INLIST>)
  {
    chomp;
    next LINE if (/^$/); # ignore null lines
    next LINE if (/^\s*#/); # ignore comment lines

    @data = split(); # Extract line data

    $histdata3 = $data[$column3];
    $histdata4 = $data[$column4];
    $ydata = $data[$column5];
    if($proces  && $yvals) {
       $ydata = $data[$column3]/$data[$column5];
       print "$ydata  $data[$column3] $data[$column5]\n";
    }
    if ( $visual ) {
       print "$data[$column] $data[$column2] $data[$column3] $data[$column4] \n";
    } 

    if ( $dsample ) {
      next LINE if ($histdata3 < $samplelow || $histdata3 > $sampleup)
    }

#    if ( $sub || $proces ) {
    if ( $sub ) {
	$histdata = @data[$column];
	$histdata2 = @data[$column2];
      if ( $sub ) {

         if ( $whattodo eq '-' ) {
            $histdata = abs($histdata-$histdata2);
         }
         if ( $whattodo eq '+' ) {
            $histdata = abs($histdata+$histdata2);
         }         
         if ( $whattodo eq 'x' ) {
            $histdata = abs($histdata*$histdata2);
         }
         if ( $whattodo  eq'/' ) {
            if ( $histdata2 == 0.0 ) {
               $histdata=0.0;
            } else
            {
              $histdata = abs($histdata/$histdata2);
            }
            print "$histdata\n";
         }         
       if ( $log ) {
	 $histdata = log10($histdata);
       }
     }

# This is just experimental part
#      if ( $proces ) {
#        if ( $histdata4 < $histdata3 && $histdata4 > -9000.0 ) {
#          print "$histdata3 $histdata4\n";
#          $histdata = $histdata2;
#        }
#      }
# Change previous if needed

    } else {
      if ( $log ) {
	$histdata = log10(@data[$column]);
      } else {
	$histdata = @data[$column];
      } 
    }

#      print "$data[$column3]\n";
	$histcolumn = ($histdata-$limits[0])/$step;
    if ($histcolumn < 0){$outlow++;}
    if ($histcolumn >= $nofsteps){$outup++;}
    if ($histcolumn >= 0 && $histcolumn < $nofsteps){ 
      if ( $yvals ) {
      $histogram[$histcolumn] = $histogram[$histcolumn] + $ydata;
      $histogramn[$histcolumn] = $histogramn[$histcolumn] + 1;
      $ysum = $ysum + $ydata;
      } else
      {
          $histogram[$histcolumn] = $histogram[$histcolumn] + 1;
      }
      $valsum = $valsum + @data[$column];
      $i++;
    }
    next LINE;
  }
if ($i == 0) {die "CHANGE LIMITS: NO VALUES INSIDE LIMITS\n";}

# ============================================================================
# HISTOGRAM IS CALCULATED => PERFORM OUTPUT OPERATIONS
# ============================================================================

if ( $yvals ) {
   $norm = $ysum;
} else
{
   $norm = $i;
}

# MAIN OUTPUT
for ($j=0;$j<$nofsteps;$j++)
  {
    if ( $log ) {
       $llim = 10**($limits[0]+$j*$step);
       $ulim = 10**($limits[0]+$j*$step+($step/2));
       $del10 = $ulim - $llim;
       if( $histogramn[$j] > 0 ) {
	   $frac = $histogram[$j]/$histogramn[$j];
       } else {
	   $frac = 0;
       }
       if( $yvals ) {
	   print OUTLIST1 $limits[0]+$j*$step+($step/2)."  ".$histogram[$j]*$times/$step.
	       "  ".$histogram[$j]*$times/$del10."  ".$histogram[$j]*$times/$norm.
	       "  ".10**($limits[0]+$j*$step+($step/2))."  ".$histogram[$j]."  ".$frac."\n";       
       }
       else {
	   print OUTLIST1 $limits[0]+$j*$step+($step/2)."  ".$histogram[$j]*$times/$step.
	       "  ".$histogram[$j]*$times/$del10."  ".$histogram[$j]*$times/$norm.
	       "  ".10**($limits[0]+$j*$step+($step/2))."  ".$histogram[$j]."\n";      
       }
   }
    else {
	if ( $yvals ) {
	    print OUTLIST1 $limits[0]+$j*$step+($step/2)."  ".$histogram[$j]*$times.
		"  ".$histogram[$j]*$times/$step."  ".$histogram[$j]*$times/$norm."  ".$histogram[$j]."  ".$frac."\n";
	} 
	else {
	    print OUTLIST1 $limits[0]+$j*$step+($step/2)."  ".$histogram[$j]*$times.
		"  ".$histogram[$j]*$times/$step."  ".$histogram[$j]*$times/$norm."  ".$histogram[$j]."\n";
	}
    }
  }
print "N of records in histogram $i\n";
print "N of records outside lower limits $outlow\n";
print "N of records outside upper limits $outup\n";
print "Sum of all values in the histogram: $valsum \n";
if ( $yvals ) {
   print "Total sum of y-values (chosen now from another column): $norm \n";
}
# ================================================================================
# IF CUMULATIVE DISTRIBUTION => CALCULATE CUMULATIVE DISTRIBUTION OF THE HISTOGRAM
# ================================================================================

if ($cumul) {
  open(OUTLIST3, ">$inputfile".$user.".cumhist")  || die "\n $0 Cannot open $! \n";
  for ($j=0;$j<$nofsteps+1;$j++)
  {
      if ( $log ) {
	  $llim = 10**($limits[1]-$j*$step);
	  $ulim = 10**($limits[1]-$j*$step-($step/2));
	  $del10 = $ulim - $llim;
	  $mid = ($ulim + $llim)/2.0;
      }
      else
      {
	  $llim = $limits[1]-$j*$step;
	  $ulim = $limits[1]-$j*$step-($step/2);
	  $mid = ($ulim + $llim)/2.0;
      }
      $cumhistogram[$nofsteps-$j] = $cumhistogram[$nofsteps-$j+1] + $histogram[$nofsteps-$j];
      $cumhistogram2[$nofsteps-$j] = $cumhistogram2[$nofsteps-$j+1] + $histogram[$nofsteps-$j]*$mid;
  }
  
  for ($j=0;$j<$nofsteps;$j++)
  {
      if ( $log ) {   
	  print OUTLIST3 $limits[0]+$j*$step+($step/2)."  ".$cumhistogram[$j]*$times."  ".
	      $cumhistogram2[$j]*$times."  ".$cumhistogram2[$j]*$times/$i."  ".
	      10**($limits[0]+$j*$step+($step/2))."\n";
      }
      else
      {	  
	  print OUTLIST3 $limits[0]+$j*$step+($step/2)."  ".$cumhistogram[$j]*$times."  ".
	      $cumhistogram2[$j]*$times."  ".$cumhistogram2[$j]*$times/$i."\n";
      }
  }

}

if ($visual) {
   &plot_data (\@xvalues, \@histogram );
}

            ########################
            # SUBROUTINES ARE BELOW
            ########################

sub min_max {
		my($max) = pop(@_);
		my($min) = pop(@_);
		foreach $foo (@_) {
			$max = $foo if $max < $foo;
			$min = $foo if $min > $foo;
		}
		return ($min,$max);
	}

sub plot_data {
    my(@xdata) = @{$_[0]};
    my(@ydata) = @{$_[1]};
    use strict;
    use Tk;
  my($o, $s) = (400, 20);
  my($pi, $x, $y, $xmaxv, $xminv, $ymaxv, $yminv, $xdelta, 
     $ydelta, $value, $i, $j, $x1, $y1) = (3.1415926, 0);
  my $mw = MainWindow->new;
  my $c = $mw->Canvas(-width => 500, -height => 500);
  my $padding = 10;
  $c->pack;
  $c->createLine(50, 450, 450, 450);
  $c->createText(250, 470, -fill => 'blue', -text => 'X');
  $c->createLine(50, 50, 50, 450);
  $c->createText(25, 250, -fill => 'blue', -text => 'Y');
  
# Configure the side pane:
  my $pane = $mw->Frame;
  $pane->pack( -side=>"left", -padx=>$padding, -pady=>$padding, -fill=>"both" );
  $pane->Button( -text=>"Exit", -command=>sub{ exit } )->pack( -side=>"bottom",
							       -pady=>$padding );

  ($xminv, $xmaxv) =  &min_max(@xdata);
  ($yminv, $ymaxv) =  &min_max(@ydata);
    $xdelta = $xmaxv - $xminv;
    $ydelta = $ymaxv - $yminv;
    print "xmin: $xminv xmax: $xmaxv\n";
    print "ymin: $yminv ymax: $ymaxv\n";

    $i = 0;
    $j = 0;
    foreach $value (@xdata)
    {
	$x = ($value - $xminv)/$xdelta;
	$y = (@ydata[$i] - $yminv)/$ydelta;
	$c->createText( $x*$o+50, 450-$y*$o, -fill => 'red', -text => 'x');
	if ($j > 1) {
	    $c->createLine($x1*$o+50, 450-$y1*$o, $x*$o+50, 450-$y*$o);
	} 
	$j++;
	$x1 = $x;
	$y1 = $y;
	$i++;
    }

  MainLoop;
}



__END__
#################################  IGNORE ALL BELOW THIS ############################

