#!/usr/bin/perl
# commandline: hp_hist file_name Ts_min Ts_max
# reads file_bame.scatter, which has
# first the histogram of the RR-intervals in 2 columns,
# and after that, in four columns:
# 1. coupling interval[s]
# 2. V-V time interval[s]
# 3. NIBs
# 4. corresponding RR-interval[s]
#     which might be averaged over some beats 
# The input file may be generated with hp_scatter.c
#############################################################

unless (@ARGV){
 die "\n usage:  hp_hist  file_name Ts_min Ts_max binsize\n 
      (Ts_max-Ts_min)/binsize = integer !!! \n\n";
} 

$filename=shift(@ARGV);
$min=shift(@ARGV);
$max=shift(@ARGV);
$binsize=shift(@ARGV);

if ($binsize <= 0) {
    $binsize = 1;	# avoid dividing by zero below
}

$Ts_steps=&round(($max-$min)/$binsize);

########### Tv #############
  $max_Tv=10;
  $bintime=10*$binsize;
  $Tv_steps=&round($max_Tv/$bintime);

for ($t=0;$t<=$Tv_steps;$t++) {
    for ($i=0;$i<=$Ts_steps;$i++) {
	$rrTv[$t][$i]=0;
    }
}

########### NIB #############
$max_nib=16;
for ($t=0;$t<=$max_nib;$t++) {
    for ($i=0;$i<=$Ts_steps;$i++) {
	$rrnib[$t][$i]=0;
    }
}

########### coup #############
  $max_coup=1.25;
  $min_coup=0.25;
  $coup_steps=&round(($max_coup-$min_coup)/$binsize);
for ($t=0;$t<=$coup_steps;$t++){
    for ($i=0;$i<=$Ts_steps;$i++)
	{
	    $rrcoup[$t][$i]=0;
	}
}
#############################


unless(open(file1,"$filename")){
    die "Cannot open file: $!";
}

while (<file1>) {
    @F=split;
    if (@F==2) 	{
	$hist_rr[&round(($F[0]-$min)/$binsize)]=$F[1];
    }
    elsif (($F[3]>=$min)&&($F[3]<=$max))  {
	$sinusrr=&round(($F[3]-$min)/$binsize);	
	if ($F[1]<=$max_Tv+$binsize/2) {
	    $t=&round($F[1]/$bintime);
	    $rrTv[$t][$sinusrr]++;
            $k=$rrTv[$t][$sinusrr];
	}
	if ($F[2]<=$max_nib){
	    $t=$F[2];
	    $rrnib[$t][$sinusrr]++;
	}
	if (($F[0]<=$max_coup)&&($F[0]>=$min_coup)) {
	    $t=&round(($F[0]-$min_coup)/$binsize);
	    $rrcoup[$t][$sinusrr]++;
	}
    }
}
close(file1);


$file=$filename.".hist_Tv";
&printhist($file, $Ts_steps, $Tv_steps, \@rrTv, 1); 

$file=$filename.".rrTv";
&printrr($file, $Ts_steps, $Tv_steps, \@rrTv, \@hist_rr, 10); 


$file=$filename.".hist_nib";
&printhist($file, $Ts_steps, $max_nib, \@rrnib, 2); 

$file=$filename.".rrnib";
&printrr($file, $Ts_steps, $max_nib, \@rrnib, \@hist_rr, 0); 


$file=$filename.".hist_coup";
&printhist($file, $Ts_steps, $coup_steps, \@rrcoup, 3); 

$file=$filename.".rrcoup";
&printrr($file, $Ts_steps, $coup_steps, \@rrcoup, \@hist_rr, 10); 


$file=$filename.".hist_rr";
unless(open(file1,">$file")){
    die "Cannot open file: $!";
}
for ($i=0;$i<=$Ts_steps;$i++){
    printf (file1 "%f %d \n", $min+$i*$binsize, $hist_rr[$i]);            
}
close(file1);


#################################
########## subroutines ##########
#################################
sub round {
    my $x= shift(@_);
    return  sprintf("%.0f",$x);
}

#################################
# prints greyscale plot
# call as printrr(filename, Ts_steps, t_steps, ref_rr_x, ref_rr_hist) 

sub printrr {
    my $file = shift(@_);
    my $Ts_steps = shift( @_);
    my $t_steps  = shift( @_);
    my ($ref_rr_x) = shift( @_);
    my ($ref_rr_hist) = shift( @_);
    my $shift = shift( @_);
    my $maximum = 0;

    unless(open(file1,">$file")){
	die "Cannot open file: $!";
    }
    print file1 $Ts_steps," ", $t_steps, "\n";

    for ($t=0;$t<$t_steps;$t++)
    {
	for ($i=0;$i<$Ts_steps;$i++){
	    $auxRRhist = ($ref_rr_hist->[$i]);
	    $aux=$ref_rr_x->[$t][$i];
	    if ($auxRRhist > $Ts_steps/1000) {  # was "> 100"
		$aux /= $auxRRhist;
		if (($t > $shift)&&($aux > $maximum)){
		    $maximum = $aux;
		}
	    }
	    else {
		$aux = 0;      
	    }
	    $ref_rr_x->[$t][$i] = $aux;
	}
    }

    if ($maximum <= 0) {
	$maximum = 1;	# avoid dividing by zero below
    }
    print $maximum," \n";

    for ($t=0;$t<$t_steps;$t++)
    {
	for ($i=0;$i<$Ts_steps;$i++){
	    $aux=($ref_rr_x->[$t][$i]/$maximum);
	    print file1 $aux," ";
	}
	printf ( file1 "\n");
    }
    close(file1);
}

#################################
# prints histogram, call as:
# printhist(filename, Ts_steps, t_steps, ref_rr_x, flag=1 Tv, =2: nib =3:coup) 

sub printhist {
    my $file = shift(@_);
    my $a_steps = shift( @_);
    my $t_steps  = shift( @_);
    my ($ref_rr_x) = shift( @_);
    my $flag       = shift(@_);

    unless(open(file1,">$file")){
	die "Cannot open file: $!";
    }
    for ($t=0;$t<=$t_steps;$t++)
	{
	    $x=0;
	    for ($i=0;$i<=$a_steps;$i++){
		$x+=$ref_rr_x->[$t][$i];
	    }
	    if ($flag==1){
		printf (file1 "%f %d \n", $t*$bintime, $x);            
	    }
	    elsif ($flag==2){
		printf (file1 "%f %d \n", $t, $x);            
	    }
	    else {
		printf (file1 "%f %d \n", $t*$binsize+$min_coup, $x);                     
	    }
	}
    close(file1);
}