#!/usr/bin/perl

# Calculate minimum significant level of agreement for real variation units.

# Copyleft, T. J. Finney, 2002-11-17. Free for non-commercial use.

# Set confidence level.
$c = 0.95;

# Set variants file name.
# Each line of this file must have the following format:
# <label> <TAB> <number> <TAB> <number> ... <TAB> <number>
$input_file = "variants.txt";

if (($c < 0) or (1 < $c)) {
    die "Probabilities must be between 0 and 1"; 
}

@weights = ();
# Parse variants file.
open(INPUT, $input_file) or die "Cannot open $input_file";
@list = <INPUT>;
@list = grep(/^\S/, @list);
foreach $item (@list) {
    chomp($item);
    @record = split(/\t/, $item);
    shift(@record);
    foreach $entry (@record) {
	unless ($entry =~ /\d+/) {
	    $s = join("\t", @record);
	    die "Format error in $input_file: $s\n";
	}
    }
    push(@weights, [@record]);
}

# Set number of trials.
$n = 1000;

# Create array of agreements.
srand();
@agreements = ();
# Create first MS.
@MS1 = &create_MS(@weights);
# Loop N times.
for ($i = 0; $i < $n; $i++) {
    # Create second MS.
    @MS2 = &create_MS(@weights);
    # Measure level of agreement with first MS.
    $level = &compare_MSS(@MS1, @MS2);
    # Add to array of agreements.
    push(@agreements, $level);
    # Create first MS.
    @MS1 = @MS2;
    # print join("", @MS1), "\n";
}

# Find minimum significant level of agreement.
@agreements = sort { $a <=> $b } @agreements;
$index = $c * $n - 1;
$value = $agreements[$index];
for ($i = $index + 1; $i < $n; $i++) {
    if ($agreements[$i] > $value) {
	last;
    }
}

# Print result.
printf("Minimum significant level of agreement (confidence level = %d\%):\n", $c * 100);
if ($i < $n) {
    $msl = $agreements[$i];
    $total = $#weights + 1;
    printf("%d/%d (%.1f\%)\n", $msl, $total, $msl / $total * 100);
}
else {
    print "none exists\n";
}

sub create_MS {
    my @weights = @_;
    my $n = $#weights + 1;
    my @record = ();
    my @MS = ();
    for (my $i = 0; $i < $n; $i++) {
	@record = @{@weights[$i]};
	push(@MS, &create_reading(@record));
    }
    return @MS;
}

sub create_reading {
    my @record = @_;
    my $n = $#record + 1;
    my @list;
    for (my $i = 0; $i < $n; $i++) {
	my $m = $record[$i];
	for (my $j = 0; $j < $m; $j++) {
	    push(@list, $i + 1);
	}
    }
    my $index = int(rand($#list + 1));
    return $list[$index];
}

sub compare_MSS {
    my @list = @_;
    my $n = int(($#list + 1) / 2);
    my $level = 0;
    for (my $i = 0; $i < $n; $i++) {
	if ($list[$i] == $list[$i + $n]) {
	    $level++;
	}
    }
    return $level;
}