#!/usr/local/bin/perl
# Paul Pavlidis. For no replicates only, two factors.
$usage = "anova-twoway-norep [-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";

# = category title
# % category name
# list of columns that apply
my ($numcat1, $numcat2, $numdup);
$/= "=";
<IN>;
$firstcat = 1;
while (<IN>) {
  chomp;
  # first category.
  ($category, @dat) = split /\n/, $_;
  $category =~ s/=//;
  print STDERR "$category";
  $incat=0;
  $numcat1++;
  
  foreach $m (@dat) {
    if($m=~/\%(.+)/) {
      $incat=1;
      $catname = $1;
      print STDERR "\n$catname ";
      if (!$firstcat) {
	$numcat2++;
      }
      next;
    }
    if ($incat) {
      push @{$cat{$category}->{$catname}}, $m;
      print STDERR "$m\t";
      if ($firstcat) {
	$n++;
      }
    }
  }
  $firstcat = 0;
  print STDERR "\n";
}
$minimum = $numcat1 * $numcat2;
$numrep = $n / $minimum;
die "You should use another program since you have replication\n" if $numrep > 1;
print STDERR "N=$n CAT1=$numcat1 CAT2=$numcat2 NUMREP=$numrep\n";
close IN;

# construct the table of replicates.
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
print STDERR "Check:\n";
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;
$remainderdf = $ffdf * $sfdf;

print STDERR "DF: Total $totaldf  Remainder $remainderdf First $ffdf Second $sfdf\n";
print STDERR "Calculate p-values as F, df1=factordf  df2=errordf\n";

if ($verbose) {
  print "label\tffdf\tsfdf\tremainderdf\tfff\tsff\tpff\tpsf\n";
} else {
  print "label\tfff\tsff\tpff\tpsf\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;

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

  # remainder
  $remainderssq = $totalssq - $ffssq - $sfssq;
  $remainderms = $remainderssq / $remainderdf;

  # F(1)
  $fff = $ffmsq/$remainderms;
  if ($fff < 0) {
    print STDERR "Warning: negative value found for f=$fff (probably roundoff error, using 0.0)\n";
    $f = 0.0
  }
  $pff = fprob($fff, $ffdf, $errordf);
  # F(2)
  $sff = $sfmsq/$remainderms;
  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);

#  print "$label\t$totalssq\t$totaldf\t$ffssq\t$ffdf\t$sfssq\t$sfdf\t$remainderssq\t$remainderdf\t$fff\t$sff\n";
  if ($verbose) {
    print "$label\t$ffdf\t$sfdf\t$remainderdf\t$fff\t$sff\t$pff\t$psf\n";
  } else {
    print "$label\t$fff\t$sff\t$pff\t$psf\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;
  $sum*=$sum;
  $sum/=$n;
  return $sum;
}

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;
  @cats = keys %$cat;
  $category = $cats[$catnum-1];
  my $grandtotal = 0;
  foreach $subcat (keys %{$cat{$category}}) {
    $sum = 0;
    foreach $trial ( @{$cat{$category}->{$subcat}}) {
      $sum+=$$data[$trial];
    }
    $grandtotal+=$sum*$sum;
  }
  $grandtotal/= $nothercat;
  $grandtotal-=$c;

  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;
}

##--