docket_compare

#!/bin/env perl
$|=1;
use strict;
use JSON;
my(@files) = @ARGV;
my $kg = "knowledgegraphs";
my $scaling_factor = 100;

# Collect dockets to be compared
my @dockets;
my %graphFiles;
foreach my $file (@files) {
	if ($file =~ /^@/) {
		open F, $file;
		while (<F>) {
			chomp;
			my($d) = split /\t/;
			if (looks_like_a_docket($d)) {
				push @dockets, $d;
				foreach (fulldirlist("$d/$kg")) {
					$graphFiles{$1}{$d}++ if /(.+)\.json(\.gz)?$/;
				}
			} else {
				print "#doesn't look like a docket: $d\n";
			}
		}
		close F;
	} else {
		if (looks_like_a_docket($file)) {
			push @dockets, $file;
			foreach (fulldirlist("$file/$kg")) {
				$graphFiles{$1}{$file}++ if /(.+)\.json(\.gz)?$/;
			}
		} else {
			print "#doesn't look like a docket: $file\n";
		}
	}
}
die "Not enough dockets to compare\n" unless scalar @dockets>1;
@dockets = sort @dockets;

# Combine, model and report
print join("\t", qw/A B rho_avg confidence dockets minN maxN rho_slope rho_int pval_slope pval_int/), "\n";
foreach my $graph (sort keys %graphFiles) {
	my @docketsWithGraph = sort keys %{$graphFiles{$graph}};
	next unless scalar @docketsWithGraph>1;
	my %info;
	my %stat;
	foreach my $d (@docketsWithGraph) {
		my $json = read_json("$d/$kg/$graph.json.gz");
		foreach my $key (keys %$json) {
			if ($key eq 'tests') {
				# ingest all test results
				foreach my $blob (@{$json->{$key}}) {
					my $pair = join("\t", $blob->{'A'}, $blob->{'B'});
					next if $pair =~ /days/;
					foreach my $stat (keys %$blob) {
						next if $stat eq 'A' || $stat eq 'B';
						push @{$stat{$pair}{$stat}}, $stat eq 'pval' ? log(1e-100+$blob->{$stat}) : $blob->{$stat};
						#$stat{$pair}{$stat}{$d} = $blob->{$stat};
					}
				}
			} else {
				$info{$key}{$json->{$key}}++;
			}
		}
	}
	
	foreach my $pair (sort keys %stat) {
		next unless scalar @{$stat{$pair}{'N'}}>2;
		my($minN, $maxN) = minmax($stat{$pair}{'N'});
		my($m1, $b1, $avg1) = lsr_cols($stat{$pair}{'N'}, $stat{$pair}{'rho'});
		unless (defined $m1) {
			print "#weird pair $pair\n";
			next;
		}
		
		my($m2, $b2, $avg2) = lsr_cols($stat{$pair}{'N'}, $stat{$pair}{'pval'});
		if ($avg1<0) {
			$m1 = -$m1;
			$b1 = -$b1;
		}
		
		my $ksm = $m1*$scaling_factor - $m2;
		print join("\t", $pair, sprintf("%.3f",$avg1), sprintf("%.3f", $ksm),
			scalar @{$stat{$pair}{'N'}}, $minN, $maxN,
			sprintf("%.3f", $m1), sprintf("%.3f", $b1),
			sprintf("%.3f", -$m2), sprintf("%.3f", -$b2),
			), "\n";
		
	}
}

#########################
sub read_json {
	my($file) = @_;
	my $json;
	open J, "gunzip -c $file |";
	while (<J>) {
		chomp;
		$json .= $_;
	}
	close J;
	return decode_json($json);
}


sub looks_like_a_docket {
	my($dir) = @_;
	
	return 1 if
		-d $dir && 
		-d "$dir/analyses" && 
		-d "$dir/fingerprints" && 
		-d "$dir/comparisons" && 
		-d "$dir/visualizations";
}

sub avgstd {
	my($values) = @_;
	my($sum, $devsqsum);

	my $n = scalar @$values;
	return unless $n>1;
	foreach (@$values) { $sum += $_ }
	my $avg = $sum / $n;
	foreach (@$values) { $devsqsum += ($_-$avg)**2 }
	my $std = sqrt($devsqsum/($n-1));
	return $avg, $std;
}

sub minmax {
	my($values) = @_;
	my($min, $max);
	foreach (@$values) {
		$max = $_ if $_>$max;
		$min = $_ if $_<$min || !defined $min;
	}
	return $min, $max;
}


sub fulldirlist {
	my($dir) = @_;
	opendir (DIR, $dir);
	my @files = grep /^[^.]/, readdir DIR;
	closedir DIR;
	return @files;
}

sub lsr_cols {
	my($xdata, $ydata) = @_;
	
	my($n, $sx, $sy, $sxy, $sxsq, $x, $y);
	foreach my $i (0..$#$xdata) {
		$x = $xdata->[$i];
		$y = $ydata->[$i];
		$n++;
		$sx += $x;
		$sy += $y;
		$sxy += $x*$y;
		$sxsq += $x**2;
	}
	return unless $n*$sxsq-$sx**2;
	my $m = ($n*$sxy-$sx*$sy)/($n*$sxsq-$sx**2);
	my $b = ($sy-$m*$sx)/$n;
	return $m, $b, $sy/$n;
}