#!/usr/bin/perl

if (scalar(@ARGV) <= 0) {
    print "perl analyze.pl <cmd> <args>\n";
    exit(1);
}

my $cmd = shift(@ARGV);
if ($cmd eq "select") {
  &globalSelect(@ARGV);
}
if ($cmd eq "idlist") {
  &analyzeIDlist(@ARGV);
}

sub getHash {
  my ($file, $index, $spstr, $spindex, $rest) = @_;
  if (!defined $index) {
    $index = 0;
  }
  my $hash = {};
  my $fh;
  if ($file eq "-") {
    $fh = \*STDIN;
  }
  else {
    open($fh, "<$file") || die "Can't open $file\n";
  }
  while (my $in = <$fh>) {
    $in =~ s/[\r\n]//g;
    my @list = split("\t", $in, -1);
    my $id = $list[$index];
    if (defined $spstr) {
      $spindex = 0 if (!defined $spindex);
      my @l = split($spstr, $list[$index]);
      $id = $l[$spindex];
      $id =~ s/\s//g;
    }
    $hash->{$id} = [@list];
  }
  if ($file ne "-") {
    close($fh);
  }
  return $hash;
}

sub mse {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = 0.0;
  if ($start > $end) {
    return $result;
  }
  my $m = &mean($arrayref, $start, $end);
  my $result = 0.0;
  for (my $i = $start; $i <= $end; $i++) {
    $result += ($arrayref->[$i] - $m) * ($arrayref->[$i] - $m);
  }
  return $result;
}

sub max {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = $arrayref->[$start];
  if ($start > $end) {
    return $result;
  }
  for (my $i = $start; $i <= $end; $i++) {
    if ($result < $arrayref->[$i]) {
      $result = $arrayref->[$i];
    }
  }
  return $result;
}

sub min {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = $arrayref->[$start];
  if ($start > $end) {
    return $result;
  }
  for (my $i = $start; $i <= $end; $i++) {
    if ($result > $arrayref->[$i]) {
      $result = $arrayref->[$i];
    }
  }
  return $result;
}

sub median {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = 0.0;
  if ($start > $end) {
    return $result;
  }
  $result = $arrayref->[int(($start+$end)/2)];
  return $result;
}

sub mean {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = 0.0;
  if ($start > $end) {
    return $result;
  }
  for (my $i = $start; $i <= $end; $i++) {
    $result += $arrayref->[$i];
  }
  return $result / ($end - $start + 1);
}

sub stddev {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $sum = 0.0;
  if ($start >= $end) {
    return $sum;
  }
  my $sumsq = 0;
  for (my $i = $start; $i <= $end; $i++) {
    $sum += $arrayref->[$i];
    $sumsq += $arrayref->[$i] * $arrayref->[$i];
  }
  my $m = $sum / ($end - $start + 1);
  my $msq = $sumsq / ($end - $start + 1);
  my $factor = ($end - $start + 1) / ($end - $start);
  my $res = $factor * ($msq - $m * $m);
  if ($res < 0) { $res = 0; }
  return sqrt($res);
}

sub variance {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $sum = 0.0;
  if ($start >= $end) {
    return $sum;
  }
  my $sumsq = 0;
  for (my $i = $start; $i <= $end; $i++) {
    $sum += $arrayref->[$i];
    $sumsq += $arrayref->[$i] * $arrayref->[$i];
  }
  my $m = $sum / ($end - $start + 1);
  my $msq = $sumsq / ($end - $start + 1);
  my $factor = ($end - $start + 1) / ($end - $start);
  return ($msq - $m * $m);
}

sub sum {
  my ($arrayref, $start, $end) = @_;
  if (!defined $start) {
    $start = 0;
  }
  if (!defined $end) {
    $end = scalar(@$arrayref) - 1;
  }
  my $result = 0.0;
  for (my $i = $start; $i <= $end; $i++) {
    $result += $arrayref->[$i];
  }
  return $result;
}

