#!/usr/bin/env perl
use strict;
use warnings;
use Data::Dumper;
use PDL;
use PDL::NiceSlice;
use Storable qw/dclone/;

my $CountMatrix = pdl([
[ 4372,  121,  113,  162,  96,  140,  324,  387,   40,  136,  212,  219,  72,   50,  213,  596,  300,   6,   39,  428],
[  121, 2748,  100,   70,  20,  188,  170,   88,   69,   48,   85,  679,  38,   26,   46,  129,  110,  10,   30,   59],
[  113,  100, 1727,  322,  16,   92,  166,  190,   89,   37,   42,  209,  17,   18,   38,  207,  160,   4,   40,   45],
[  162,   70,  322, 3092,   8,  127,  732,  154,   68,   18,   36,  176,  11,   14,   67,  204,  120,   4,   24,   28],
[   96,   20,   16,    8, 680,    6,    6,   20,    6,   30,   30,   12,  10,   22,    6,   61,   54,   7,   16,   83],
[  140,  188,   92,  127,   6, 1367,  346,   58,   70,   35,  100,  310,  52,   18,   48,  128,  103,   6,   30,   66],
[  324,  170,  166,  732,   6,  346, 3378,  138,   63,   40,   68,  411,  29,   13,   93,  212,  192,   6,   28,  104],
[  387,   88,  190,  154,  20,   58,  138, 5788,   30,   19,   36,  136,  14,   17,   46,  230,   78,   8,   12,   40],
[   40,   69,   89,   68,   6,   70,   63,   30, 1059,   14,   39,   92,  16,   54,   34,   48,   40,  13,   96,   27],
[  136,   48,   37,   18,  30,   35,   40,   19,   14, 3217,  858,   72, 198,  144,   24,   55,  156,  12,   37, 1300],
[  212,   85,   42,   36,  30,  100,   68,   36,   39,  858, 5061,  124, 410,  331,   62,  104,  133,  20,   92,  593],
[  219,  679,  209,  176,  12,  310,  411,  136,   92,   72,  124, 3054,  48,   23,   84,  216,  172,   4,   32,   84],
[   72,   38,   17,   11,  10,   52,   29,   14,   16,  198,  410,   48, 923,   82,   12,   38,   64,   5,   24,  134],
[   50,   26,   18,   14,  22,   18,   13,   17,   54,  144,  331,   23,  82, 2191,   14,   38,   39,  59,  367,  101],
[  213,   46,   38,   67,   6,   48,   93,   46,   34,   24,   62,   84,  12,   14, 2757,  147,   66,   2,   10,   54],
[  596,  129,  207,  204,  61,  128,  212,  230,   48,   55,  104,  216,  38,   38,  147, 2326,  490,  13,   28,  100],
[  300,  110,  160,  120,  54,  103,  192,   78,   40,  156,  133,  172,  64,   39,   66,  490, 2808,   4,   31,  299],
[    6,   10,    4,    4,   7,    6,    6,    8,   13,   12,   20,    4,   5,   59,    2,   13,    4, 582,   48,   14],
[   39,   30,   40,   24,  16,   30,   28,   12,   96,   37,   92,   32,  24,  367,   10,   28,   31,  48, 1612,   42],
[  428,   59,   45,   28,  83,   66,  104,   40,   27, 1300,  593,   84, 134,  101,   54,  100,  299,  14,   42, 3784],
]);

# $r1 is a row vector of 1 - M[i,i]
# $r2 is a column vector of values from $F (actually from $Freq)
# $r1 x $r2 is the product which and it's a scalar (actually a 1x1 matrix, so that's
# why we're using ->at(0,0) to pull out that number from that 1x1 matrix).
sub RateMutation {
  my ($M,$F) = @_;
  my $Mdim = shape($M)->at(1,0);
  my $r1 = pdl(ones($Mdim) - $M->diagonal(0,1));
  my $r2 = $F->(*1);
  return pdl($r1 x $r2)->at(0,0);
}

# $Tot holds the total of all elements in the matrix $M
# $Mdim is the dimension of the matrix minus one (in our case 20 since we have a 20 x 20 square matrix)
# That's because indexing starts at 0 for piddles (the primite PDL data structures)

sub ComputeFrequencies {
  my ($M) = @_;
  my $Tot = sum($M);
  my $Mdim = shape($M)->at(1,0) - 1;
  my $F = pdl(map { $M->slice(",$_")->sum /$Tot } (0..$Mdim));
  return $F;
}

sub ComputeMutationMatrix {
  my ($C) = @_;
  my $N = zeroes(20,20);
  for(0..19) {
    $N->set($_,$_,1/$C->slice(":,$_")->sum );
  };
  my $M = $C x $N;
  return $M;
};

my $MutationMatrix = ComputeMutationMatrix($CountMatrix);
my $Freq = ComputeFrequencies($CountMatrix);

my $DBL_EPSILON;   # e is equal to the DBL_EPSILON constant
$DBL_EPSILON=1;$DBL_EPSILON/=2 while $DBL_EPSILON/2+1>1;

print "dbl_epsilon=$DBL_EPSILON\n";

my $PAM1   = $MutationMatrix->copy();
my $lnPAM1 = log($PAM1);
my $k = 1;

my $r = RateMutation($PAM1,$Freq);
while(abs($r -0.01) > $DBL_EPSILON) {
  $k *= $r * 100;
  $lnPAM1 /= $r * 100;
  $PAM1 = $lnPAM1->exp;
  $r = RateMutation($PAM1,$Freq);
  print "r=$r\n";
};

# print the 10 x 10 upper left sub-matrix of $PAM1
# not sure if this is the same as 
# the Darwin code ==>   PrintMatrix(transpose(transpose(PAM1[1..10])[1..10]),'% .5f ');
print $PAM1->slice('0:9,0:9');