$left_idx++ if $store->[$left_idx] < $start; 

##</code><code>##

#  $left_idx++ if $store->[$left_idx] < $start; #normal
  $left_idx++ if ($store->[$left_idx]+1) < ($start+1); #0.6 bug workaround

##</code><code>##

1.6     [8]:    0.8     1       1.2     1.4     1.8     2       2.2     2.4     2.6
2.2     [11]:   1.4     1.6     1.8     2       2.4     2.6     2.8     3       3.2
4.4     [22]:   3.6     3.8     4       4.2     4.6     4.8     5       5.2     5.4
8.8     [44]:   8       8.2     8.4     8.6     9       9.2     9.4     9.6     9.8
16.6    [83]:   15.8    16      16.2    16.4    16.8    17      17.2    17.4    17.6
32.2    [161]:  31.4    31.6    31.8    32      32.4    32.6    32.8    33      33.2
64.4    [322]:  63.6    63.8    64      64.2    64.6    64.8    65      65.2    65.4

##</code><code>##

  DB<25> x ($store->[3]) == $start
0  ''
  DB<26> x ($store->[3]+0) == $start
0  ''
  DB<27> x $start
0  0.6
  DB<28> x $store->[3]
0  0.6
  DB<29> x 0.6 == $start
0  ''
  DB<30> x $start == 0.6
0  ''
  DB<31> x 0.6 == 0.6
0  1
  DB<32> x ($store->[3]+0) == ($start+0)
0  ''
  DB<33> x ($store->[3]+1) == ($start+1)
0  1

##</code><code>##

#call:
use Spikehisto;
$histo=autocorr([map{$_/5}(0..1000)],2,100,2000); 
print join "\t", grep {$_} @$histo;

##</code><code>##

#!/usr/bin/perl -w

package Spikehisto;

use strict;
use warnings;
use diagnostics;

use Math::Round;

use Exporter;
use vars qw( @EXPORT @ISA @VERSION);
@VERSION = 0.01;
@ISA = qw( Exporter );
@EXPORT = qw( &autocorr &crosscorr &interval);


=head1 $a = crosscorr( $source, $compare, $interval, $number_bins, $max_count)


Calculates crosscorrelation histogram of $compare (arrayref) to the source
(arrayref). The $interval gives the interval around the source event about
which the crosscorrelation is calculated. The $num_bins denotes the number
of bins of the interval. If the number of bins is odd, the boundaries will
be lower. The $max_count gives the maximum number of spikes to count.


=cut


my $PDL = <<'_PDL';
_crosscorr_
extract arguments
prepare
return empty histogram unless the source has more than two items
skip the items at the beginning of the source for a length of half an interval
skip the items at the end of the source for a length of half an interval
walk over the remaining items
  fill the buffer
  if the buffer has one or more items
    add the buffer to the histogram
    increase the count
  exit if the count reaches the desired count
return the histogram
_PDL