sub fitStep {
  my ($data, $start, $end, $times) = @_;
  my $count = $end - $start + 1;
  if ($count <= 0) {
    return [0, 0, 0, 0, 0, 0, 0, 0];
  }
  if (!defined $times) {
    $times = [0 .. $end];
  }
  else {
    my @sk = sort { $times->[$a] <=> $times->[$b] } $start .. $end;
    $data = [map { $data->[$_] } @sk];
    $times = [map { $times->[$_] } @sk];
    $start = 0;
    $end = $count - 1;
  }
  my @sseArray = map { 0.0 } 0 .. ($count - 1);
  my $sum = &sum($data, $start, $end);
  my $mean = &mean($data, $start, $end);
  my $sstot = &mse($data, $start, $end);
  my $sum1 = 0.0;
  my $count1 = 0;
  my $m1 = 0.0;
  my $sum2 = $sum;
  my $count2 = $count;
  my $m2 = ($sum/$count);
  my $sum1sq = 0.0;
  my $sum2sq = $sstot;
  my $sse = $sum1sq + $sum2sq;
  for (my $i = 0; $i < $count; $i++) {
    my $entry = $data->[$i + $start];
    if (!defined $entry || $entry eq "") {
      $sseArray[$i] = $sse;
      next;
    }
    $count1 ++;
    $count2 --;
    if ($count2 == 0) {
      $sseArray[$i] = $sstot;
      next;
    }
    my $tmp = ($mean - ($entry + $sum1)/$count1);
    $sum1sq = $sum1sq + ($entry - $mean) * ($entry - $mean) - 
      $tmp * $tmp * $count1 + ($count1 - 1) * ($mean - $m1) * ($mean - $m1);
    $tmp = ($mean - ($sum2 - $entry)/$count2);
    $sum2sq = $sum2sq - ($entry - $mean) * ($entry - $mean) - 
      $tmp * $tmp * $count2 + ($count2 + 1) * ($mean - $m2) * ($mean - $m2);
    $sum1 += $entry;
    $sum2 -= $entry;
    $m1 = $sum1/$count1;
    $m2 = $sum2/$count2;
    $sse = $sum1sq + $sum2sq;
    $sseArray[$i] = $sse;
  }
  #print join(",", map { sprintf("%.2f", $_) } @sseArray), "\n-----------\n";
  my $bestSse;
  my $bestIndex = 0;
  my $dof = 0;
  for (my $i = 0; $i < $count ; $i++) {
    if ($i < ($count - 1)) {
       next if ($times->[$i+$start] eq $times->[$i+$start+1]);
    }
    $dof++;
    my $index = $i + $start;
    if (!defined $bestSse) {
      $bestSse = $sseArray[$i];
      $bestIndex = $index;
    }
    if ($sseArray[$i] < $bestSse) {
      $bestSse = $sseArray[$i];
      $bestIndex = $index;
    }
    #{# Debug code start
    #  $entry = $data->[$index];
    #  $sum1sq = &mse($data, $start, $index);
    #  $sum2sq = &mse($data, $index + 1, $end);
    #  $sse = $sum1sq + $sum2sq;
    #  printf "%.2f\t%.2f\t%.2f\t%.2f\n", $sseArray[$i], $entry, $sse, $sum1sq;
    #}# Debug code end
  }
  #printf "bestSse = %.2f\t bestIndex=%d\n", $bestSse, $bestIndex;
  $m1 = &mean($data, $start, $bestIndex);
  $m2 = &mean($data, $bestIndex + 1, $end);
  my $thr = ($m1 + $m2)/2.0;
  #{# Debug code start
  #  $sum1sq = &mse($data, $start, $bestIndex);
  #  $sum2sq = &mse($data, $bestIndex + 1, $end);
  #  $sse = $sum1sq + $sum2sq;
  #  if ($sse != $bestSse) {
  #    print STDERR "SSE calculation is wrong :", $sse, ":", $bestSse, "\n";
  #  }
  #}# Debug code end
  my $label = 0;
  if ($m1 < $m2) {
    $label = 1;
  }
  else {
    $label = 2;
  }
  my $statistic = 0;
  if ($bestSse > 0) {
    if ($dof > 4) {
      $statistic = ($sstot - $bestSse)/3/($bestSse/($dof - 4));
    }
    else {
      $statistic = ($sstot - $bestSse)/2/$bestSse;
    }
  }
  return [$bestIndex, $bestSse, $sstot, $statistic, $m1, $m2, $thr, $label];
}

sub globalSelect {
  my ($file, $rest) = @_;
  my $hash = &getHash($file, 2);
  my $ids = [keys %{$hash}];
  my @data = sort { $a <=> $b } map {$hash->{$_}->[0]} @$ids;
  my $res = &fitStep(\@data, 0, $#data);
  my $thr = $res->[6];
  print STDERR $thr, "\n";
  $ids = [sort { $hash->{$b}->[1] <=> $hash->{$a}->[1]} grep {$hash->{$_}->[0] > $thr} @$ids];
  foreach my $id (@$ids) {
    print join("\t", @{$hash->{$id}}), "\n";
  }
}

sub analyzeIDlist {
  my ($pclfile, $listfile, $rest) = @_;
  open(FL, "<$listfile") || die "Can't open $listfile\n";
  my $hash = {};
  while (<FL>) {
    s/[\r\n]//g;
    my ($id, $rest) = split("\t");
    $id =~ s/\s//g;
    $hash->{$id} = 1;
  }
  close(FL);
  open(FL, "<$pclfile") || die "Can't open $pclfile\n";
  while (<FL>) {
    s/[\r\n]//g;
    my $line = $_;
    my ($id, $rest) = split("\t");
    $id =~ s/\s//g;
    if (defined $hash->{$id}) {
      print "$line\n";
    }
  }
  close(FL);
}

