#!/usr/local/bin/perl5 -s
 
require 5.000;

open (MARKERS, "<clustered_markers.txt") or die "cannot open MARKERS\n";
open (RSCAN, ">r-scan.txt") or die "cannot open R-SCAN\n";

@r_values = (1, 2, 3, 4, 5, 7, 10);
$factorial[0] = 1;
for ($i = 1; $i < 11; $i++) {
$factorial[$i] = $factorial[$i-1] * $i;
}

##	parse markers
while (<MARKERS>) {
if (/(\d+)/) {
	$new_marker = $1;
	push @markers, $new_marker;
}
}

		# select unique markers
		%seen = ();
		foreach $item (@markers) {
    			$seen{$item}++;
			}
		@markers = keys %seen;

$n = @markers;
@ordered_markers = sort {$a <=> $b} @markers;
$big_n = $ordered_markers[$n-1] - $ordered_markers[0];

print RSCAN "n: \t$n\nN: \t$big_n\nr\tk\tmax/min\tUr\tx\t\tp\n";

foreach $r (@r_values) {
@u = ();

    for ($i = 0; $i < ($n - $r); $i++) {
	$big_r = $ordered_markers[$i+$r] - $ordered_markers[$i];
	push @u, $big_r;
	}

    @ordered_u = sort {$a <=> $b} @u;
    @reversed_u = reverse @ordered_u;

$k_value = 10;
    for ($k = 1; $k <= $k_value; $k++) {

	$x = $ordered_u[$k] * ($n ** (1.0 + (1.0/$r))) / $big_n;
	$pr_small_m = &prob_lt($x, $r, $k);
	#if ($pr_small_m < 0.05) {
	print RSCAN "$r\t$k\tmin\t$ordered_u[$k]\t$x\t$pr_small_m\n";
	#}

	$x = ($reversed_u[$k] * $n - log ($n) - ($r - 1.0) * log ( log ($n))) / $big_n;
	$pr_big_m = &prob_gt($x, $r, $k);
	#if ($pr_small_m < 0.05) {
	print RSCAN "$r\t$k\tmax\t$reversed_u[$k]\t$x\t$pr_big_m\n";
	#}

	}
} #r_values loop

close (MARKERS) || die "MARKERS file was not closed\n";
close (RSCAN) || die "RSCAN file was not closed\n";
exit (0);

sub prob_lt($x, $r, $k)
{
#print RSCAN "in prob_lt:\tr-fragment= $x, r= $r, k= $k\n";
$lambda = ($x ** $r)/$factorial[$r];
$sum = 0;
for  ($i = 0; $i < $k; $i++) {
$sum = $sum + (($lambda ** $i) / $factorial[$i]) * exp (-1.0 * $lambda);
}
$prob = 1.0 - $sum;
return $prob;
}

sub prob_gt($x, $r, $k)
{
#print RSCAN "in prob_gt:\tr-fragment= $x, r= $r, k= $k\n";
$mu = (exp (-1 * $x)) / $factorial[$r-1];
$sum = 0;
for  ($i = 0; $i < $k; $i++) {
$sum = $sum + (($mu ** $i) / $factorial[$i]) * exp (-1.0 * $mu);
}
$prob = 1.0 - $sum;
return $prob;
}

__END__

	#print RSCAN "m(r)k = $ordered_u[$k]\n";
	#print RSCAN "Pr{m(r)k < $ordered_u[$k]} = $pr_small_m\n";
	#print RSCAN "M(r)k = $reversed_u[$k]\n";
	#print RSCAN "Pr{M(r)k > $reversed_u[$k]} = $pr_big_m\n";