#!/usr/local/bin/perl
# Paul Pavlidis. For equal numbers of replicates only, two factors.
use Stats;
$usage = "anova-twoway-withrep  [-l: log transform; -v: verbose output; -r: format line needs to be removed] <datafile> <layoutfile>\n";
die $usage unless @ARGV > 1;
while ($ARGV[0] =~ /^-/) {
  $opt = shift @ARGV;
  if ($opt eq "-r") {
    $rdb++;
  } elsif ($opt eq "-l") {
    $log++;
  } elsif ($opt eq "-v") {
    $verbose++;
  } else {
    die "Illegal option\n";
  }
}


($data, $layout) = @ARGV;

open (IN, "<$layout") or die "Couldn't open layout $layout\n";

# format of the layout file:
# = category title
# % category name
# list of columns that apply
# etc.
my ($numcat1, $numcat2, $numdup);
$/= "=";
<IN>;

$firstcat = 1;
while (<IN>) {
  chomp;
  s/\cM//g;
  ($category, @dat) = split /\n/, $_;
  $category =~ s/=//;
  print STDERR "$category = ";
  push @maincats, $category;
  $incat=0;
  $numcat1++; # number of items in category 1.
  $k = scalar @dat;
  foreach $m (@dat) {
    if($m=~/\%(.+)/) { # start a new category.
      $incat=1;
      $catname = $1;
      print STDERR "\n$catname:\t";
      if (!$firstcat) {
	$numcat2++; # number of things in the second category.
      }
      next;
    }
    if ($incat) { # subcategories.
      print STDERR "$m\t";
      push @{$cat{$category}->{$catname}}, $m;
      if ($firstcat) {
	$n++;
      }
    }
  }
  $firstcat = 0;
  print STDERR "\n";
}
$minimum = $numcat1 * $numcat2;
$numrep = $n / $minimum;
die "You should use another program since you have no replication\n" if $numrep < 2;
print STDERR "N=$n CAT1($maincats[0])=$numcat1 CAT2($maincats[1])=$numcat2 NUMREP=$numrep\n";
#print STDERR "Calculate p-values as F, df1=factordf  df2=errordf\n";
close IN;

# construct the table of replicates. This is a bit tricky.
foreach $m (keys %cat) {
  foreach $k (keys %{$cat{$m}}) {
    foreach $q (@{$cat{$m}->{$k}}) {
      push @{$numwithcats{$q}}, $k; # associate categories with this trial.
    }
  }
}
foreach $k (sort keys %numwithcats) {
  $cat = join " ", @{$numwithcats{$k}};
  push @{$replicate{$cat}}, $k; # reverse the hash...
}

# check
foreach $k (sort keys %replicate) {
  foreach $m (@{$replicate{$k}}) {
#    print STDERR "$m $k\n";
  }
} 

# whsehw.
open (IN, "<$data") or die "couldn't open data\n";

$/="\n";
<IN>;
if ($rdb) {
  <IN>; # if rdb
  print STDERR "Removed format line\n";
} else {
  print STDERR "Assuming this is NOT rdb format!\n";
}

$totaldf = $n-1;
$celldf = $numcat1 * $numcat2 - 1;
$ffdf = $numcat1 - 1;
$sfdf = $numcat2 - 1;
$ffxsfdf = $ffdf * $sfdf;
$errordf = $numcat1*$numcat2*($numrep-1);
print STDERR "DF: Total $totaldf  Cell $celldf First $ffdf Second $sfdf FFXSF $ffxsfdf ERROR $errordf\n";
#print "label\ttotalssq\ttotaldf\tcellssq\tcelldf\tffssq\tffdf\tsfssq\tsfdf\tffxsfssq\tffxsfdf\terrorssq\terrordf\tfcell\tfff\tsff\tfxsf\n";

if ($verbose) {
  print "label\tffdf\tsfdf\tffxsfdf\terrordf\tfff\tsff\tfxsf\tpff\tpsf\tpsxff\n";
} else {
  print "label\tfff\tsff\tfxsf\tpff\tpsf\tpsxff\n";
}

