#!/usr/bin/perl
#
#   coords_2_dist_array.pl <FName>
#
#   Given  a list of coordinates in file <FName>, compute the array of point-to-point
#   distances.  If FName isn't given, use DATA
#
use strict;
use warnings;
use autodie;

my $IFH;
my $FName = shift;
if (defined $FName) {
    open $IFH, '<', $FName;
}
else {
    $IFH = \*DATA;
}

my @PTs;
while (<$IFH>) {
    my @coords = split /\s+/,$_;
    push @PTs, [ @coords ] if @coords;
}

for my $i (0 .. $#PTs) {
    my @vals=();
    for my $j (0 .. $#PTs) {
        push @vals, sprintf "%.2f", dist($PTs[$i], $PTs[$j]);
    }
    print join("\t", @vals), "\n";
}

sub dist {
    my ($r1, $r2) = @_;
#    print ref($r1), " ", ref($r2),": ",join(",",@$r1), "; ", join(",",@$r2), "\n";
    die "MISMATCH! $#$r1, $#$r2\n" unless $#$r1 == $#$r2;
    my $sum = 0;
    for my $i (0 .. $#$r1) {
        $sum += ($$r1[$i]-$$r2[$i])**2;
    }
    return sqrt($sum);
}

__DATA__
0 0
300 -100
550 550
800 -200
1100 200
1250 0

##</code><code>##

#!/usr/bin/perl
#
#   trv_slsmn_triangulate.pl <FName>
#
#   Given the distance matrix, compute a set of 2D locations for the points that yields
#   the specified distance map.
#
#   Since any coordinates will do, we'll first set point A at (0,0), and point B at (aDb,0).
#   Finally, we'll make sure the angle of point C is counterclockwise from the line AB.
#   Obviously, we have to be sure that aDc+bDc > aDb, or it won't make a triangle, so we
#   have to use a bit of care to select our first points.
#
use strict;
use warnings;
use List::Util qw( shuffle );

# Prepare data source
my $FName = shift;
my $IFH;
if (defined $FName) {
    open $IFH, '<', $FName;
}
else {
    $IFH = \*DATA;
}

# Read data
my %D;
my ($R, $C) = (0);
my $max_dist = [ 0 ];
my @LOCs;
while (<$IFH>) {
    my @t = split /\s+/, $_;
    $C=0;
    $LOCs[$R] = [ ];
    while (@t) {
        my ($dist, $a, $b);
        $dist = shift @t;
        ($a,$b) = ($R>$C) ? ($C, $R) : ($R, $C);
        if ($a != $b) {
            if (exists $D{$a}{$b}) {
                die "Mismatch ($R,$C)=$dist, but ($C,$R)=$D{$a}{$b}\n" if $dist != $D{$a}{$b};
            }
            $D{$a}{$b} = $dist;
            if ($dist > $$max_dist[0]) {
                $max_dist = [ $dist, $a, $b ];
            }
        }
        ++$C;
    }
    ++$R;
}

# OK, points A and B are the points in max_dist, set their initial locations:
my ($ptA, $ptB) = ($max_dist->[1], $max_dist->[2]);
$LOCs[$ptA] = [0, 0];
$LOCs[$ptB] = [$max_dist->[0], 0];
print "Point A is $ptA (0, 0); point B is $ptB ($max_dist->[0], 0)\n";

# For point C, we'll find the point giving the triangle with the largest area
$max_dist = [ 0 ];
for my $i (1 .. $#LOCs) {
    next if defined $LOCs[$i][0]; # skip A and B
    my ($dCA, $dCB, $dAB) = (dist($i,$ptA), dist($i,$ptB), dist($ptA,$ptB));
    my $s = ($dCA+$dCB+$dAB)/2.0;
    #my $dsq = $dA*$dA + $dB*$dB;
    #my $dsq = abs($dA-$dB) / ($dA*$dA+$dB*$dB);
    my $dsq = sqrt( $s * ($s-$dCA) * ($s-$dCB) * ($s-$dAB) );
#    print "pt $i: dCA=$dCA, dCB=$dCB, dAB=$dAB, s=$s area=$dsq\n";
    if ($dsq > $max_dist->[0]) {
#        print "\tnew ptC candidate\n";
        $max_dist = [$dsq, $i];
    }
}
my $ptC = $max_dist->[1];
print "Point C is $ptC\n";

# Find position of c via law of cosines (CAB) and sin(arccos x)=sqrt(1-x**2)
my ($dCA, $dBA, $dCB) = (dist($ptC,$ptA), dist($ptB,$ptA), dist($ptC,$ptB));
print "dCA=$dCA, dBA=$dBA, dCB=$dCB\n";
my $cosCAB = ($dCA*$dCA + $dBA*$dBA - $dCB*$dCB) / (2*$dCA*$dBA);
my $sinCAB = sqrt(1 - $cosCAB**2);
$LOCs[$ptC] = [ $dCA * $sinCAB, $dCA * $cosCAB ];
print "\t($LOCs[$ptC][0], $LOCs[$ptC][1])\n";

