#!/usr/bin/perl -w
# Construction of phylogenetic trees using UPGMA method
# ------------------------------------------------------------

sub inputfile

# Read the distance matrix ($matrixfilename)

        {

        my($matfilename)=$_[0];
        my ($i,$j);
        open MAT,$matfilename or die "ERROR: can't open the file:  $matfilename !\n";

# Get the number of the organisms

        $LINE=<MAT>;
        chomp $LINE;
        @element=split /\s+/,$LINE;
        $n=$#element+1;

   if ($n<3){
   die "ERROR: In order to build a phylogenetic tree we need at least 3 species!\n"
   }

# Distances matrix initialization

        for ($i=0;$i<=$n;$i++)
          {
          for ($j=0;$j<=$n;$j++)
             {
             $dist[$i][$j]=1000;
             }
          }
		  

# Fill distmat (dist) with the distances

        for ($i=2;$i<=$n;$i++){
        $dist[$i][1]=$element[$i-1];
        }

        for ($j=2;$j<=$n;$j++){
        $LINE=<MAT>;
        @element=split /\s/,$LINE;
          for ($i=$j+1;$i<=$n;$i++){
          $dist[$i][$j]=$element[$i-1];
          }
        }
        
		#TEST############
		print "Matrix snapshot:\n\n";
		for($z=0;$z<=$n;$z++){
			for($yz=0;$yz<=$n;$yz++){
				print "\t$dist[$z][$yz]";
			}
			print "\n";
		}
		print "\n\n";
		#END DEBUG #######
		
# OUTPUT :		
# Matrix snapshot:
#
#        1000    1000    1000    1000    1000    1000    1000    1000
#        1000    1000    1000    1000    1000    1000    1000    1000
#        1000    19.00   1000    1000    1000    1000    1000    1000
#        1000    27.00   31.00   1000    1000    1000    1000    1000
#        1000    8.00    18.00   26.00   1000    1000    1000    1000
#        1000    33.00   36.00   41.00   31.00   1000    1000    1000
#        1000    18.00   1.00    32.00   17.00   35.00   1000    1000
#        1000    13.00   13.00   29.00   14.00   28.00   12.00   1000		
        
        return @dist;
}

# ---------------------------------------------------------

sub inputmatrix

# Data read from STDIN

{
        my ($i,$j);

        print "Enter the number of the species\n";
        $n=<STDIN>;
        chop $n;

   if ($n<3){
   die "ERROR: In order to build a phylogenetic tree we need at least 3 species!!\n"
   }

# Distances matrix initialization

        for ($i=0;$i<=$n;$i++)
          {
          for ($j=0;$j<=$n;$j++)
             {
             $dist[$i][$j]=1000;
             }
          }

# Read the data

        for ($i=2;$i<=$n;$i++)
          {
          for ($j=1;$j<=$i-1;$j++)
             {
             print "Enter the distance ($i,$j)\n";
             $dist[$i][$j]=<STDIN>;
             chop $dist[$i][$j];
             }
          }


        return @dist;
}

# ---------------------------------------------------------

# newmatrix: This function read the content of the global variable @dist (a variable containing the current snapshot of the
#            distance matrix), update the GLOBAL variables $min (minimum distance), $p (row index) and $q (column index) at which
#            the minimun was found and then "neutralize" (fill with 1000) the cells in the corresponding row/column pair. 
#            The function then update a GLOBAL string used to build progressively a string describing the phylogenetic tree.
#            NB: This function DOES NOT include the computed distances in the STRING describing the tree. This informatio, at
#                this point has been already computed but it is not used in the tree building process ...
#            The core of this exercise is this function 
# ---------------------------------------------------------
sub newmatrix
{
    my ($i,$j);
    mini(@dist);

    print "Grouping $q and $p \n";
# In this script we have a matrix ($n x $n) in which we "neutralize" the columns and rows corresponding to the merged species. Another possible way would be to code a script in which the size of the matrix change at runtime (i.e. by removing the columns and the rows).   
# Matrix update

    for ($i=2;$i<=$n;$i++)
       {
       for ($j=1;$j<=$i-1;$j++)
           {
           if ($i==$p || $j==$p)
               { # condition 1:
	       ########## set the value of the current cell to 1000 ########### FILL IN ###########
						print" ($i , $j : 1)";
               }
           elsif ($i!=$q && $i!=$p && $j!=$q && $j!=$p)
               { # condition 2:
               ########## set the value of the current cell to the value stored in the cell ########### FILL IN ###########
						print" ($i , $j : 2)";
               }
           elsif ($i==$q && $j<$q)
               {  # condition 3:
                $dist[$i][$j]=(($clu[$q])*$dist[$q][$j]+($clu[$p])*$dist[$p][$j])/($clu[$q]+$clu[$p]);
						print" ($i , $j : 3)";
               }
           elsif ($j==$q && $i>$q && $i<$p)
               {  # condition 4:
               ########## VCC= (clust_q*cella_iq  + clust_p*cella_pi) / (clust_q+clust_p) ### FILL IN ##################
						print" ($i , $j : 4)";
               }
           elsif ($j==$q && $i>$p)
               {  # condition 5:
               ########## VCC= (clust_q*cella_iq  + clust_p*cella_ip) / (clust_q+clust_p) ### FILL IN #################
						print" ($i , $j : 5)";
               }
           else
               {
               print "Bug!i=$i,j=$j,p=$p,q=$q\n";
               }
           }
       }
#Clustering variables update
    $clu[$q]=$clu[$q]+$clu[$p];
    $clustr[$q]=join(",","($clustr[$q]","$clustr[$p])");
#Print the current clade
    print "\nNew Tree => $clustr[$q]  => distance=$min\n\n";
	
		#DEBUG############
		print "Matrix snapshot :\n\n";
		for($z=0;$z<=$n;$z++){
			for($yz=0;$yz<=$n;$yz++){
				print "\t$dist[$z][$yz]";
			}
			print "\n";
		}
		print "\n\n";
		#END DEBUG #######

        return @dist;
}