while (<IN>) {
  chomp;
  s/\cM//;
  if (!$_) {
    print STDERR "Skipping blank line\n";
    next;
  }

  ($label, @data) = split /\t/, $_;

  if (scalar @data != $n) {
    print STDERR "$label: Skipping because it is missing data\n";
    next;
  }

  if ($log) {
    for ($i=0; $i<scalar @data; $i++) {
      if ($data[$i] <= 0) {
	die "You can't use the '-l' option when there are values less than or equal to zero in the data. (at $label)\n";
      }
      $data[$i] = log($data[$i]);
    }
  }

  # do anova on it.
  #calculate C;
  $C = calcC(\@data);

  # total SS
  $totalssq = calcTotalSSQ(\@data, $C);
  $totalmsq = $totalssq/$totaldf;

  #cells SS
  $cellssq = calcCellSSQ(\@data, $C, \%replicate);
  $cellmsq = $cellssq/$celldf;

  #first factor S
  $denom =  $numrep*$numcat2;
  $ffssq = calcFactSSQ(\@data, $C, \%cat, 1, $denom);
  $ffmsq = $ffssq/$ffdf;

  #second factor SS
  $denom =  $numrep*$numcat1;
  $sfssq = calcFactSSQ(\@data, $C, \%cat, 0, $denom);  
  $sfmsq = $sfssq/$sfdf;

  # 1x2 SS
  $ffxsfssq = $cellssq - $ffssq - $sfssq;
  $ffxsfmsq = $ffxsfssq/$ffxsfdf;

  # error
  $errorssq = $totalssq - $cellssq;
  $errormsq = $errorssq/$errordf;

  # F(cells)
  $fcell = $cellmsq/$errormsq;

  # F(1)
  $fff = $ffmsq/$errormsq;
  if ($fff < 0) {
    print STDERR "Warning: negative value found for f=$fff (probably roundoff error, using 0.0)\n";
    $fff = 0.0
  }
  $pff = fprob($fff, $ffdf, $errordf);

  # F(2)
  $sff = $sfmsq/$errormsq;
  if ($sff < 0.0) {
    print STDERR "Warning: negative value found for f=$sff (probably roundoff error, using 0.0\n";
    $sff = 0.0;
  }
  $psf = fprob($sff, $sfdf, $errordf);

  # F(1x2)
  $fxsf = $ffxsfmsq/$errormsq;
  if ($fxsf < 0.0) {
    print STDERR "Warning: negative value found for f=$fxsf (probably roundoff error, using 0.0\n";
    $fxsf = 0.0;
  }
  $pfxs = fprob($fxsf, $ffxsfdf, $errordf);

#  print "$label\t$totalssq\t$totaldf\t$cellssq\t$celldf\t$ffssq\t$ffdf\t$sfssq\t$sfdf\t$ffxsfssq\t$ffxsfdf\t$errorssq\t$errordf\t$fcell\t$fff\t$sff\t$fxsf\n";
  if ($verbose) {
    print "$label\t$ffdf\t$sfdf\t$ffxsfdf\t$errordf\t$fff\t$sff\t$fxsf\t$pff\t$psf\t$pfxs\n";
  } else {
    print "$label\t$fff\t$sff\t$fxsf\t$pff\t$psf\t$pfxs\n";
  }
}


##################################
# Subs
##################################
sub calcC {
  my ($data) = @_;
  my $n=0;
  my $sum = 0;
  foreach $m (@$data) {
    $sum+=$m;
    $n++;
  }
  die "Hey!!!! No zero divides, please.\n" if $n == 0;
  return ($sum * $sum) / $n;
}

sub calcTotalSSQ {
  my ($data, $c) = @_;
  my $sum = 0;
  foreach $m (@$data) {
    $sum+=$m*$m;
  }
  $sum-=$c;
  return $sum;
}

sub calcFactSSQ {
  my ($data, $c, $cat, $catnum, $nothercat) = @_;
  my $sum;
  my @cats = keys %$cat;
  my $category = $cats[${catnum}-1];
  my $trial;
  my $grandtotal = 0;
  foreach $subcat (keys %{$cat{$category}}) {
    $sum = 0;
    foreach $trial ( @{$cat{$category}->{$subcat}}) {
      $sum+=$$data[$trial];
    }
    $grandtotal+=$sum*$sum;
  }
#  print STDERR "G $grandtotal  $sum\n";
  $grandtotal/= $nothercat;
  $grandtotal-=$c;
#  print STDERR "G $grandtotal  $sum\n";

  return $grandtotal;
}

sub calcCellSSQ {
  my ($data, $c, $cat) = @_;
  my $sum = 0;
  my $grandtotal = 0;
  foreach $m (sort keys %$cat) { # for each cell
    $sum=0;
    foreach $k (@{$$cat{$m}}) { # add up all the replicates
      $sum+=$$data[$k];
    }
    $sum*=$sum;
    $grandtotal+=$sum;
  }
  $grandtotal/=$numrep;
  $grandtotal-=$c;
  return $grandtotal;
}

##--