# remaining locations can be found by triangulation from the three known points
for my $i (1 .. $#LOCs) {
    next if $i==$ptA or $i==$ptB or $i==$ptC;
    my ($dIA, $dIB, $dIC) = (dist($i,$ptA), dist($i,$ptB), dist($i,$ptC));
    print "$i: ($dIA, $dIB, $dIC)\n";
    my $cosIAB = ($dIA*$dIA + $dBA*$dBA - $dIB*$dIB) / (2*$dIA*$dBA);
    my $sinIAB = sqrt(1 - $cosIAB**2);
    my ($x, $y) = ($dIA*$cosIAB, $dIA*$sinIAB);
    my $check = sqrt( ($LOCs[$ptC][0]-$x)**2 + ($LOCs[$ptC][1]-$y)**2 );
    my $err = abs($check-$dIC);
    my $err_pct = 100*$err / $dIC;
    if ($err < $dIC*0.01) {
        print "\t($x, $y) <case 1: $err, $err_pct>\n";
        next;
    }
    my ($x2, $y2) = ($dIA*$cosIAB, -$dIA*$sinIAB);
    my $check2 = sqrt( ($LOCs[$ptC][0]-$x2)**2 + ($LOCs[$ptC][1]-$y2)**2 );
    $err = abs($check2-$dIC);
    $err_pct = 100*$err / $dIC;
    if ($err < $dIC*0.01) {
        print "\t($x2, $y2) <case 1: $err, $err_pct>\n";
        next;
    }
    print "\t($x,$y),  ($x2,$y2)\n";
    print "\tcheck=$check, check2=$check2, dIC=$dIC\n";
}

sub dist {
    my ($r, $c) = @_;
    die "Eh?  $r == $c!\n" if $r eq $c;
    ($r, $c) = ($c, $r) if $c < $r;
    die "$r,$c entry DNE?\n" unless exists $D{$r}{$c};
    $D{$r}{$c};
}

# Matrix shows distance from pt on left to dest column
__DATA__
0  3812 13902  8619 15811  5015  5230  9615 10624 13346  7575  6170  6812  9487 18135  8030  5409 17959 12822 17136  3267 12882 11223 11078
3812     0 11527 12431 15446  8057  4519  8761 14398 12569  3764  6668 10603 11117 14805  5154  8276 18175  9367 14840  7056  9698 13603  7266
13902 11527     0 13638 10220 18405  8675  4611 12993 11970  8798  8226 15087 16591  6859  6381 11223  7602  3236  3457 14535  8830 14748  6655
8619 12431 13638     0  9965  5555 11256 11549  2157 10917 16194  9609  1926  7111 12565 14906  5737  9378 16868 11899  5402 15921  5765 19675
15811 15446 10220  9965     0 11122 18683 14599  7873  2940 13119 17793 11057  6374  3517 14973 15565  3627 10307  7077 14478  6475  5155 10235
5015  8057 18405  5555 11122     0 10109 13856  6744  9617 11348 10270  3811  4858 14594 12975  7186 13212 17301 17106  3881 12658  6210 14138
5230  4519  8675 11256 18683 10109     0  4584 13284 16923  6018  2299 10267 14709 15205  3741  5549 15863  8377 11981  6865 12543 16177  8600
9615  8761  4611 11549 14599 13856  4584     0 12503 16548  8240  3619 11888 18524 11444  4627  6948 11282  6102  7523  9992 12321 16782  8550
10624 14398 12993  2157  7873  6744 13284 12503     0  9231 18070 11405  3813  6447 10427 16699  7735  7340 15982 10421  7492 14151  4374 18092
13346 12569 11970 10917  2940  9617 16923 16548  9231     0 11052 19220 11199  4874  5269 14093 16358  6543 10709  9458 13488  5008  5170  9146
7575  3764  8798 16194 13119 11348  6018  8240 18070 11052     0  8220 14354 12490 11186  3614 11306 14452  5954 11533 10801  6649 14987  3508
6170  6668  8226  9609 17793 10270  2299  3619 11405 19220  8220     0  9124 15115 15059  5330  3975 14165  9100 10994  6469 14453 15268 10326
6812 10603 15087  1926 11057  3811 10267 11888  3813 11199 14354  9124     0  6645 14126 13999  5168 11153 17968 13808  3728 15782  6168 17749
9487 11117 16591  7111  6374  4858 14709 18524  6447  4874 12490 15115  6645     0  9754 15942 11588  9288 15400 13331  8652  9157  2621 12782
18135 14805  6859 12565  3517 14594 15205 11444 10427  5269 11186 15059 14126  9754     0 11725 16484  3379  6883  4250 17835  5311  8643  7729
8030  5154  6381 14906 14973 12975  3741  4627 16699 14093  3614  5330 13999 15942 11725     0  9173 13841  4874  9780 10451  9122 18555  5003
5409  8276 11223  5737 15565  7186  5549  6948  7735 16358 11306  3975  5168 11588 16484  9173     0 13413 12955 12642  3450 17931 11293 14142
17959 18175  7602  9378  3627 13212 15863 11282  7340  6543 14452 14165 11153  9288  3379 13841 13413     0  9181  4145 14775  8581  7147 10945
12822  9367  3236 16868 10307 17301  8377  6102 15982 10709  5954  9100 17968 15400  6883  4874 12955  9181     0  5593 15242  6278 15447  3419
17136 14840  3457 11899  7077 17106 11981  7523 10421  9458 11533 10994 13808 13331  4250  9780 12642  4145  5593     0 15914  8478 11291  8453
3267  7056 14535  5402 14478  3881  6865  9992  7492 13488 10801  6469  3728  8652 17835 10451  3450 14775 15242 15914     0 15624  9326 14308
12882  9698  8830 15921  6475 12658 12543 12321 14151  5008  6649 14453 15782  9157  5311  9122 17931  8581  6278  8478 15624     0 10157  4139
11223 13603 14748  5765  5155  6210 16177 16782  4374  5170 14987 15268  6168  2621  8643 18555 11293  7147 15447 11291  9326 10157     0 14265
11078  7266  6655 19675 10235 14138  8600  8550 18092  9146  3508 10326 17749 12782  7729  5003 14142 10945  3419  8453 14308  4139 14265     0