#!/usr/bin/env perl

# spearman.perl - calculating spearman's rank correlation value
#                 from s1.txt and s2.txt
#
# September 29, 2006 TOMINAGA Daisuke
#
# Input: {s1,s2}.txt
# Output: stdout
#         tab separated text. the second column is the calculated result.
# ref. numerical receipes in fortran 77 second ed.

open(S1, "s1.txt"); @S1 = <S1>; close(S1); $s1 = @S1;
open(S2, "s2.txt"); @S2 = <S2>; close(S2); $s2 = @S2;
if ($s1 != $s2) { print "contens of s1.txt and s2.txt differ.\n"; exit; }

for ($i = 0; $i < @S1; $i++) {
	@elm = split(/\s+/, $S1[$i]); $ps1{$elm[0]} = $elm[1];
	@elm = split(/\s+/, $S2[$i]); $ps2{$elm[0]} = $elm[1];
}

@ss1 = sort {$ps1{$a} <=> $ps1{$b}} keys %ps1;
@ss2 = sort {$ps2{$a} <=> $ps2{$b}} keys %ps2;

$prep = $m = 0;
for ($i = 0; $i < @ss1; $i++) {
	if ($ps1{$ss1[$i]} == $prep) { $m++; }
	else {
		$rs1{$ss1[$i]} = $i;
		$prep = $ps1{$ss1[$i]};
		if ($m != 0) {
			$r = $i - 1 - $m / 2;
			for ($j = $i - $m - 1; $j < $i; $j++) { $rs1{$ss1[$j]} = $r; }
			@ts1 = (@ts1, $m); $m = 0;
		}
	}
}
$i = @ss1;
if ($m != 0) {
	$r = $i - 1 - $m / 2;
	@ts1 = (@ts1, $m);
	for ($j = $i - $m - 1; $j < $i; $j++) { $rs1{$ss1[$j]} = $r; }
}
# same process as the above, s1
$prep = $m = 0;
for ($i = 0; $i < @ss2; $i++) {
	if ($ps2{$ss2[$i]} == $prep) { $m++; }
	else {
		$rs2{$ss2[$i]} = $i;
		$prep = $ps2{$ss2[$i]};
		if ($m != 0) {
			$r = $i - 1 - $m / 2;
			for ($j = $i - $m - 1; $j < $i; $j++) { $rs2{$ss2[$j]} = $r; }
			@ts2 = (@ts2, $m); $m = 0;
		}
	}
}
if ($m != 0) {
	$r = $i - 1 - $m / 2;
	@ts2 = (@ts2, $m); 
	for ($j = $i - $m - 1; $j < $i; $j++) { $rs2{$ss2[$j]} = $r; }
}

$D = 0;
foreach $key (keys %rs2) { $D += ($rs1{$key} - $rs2{$key})**2; }

$N = @ss1; $Ni = $N**3 - $N;
$tx = $ty = 0;
for ($i = 0; $i < @ts1; $i++) { $tx += ($ts1[$i]**3 - $ts1[$i]); }
for ($i = 0; $i < @ts2; $i++) { $ty += ($ts2[$i]**3 - $ts2[$i]); }
$rho = -(2*(6*$D-$Ni)+$tx+$ty)/(2*sqrt(($Ni-$tx)*($Ni-$ty)));
print "spearman's rank correlation value:\t$rho\n";