#---------------------------------------------------------

sub mini
{

# This function find the minimum distance by reading all the values stored in the matrix and saving the value, the row index and the column index

   $min=1000;
   for ($i=2;$i<=$n;$i++)
       {
       for ($j=1;$j<=$i-1;$j++)
           {
            if($min>$dist[$i][$j])
               {
               $min=$dist[$i][$j];
               $p=$i;
               $q=$j;
               }
           }
       }
    print "The minimum distance is $min  => ";
}

#----------------------------------------------------------

sub showmatrix

# Help function that prints the content of a matrix ... a quite boring stuff. 

{
  my ($i,$j,$l);
  print "\n";

  for($l=1;$l<=$n;$l++){
  $cond[$l]=0;
  $cond[$l]=0;
  }

  $k=$k+1;

  if ($k>1){

    $kk=$k-1;
    $pp[$kk]=$p;

    for($l=1;$l<=$kk;$l++){
      for ($j=1;$j<=$n;$j++){
        if ($j==$pp[$l]){
         $cond[$j]=1;
        }
     }
    }

    for($l=1;$l<=$kk;$l++){
      for ($i=1;$i<=$n;$i++){
        if ($i==$pp[$l]){
        $cond[$i]=1;
        }
      }
    }
  }


   for ($j=1;$j<=$n;$j++){
   if ($cond[$j]==0){
      for ($i=1;$i<=$n;$i++){
      if ($cond[$i]==0){
         if ($dist[$i][$j]!=1000){		# NB: values corresponding to the "Neutralization" placeholder (1000) are omitted
         printf "%10.2f",$dist[$i][$j];}
         else{
           print "         -";
         }
      }
      }
        print "\n";
   }
   }

        print "\n";
}

#----------------------------------------------------------

# main

# ---------------------------------------------------------
my ($i);


if (scalar @ARGV==1){
   ($matrixfilename) = $ARGV[0];
   inputfile($matrixfilename);
   }
   else{
     if (scalar @ARGV==0){
     inputmatrix();
     }
     else{
     die "Wrong number of arguments!\n"
     }
   }



#################################################

print "This program builds a phylogenetic tree based on a distances matrix provided as input using the UPGMA method.\n\n";
print "N. species : $n\n\n...";

   $k=0;

# Print the initial matrix

   print "\n--------------------------------------------------------------------------------------------\n";
   print "Initial matrix:";
   print "\n--------------------------------------------------------------------------------------------\n";
   $p=1001;
   showmatrix (@dist);

# Initialization of the clusters and of cluster_string (the variable used to store the solution)

for ($i=1;$i<=$n;$i++)
   {
   $clu[$i]=1;
   $clustr[$i]=$i;
   }

# Computation and printing of the new matrices. This is performed at each step of the tree construction process in order to allow you to check the correctness of the computed values. You can check the correctness of the process step by step by looking at the content of the class slides (Biologia_Computazionale_P5). I also included some rows printing useful values. These rows are marked with a comment line including the word DEBUG.

for ($i=1;$i<=$n-1;$i++)
   {
   newmatrix(@dist);
   print "\n------------------------------------------------------------------------------------------\n";
   print "New matrix (UPGMA step $i):";
   print "\n------------------------------------------------------------------------------------------\n";
   showmatrix (@dist);
   }