sub crosscorr {
  #extract arguments
  my $source = shift;
  my $comp = shift;
  my $cross_interval = shift;
  my $cross_bins = shift;
  my $count_crosscorr = shift;
  
  my $count = 0;
  my @times_source = @$source;
  my @times_comp = @$comp;
  my $start_time = $times_source[0];
  my $end_time = $times_source[ $#times_source ] - $cross_bins/2;
  my @cross_histo = (0) x $cross_bins;

  return \@cross_histo unless scalar(@times_source)>1;
  
  my $idx = 1;
  $idx++ while( $times_source[$idx] < ($start_time + $cross_interval/2) 
            and $idx < $#times_source );
  
  my $last_idx = $#times_source;
  $last_idx-- while( $times_source[$last_idx] > $end_time and $last_idx > 0 );
  
  my ($comp_begin_idx, $comp_last_idx) = (0,0);
  while( $idx <= $last_idx ) {
    my $buff = _fill_crossbuffer( \@times_comp, $times_source[$idx], $cross_interval,
       \$comp_begin_idx, \$comp_last_idx);
    if (scalar @$buff){
      _buffer2histo( $buff, \@cross_histo, $times_source[$idx] - $cross_interval/2, 
        $cross_interval, $cross_bins );
      $count++;
    }
    last if $count > $count_crosscorr;
    $idx++;
  }
  
  \@cross_histo;
}
 

=head1 $a = autocorr( $source, $interval, $number_bins, $max_count)


Calculates autocorrelation histogram of $source (arrayref). The $interval
gives the interval around the source event about which the autocorrelation
is calculated. The $num_bins denotes the number of bins of the interval. If
the nubmer of bins is odd, the boundaries will be lower. The $max_count
gives the maximum number of spikes to count.


=cut



$PDL = <<'_PDL';
_autocorr_
extract arguments
prepare
return empty histogram unless the source has more than two items
skip the items at the beginning of the source for a length of half an interval
skip the items at the start of the source for a length of half an interval
walk over the remaining times
  fill the buffer with an interval around the current spike
  if the buffer has one or more items
    add the buffer to the histogram
    increase the count
  exit if the count reaches the desired count
return the histogram
_PDL


sub autocorr {
  my $source = shift;
  my $auto_interval = shift;
  my $auto_bins = shift;
  my $count_autocorr = shift;

  my $count = 0;
  my @times = @$source;
  my $start_time = $times[0];
  my @auto_histo = (0) x $auto_bins;
  
  return \@auto_histo unless scalar(@times)>1;
  
  my $idx = 1;
  $idx++ while( $times[$idx] < ($start_time + $auto_interval/2) and $idx < $#times );
  my $last_idx = $#times;
  my $end_time = $times[ $#times ] - $auto_bins/2;
  $last_idx-- while( $times[$last_idx] > $end_time and $last_idx > 0 );
  
  while( $idx <= $last_idx ) {
    my $buff = _fill_autobuffer( \@times, $idx, $auto_interval );
    if (scalar @$buff){
      _buffer2histo( $buff, \@auto_histo, $times[$idx] - $auto_interval/2, 
        $auto_interval, $auto_bins );
      $count++;
    }
    last if $count > $count_autocorr;
    $idx++;
  }
  
  \@auto_histo;
}


=head1 $a = interval( $source, $interval, $number_bins)


Calculates inter-event interval histogram of the source (arrayref). The
$interval gives the maximum inter-event interval that should be counted. The
$number_bins denotes the number of bins of the interval. 


=cut


$PDL = <<'_PDL';
_interval_
extract arguments
prepare
return empty histogram unless the source has more than two items
get the first item of the series as a reference
walk over series
  calculate difference between item and reference
  assign current item to reference
add the buffer to the histogram
return the histogram
_PDL


sub interval {
  my $source = shift;
  my $interval = shift;
  my $interval_bins = shift;
  
  my @times = @$source;
  my $start_time = $times[0];
  my @interval_histo = (0) x $interval_bins;
  
  return \@interval_histo unless scalar(@times)>1;
  
  my $prev = shift @times;
  my @buff = map { my $time = $_ - $prev; $prev = $_; $time} @times;
  
  buffer2histo( \@buff, \@interval_histo, 0, $interval, $interval_bins ) if scalar @buff;
  
  \@interval_histo;
}


$PDL = <<'_PDL';
__buffer2histo_
extract arguments
walk over buffer
  calculate time relative to reference spike
  add time to histogram
_PDL


sub _buffer2histo{
  my $buff = shift;
  my $histo = shift;
  my $reftime = shift;
  my $interval = shift;
  my $length = shift;
  
  for (@$buff){
    my $time = $_ - $reftime;
    my $idx = round( $time *  ($length-1) / $interval );
    next if $idx >= $length or $idx < 0;
    $histo->[$idx]++;
  }
}

 
$PDL = <<'_PDL';
__fill_autobuffer_
extract arguments
determine start and end borders, a half interval before and after the reference spike
determine the index of start and end
assign the slice to the buffer, leaving out the reference spike
return the buffer
_PDL

sub _fill_autobuffer{
  my $store = shift;
  my $store_idx = shift;
  my $interval = shift;
  
  my @buff;
  my $left_idx = $store_idx;
  my $start = $store->[$store_idx] - $interval/2;
  my $end = $store->[$store_idx] + $interval/2;
  
  $left_idx-- while( $store->[$left_idx] > $start and $left_idx > 0);
  $left_idx++ if $store->[$left_idx] < $start; #normal
#  $left_idx++ if ($store->[$left_idx]+1) < ($start+1); #0.6 bug workaround
  my $right_idx = $store_idx;
  $right_idx++ while( $#{$store} > $right_idx and $store->[$right_idx] < $end );
  $right_idx-- if not defined $store->[$right_idx] or $store->[$right_idx] > $end;
  
  my $idx = $store_idx;
  @buff = @{$store}[ 
    $left_idx .. ($store_idx - 1), ($store_idx + 1) .. $right_idx 
  ];
  
  \@buff;
}


$PDL = <<'_PDL';
__fill_crossbuffer_
extract arguments
determine start and end borders, a half interval before and after the reference spike
determine the index of start and end
assign the slice to the buffer
point the index-references to the new start and end
return the buffer
_PDL

sub _fill_crossbuffer{
  my $store = shift;
  my $time_source = shift;
  my $interval = shift;
  my $left_idx_ref = shift;
  my $right_idx_ref = shift;
  
  my $left_idx = $$left_idx_ref;
  my $right_idx = $$right_idx_ref;
  my @buff;
  my $start = $time_source - $interval/2;
  my $end = $time_source + $interval/2;
  
#  $left_idx++ while( $store->[$left_idx] < $start and $left_idx < $#{$store}); #normal
  $left_idx++ while( ($store->[$left_idx]+1) < ($start+1) and $left_idx < $#{$store}); #bug workaround
  $right_idx++ while( $#{$store} > $right_idx and $store->[$right_idx] < $end );
  $right_idx-- if not defined $store->[$right_idx] or $store->[$right_idx] > $end;
  
  @buff = @{$store}[ $left_idx .. $right_idx ] ;
  
  $left_idx_ref = \$left_idx;
  $right_idx_ref = \$right_idx;
  
  \@buff;
}

1;