#!/usr/bin/perl

#Copyright (C) 2006-2008 Keio University
#(Kris Popendorf) <comp@bio.keio.ac.jp> (2006)
#
#This file is part of Murasaki.
#
#Murasaki is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#Murasaki is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with Murasaki.  If not, see <http://www.gnu.org/licenses/>.

use File::Basename;
use Getopt::Long;
use Pod::Usage;
#use Data::Dump qw{dump};

use strict;
my ($help,$man,$opt_prefix);

my ($opt_x)=("auc.tfidf");
my $samples;
my $fn;
my $format="png";
my $lwd=3;
my ($opt_log,$opt_clean,$opt_nofstats,$fstats,%avg,$drawAvg,$maxsamples,$nofn,$opt_title,$opt_names,$cleanup,$opt_noxlab,%defs,@def_order,$opt_n,$opt_bw,$opt_denopts,$opt_xlab,$opt_special);

GetOptions('help|?' => \$help, man => \$man, "x=s"=>\$opt_x,
	   'log=s'=>\$opt_log, clean=>\$opt_clean, 'format=s'=>\$format,
	   pdf=>sub{$format='pdf'},'lwd=f'=>\$lwd,
	   avg=>\$drawAvg,
	   'noavg'=>sub{$drawAvg=0},nofstats=>\$opt_nofstats,
	   'title=s'=>\$opt_title,
	   'names=s'=>\$opt_names,cleanup=>\$cleanup,noxlab=>\$opt_noxlab,
	   'def=s%'=>sub{$_[1]=~s/\W/./g;
			 print "Adding definition for $_[1]\n";
			 $defs{$_[1]}=$_[2];push(@def_order,$_[1])},
	   'n=i'>\$opt_n,'bw=f'=>\$opt_bw,'denopts=s'=>\$opt_denopts,
	   'xlab=s'=>sub{unless($_[1]){$opt_noxlab=1;return 1;}
			 $opt_xlab=$_[1]},
	   special=>\$opt_special
	  ) or pod2usage(1);
pod2usage(1) if $help or $#ARGV<0;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;

foreach my $file (@ARGV){
  my ($basename,$path,$suffix)=fileparse($file,qr/\.[^.]+/);
  $path=~s!/$!!; #kill trailing / if any
  my $basefile="$path/$basename";
  $basefile.="$suffix" if $suffix;

  print "Basefile: $basefile\n";

  my (@fields);
  open(my $fh,$basefile);
  $_=<$fh>;
  chomp;
  @fields=split(/\t/,$_);
  print "Fields available: ".join(", ",@fields)."\n";

  my @x;
  @x=split(/,/,$opt_x);
  if ($opt_x eq 'all') {
    @x=grep {$_ ne 'label'} @fields;
  }
  foreach(@x){
    s/\W/./g;
  }

  my $plotname=join("-",@x);
  print "Making histogram graph for $plotname\n";
  my $rout="$basefile.$plotname.hist.$format";
  my $rsrc="$basefile.$plotname.hist.R";
  
  my $outputter=$format ne 'pdf' ? 
    qq!bitmap(file="$rout",type="png16m",width=10,height=7,res=96)!:
      qq!pdf(file="$rout",width=10,height=7)!;

  my $densityOpts=",bw=$opt_bw" if $opt_bw;
  $densityOpts.=",n=$opt_n" if $opt_n;
  $densityOpts.=",$opt_denopts" if $opt_denopts;
  
  #do the R output
  open(my $R,">$rsrc");
  print $R <<ENDTEXT;
$outputter
dat<-read.delim('$basefile')
ENDTEXT

  if($cleanup){
    foreach my $x (@x){
    print $R <<ENDTEXT;
dat[,'$x']<-as.numeric(gsub("[^0-9.]","",dat[,'$x']));
ENDTEXT
    }
  }

  #define defs
  foreach my $term (@def_order){
    my $val=$defs{$term};
    print "$term defined as $val\n";
    foreach my $s (@fields){
      $s=~s/\W/./g;
      $val=~s/\{$s\}/dat[,'$s']/g;
    }
    print "$term refiend to $val\n";
    print $R "dat[,'$term']<-$val\n";
    push(@fields,$val);
  }

  my @legendTerms=($opt_names ? split(/,/,$opt_names):@x);

#determine ranges
  foreach my $x (@x){
    if($x eq $x[0]){
      print $R <<ENDTEXT;
x<-subset(na.omit(dat[,'$x']),T);
xd<-density(x$densityOpts);
y<-xd\$y;

xl<-c(min(xd\$x),max(xd\$x));
yl<-c(min(y),max(y));
ENDTEXT
    }else{
      print $R <<ENDTEXT;
x<-subset(na.omit(dat[,'$x']),T);
xd<-density(x$densityOpts);
y<-xd\$y;

xl<-c(min(xd\$x,xl[1]),max(xd\$x,xl[2]));
yl<-c(min(y,yl[1]),max(y,yl[2]));
ENDTEXT
    }
  }

#plot each
  my $i=2;
  foreach my $x (@x) {
    print $R <<ENDTEXT;
x<-subset(na.omit(dat[,'$x']),T);
xd<-density(x$densityOpts);
y<-xd\$y;
ENDTEXT

    my $color=$i;
    if($i==2){
      
      my $title="$basename $plotname density";
      $title=$opt_title if $opt_title;
      my $xlab_val="Value";
      unless($#x>0) { #if only 1 x...
	$color=1;
	$xlab_val=$x;
      }
      $xlab_val=$opt_xlab if $opt_xlab;
      my $xlab="xlab='$xlab_val'," unless $opt_noxlab;
      
      print $R <<ENDTEXT;
plot(xd,xlim=xl,ylim=yl,col=$color,lwd=$lwd,${xlab}ylab='Density',main='$title',log='$opt_log')
ENDTEXT
    }else{
      print $R <<ENDTEXT;
lines(xd,col=$color,lwd=$lwd)
ENDTEXT
    }
    if ($drawAvg){
      print $R <<ENDTEXT;
text(median(x)*1.01,median(y)*1.01,adj=c(0,1),labels=paste("median: ",median(x)));
ENDTEXT
    }
    my $legendTerm=$legendTerms[$i-2];

    if($opt_special){
      print $R <<ENDTEXT;
text(1.5,10-$i,col=$color,labels=paste("$legendTerm >= 1: ",100*round(length(x[x>=1])/length(x),4),"%"));
ENDTEXT
    }
  }continue{$i++}
  print $R "abline(v=1,lty=2,lwd=$lwd/2,col=1)\n" if $opt_special;


  if($#x>0){ #draw a legend
    my $legendpos="xl[1],yl[2]";
    my $legendcols=join(",",2..($#x+2));
    my $names=join(",",map {qq!"$_"!} @legendTerms);
    print $R "legend($legendpos,c($names),col=c($legendcols),lwd=$lwd);\n";
  }

  close($R);
  system("R --vanilla < $rsrc");
}


sub sum {
  my $sum=0;
  grep {$sum+=$_} @_;
  return $sum;
}

sub mean {
  my $total;
  foreach(@_){
    $total+=$_;
  }
  return $total/($#_+1);
}

sub min {
  my $best=$_[0];
  foreach(@_){
    $best=$_ if $_<$best;
  }
  return $best;
}

sub max {
  my $best=$_[0];
  foreach(@_){
    $best=$_ if $_>$best;
  }
  return $best;
}

sub pickOne {
  my ($ps1,$ps2,@opts)=@_;
  print $ps1."\n";
  print map {($_==0 ? "[$_]":" $_ ").": $opts[$_]\n"} 0..$#opts;
  my $res;
  do{
    print $ps2;
    $res=<STDIN>;
    chomp $res;
  }while($res && ($res<0 or $res>$#opts));
  return $opts[$res];
}

__END__

=head1 NAME

stat-histo.pl -- density plot of some given stat

=head1 SYNOPSIS

stat-histo.pl <input1> [input2 ...]

=head1 OPTIONS

Input file should be some stat summary (like generated by gatherstats.pl)
default stat is: auc-tfidf

 Other options:
--stat|predictor|x specifies what stat to use as a predictor
--log can apply log scale to x or y or xy axes
--fn can specify a FN count for calculating sensitivity
--avg enable plotting of averages
--noavg don't plot the line for averages
--lwd can specify line weight
--format can specify output file format (default png)
--pdf set format to pdf
--cleanup forces cleanup of values (ie: removes % and reparses as numeric)
--bw set bandwidth for density function calculation
--n samples for density function calculation
--denopts add extra opts for density function (separate with commas)
--def define custom fields (input fields from other fields can be taken
    by quoting with {} (eg: --def tough={size}/{speed})
