#!/usr/bin/perl -w

=head1 NAME

normalize_data.pl -- perform quantile normalization

=head1 SYNOPSIS

B<normalize_data.pl> [options] <raw data files OR directories>

Normalize data using quantile normalization.  
Input files should be .raw files generated by raw_file_maker.pl.

=cut

##############################################################################

=head1 REQUIRES

FileHandle,
File::Basename,
File::Spec,
Getopt::Long,
Pod::Usage,
GD (only when drawing images),

=cut

use FileHandle;
use File::Basename;
use File::Spec;
use Getopt::Long;
require Pod::Usage;
use strict;

##############################################################################

### CONSTANTS

my $PROGRAM = $0;
$PROGRAM =~ s;^.*/;;;

my $VERSION = '$Revision: 31 $';
$VERSION =~ s/^\$Revision: //;
$VERSION =~ s/\$$//;
$VERSION = 0 unless $VERSION;
$VERSION = sprintf "%3.1f", $VERSION/10;

my $FONT = find_font("/usr/X11R6/lib/X11/fonts/TTF/luxisr.ttf");
sub find_font {
    my $default = shift;
    return $default if -f $default;
    foreach my $path (qw/C:\WINDOWS\Fonts C:\WINNT\Fonts/) {
	if (-d $path) {
	    my $ft = File::Spec->catfile($path, "ARIAL.TTF");
	    return $ft if -f $ft;
	}
    }
}

my $SCALE_MEAN = 1500;

##############################################################################

=head1 DESCRIPTION

This program performs quantile normalization and mean normalization
on data from TAG3 or TAG4 chips.
Quantile normalization assumes all experiments in a group have similar
distributions, so the user should take care to include only experiments 
with similar properties in each group to be normalized.

Each tag group on the chip is normalized independently.  On the TAG4 chip, the
tags are grouped by uptag and downtag; on the TAG3 chip, the tags are grouped
by uptag:sense, uptag:antisense, downtag:sense, and downtag:antisense.

The mean normalization implemented in this script scales each tag in a tag
group by the mean of all values in the tag group.  The scaling factors have
been designed to scale each tag group to have an average value of 1500.

=head2 Input files

This script uses the raw files made by C<raw_file_maker.pl> as input.  Input 
file names are expected to have the format "expname.raw".
The data used from these files are the columns strain_tag and average; other
columns are ignored.

=head2 Distribution file

Quantile normalization calculates an ideal distribution which is then used
as the values for normalization.  The distribution for a given input set
can be printed to a file using the C<--savedistr> option.  This file is
tab-delimited with a header listing the experiments used to generate the
distribution followed by three columns:  tag_group, rank, and intensity.

The distribution file can be specified with the C<--usedistr> option to 
apply that particular distribution to the input files instead of calculating 
a new distribution.

=head2 Output files

Output files are called C<expname.norm> where C<expname> is parsed
from the name of the input file used to create the values in the output file.
A separate output file is generated for each input raw file.  
The output file is tab-delimited with columns:  strain_tag, norm_quantile,
and norm_mean.  Each row lists the normalized value for a given strain tag 
(a combination of strain and tag group).

=cut

##############################################################################

=head1 OPTIONS

=over 4

=item B<--usedistr> <file>

Use the distribution in <file>.

=item B<--savedistr> <file>

Save ideal distribution to file named <file>.

=item B<--imgdir> <dir>

Generate images (PNG format) showing both the normalized distribution
curve and raw intensity distribution curve for each cel file
and save images in the directory specified.

=item B<--od> <dir>

=item B<--outdir> <dir>

Directory to save created files.  (Default: same directory as raw data file)

=item B<--man>

Manual page.

=back

=head1 MORE OPTIONS

=over 4

=item B<--patt> <patt>

Limit input to files matching <patt>. 

=item B<--median>

Use median instead of mean to calculate ideal distribution.

=item B<--nonorm>

Do not write .norm files. 

=item B<--nomnorm>

=item B<--nomeannorm>

Do not perform mean normalization.

=item B<--noqnorm>

=item B<--noquannorm>

Do not perform quantile normalization.

=item B<--font> <fontfile>

Path to font file for creating image.

=item B<--v>

Verbose output.

=back

=cut


sub setup {
    Pod::Usage::pod2usage(-verbose => 1) unless @ARGV;

    my $man = 0;
    my $help = 0;

    Getopt::Long::Configure('pass_through');
    Getopt::Long::GetOptions('help|?' => \$help, man => \$man);
    Pod::Usage::pod2usage(-verbose => 1,
			  -exitstatus => 0) if $help;
    Pod::Usage::pod2usage(-verbose => 2,
			  -exitstatus => 0) if $man;

    my $p = {
	IMGDIR => undef,
	DISTR_FILE => undef,
	USE_DISTR => undef,
	USE_MEDIAN => 0,
	NORM => 1,
	OUTDIR => undef,
	EXT => '.raw',
	PATT => '\w',
	TAGINFO => {},
	MNORM => 1,
	QNORM => 1,
        FONT => $FONT,
        V => 0,
    };

    Getopt::Long::Configure('no_pass_through');
    Getopt::Long::GetOptions(
            'version' => sub { print "$PROGRAM, version $VERSION\n\n";
			       Pod::Usage::pod2usage(-verbose => 0,
						     -exitstatus => 0) },
	    'imgdir=s' => \$$p{IMGDIR},
	    'savedistr=s' => \$$p{DISTR_FILE},
	    'usedistr=s' => \$$p{USE_DISTR},
	    'median' => \$$p{USE_MEDIAN},
	    'norm!' => \$$p{NORM},
	    'v' => \$$p{V},
	    'mnorm|meannorm!' => \$$p{MNORM},
	    'qnorm|quannorm!' => \$$p{QNORM},
	    'od|outdir=s' => \$$p{OUTDIR},
	    'ext=s' => \$$p{EXT},
	    'patt=s' => \$$p{PATT},
        ) || die "\n";

    @ARGV || die "Need files to process.\n";
    my $infiles = get_input_files($p, @ARGV);
    unless ($infiles && %$infiles) { die "No files to process.\n"; }
    if ($$p{IMGDIR}) { require GD; }
    ($infiles, $p);
}

sub get_input_files {
    my $p = shift;
    my %infiles;

    my %suf = map {$_=>1} ($$p{EXT}, lc $$p{EXT}, uc $$p{EXT});
    my @suf = keys %suf;
    my $ext = join("|", map { quotemeta($_) } @suf);
    foreach my $fp (@_) {
        if (-d $fp) { # directory of raw files
            opendir(DIR, $fp) || die "Dir $fp: $!";
            my @files = grep(/$ext$/, readdir(DIR));
            closedir(DIR);
            foreach my $f (@files) {
                my $base = basename($f, @suf);
                $infiles{$base} = File::Spec->catfile($fp, $f);
            }
        } elsif (-f $fp && $fp =~ /$ext$/) { # raw file
            my $base = basename($fp, @suf);
            $infiles{$base} = $fp;
        }
    }
    return %infiles ? \%infiles : undef;
}

##############################################################################

my ($infiles, $p) = setup();

my $data = process_input($infiles, $p);
my $numexp = my @exps = sort keys %$data;
my ($ranks, $sorted_data) = rank_and_sort_data($data, $p);
my $distr = $$p{USE_DISTR} ?  get_precalculated_distribution($p) :
                calculate_ideal_distribution($sorted_data, \@exps, $p);
print_distribution($distr, $p) if $$p{DISTR_FILE};

exit unless $$p{NORM} || $$p{IMGDIR};
print STDERR "\nAssigning normalized values.\n";
my $num = 0;
foreach my $cel (@exps) {
    printf STDERR "  %d/%d) $cel ", ++$num, $numexp;
    my $norm = reassign_values($distr, $$data{$cel}, $$ranks{$cel}, $p);
    my $outf = outfile_name("$cel.norm", $$infiles{$cel}, $$p{OUTDIR});
    print_norm_data($outf, $norm, $p) if $$p{NORM};
    if ($$p{IMGDIR}) {
        my $pngf = outfile_name("$cel-norm.png",$$infiles{$cel},$$p{IMGDIR});
        create_distr_image($$data{$cel}, $distr, $pngf, $p);
    }
}

END { $$p{DBH}->disconnect if $$p{DBH} }

##----------------------------------------------------------------------------

sub process_input {
    my $infiles = shift;
    my $p = shift;

    my $numraw = my @exps = sort grep(/$$p{PATT}/, keys %$infiles);
    print STDERR "Got $numraw files to process.\n";
    die "No files to process.\n" unless $numraw;

    my %data;
    foreach my $exp (@exps) {
	my $file = $$infiles{$exp};
	my $hash = parse_raw_file($file, $p);
	unless ($hash && %$hash) { warn "  No data for $exp\n"; next; }
	$data{$exp} = $hash;
    }
    die "No data to process.\n" unless %data;
    (\%data);
}


sub parse_raw_file {
    my $file = shift;
    my $p = shift;

    my %data = ();
    my %datazyg = ();
    my $fh = new FileHandle $file or die "$file: $!";
    my @lines = grep(!/^#/ && /\w/ && s/[\r\n]*$//, $fh->getlines);
    my $head = shift @lines;
    my @fields = split(/ *\t */, $head);
    my %extras;
    foreach (@lines) {
        my @data = split(/ *\t */, $_);
        my %h;
        foreach my $f (@fields) {
            my $v = shift @data;
            $h{$f} = $v;
        }
	my $tag = $h{strain_tag};
	my $v = $h{raw_average};
	my $z = $h{zygosity};
	my $at = $h{above_thresh};
	next unless defined($v) && $v =~ /\d/;
        my $taggroup = $tag; $taggroup =~ s/.*:([updown]+tag:?)/$1/;
        next if $taggroup =~ /unused/ || !defined($v);
        $data{$taggroup}{RA}{$tag} = $v;
        $data{$taggroup}{ZY}{$tag} = $z;
        $data{$taggroup}{AT}{$tag} = $at;
        if ($z) {
            $taggroup .= ":$z";
            $datazyg{$taggroup}{RA}{$tag} = $v;
            $datazyg{$taggroup}{ZY}{$tag} = $z;
            $datazyg{$taggroup}{AT}{$tag} = $at;
        }
    }
    $fh->close;
    if (use_zygosity(\%datazyg)) {
        %data = %datazyg;
    }
    print STDERR "  $file " if $$p{V};
    foreach my $tg (keys %data) {
	my $num = scalar(keys %{$data{$tg}});
	print STDERR "\t$tg=$num" if $$p{V};
    }
    print STDERR " tags\n" if $$p{V};
    (\%data);
}

sub use_zygosity {
    my $hash = shift;

    my %z;
    foreach my $tg (keys %$hash) {
        my $num = $$hash{$tg}{RA} ? keys %{ $$hash{$tg}{RA} } : 0;
        my ($z) = values %{ $$hash{$tg}{ZY} };
        $z{$z}++;
#        print "$tg $num\n";
        return 0 if $num < 500;
    }
    my $numz = keys %z;
    return $numz > 1 ? 1 : 0;
}

##----------------------------------------------------------------------------

sub rank_and_sort_data {
    my $celdata = shift;
    my $p = shift;

    print STDERR "\nRanking and sorting.\n";
    my (%ranks, %sorthash);
    foreach my $cel (keys %$celdata) {
	foreach my $tg (keys %{ $$celdata{$cel} }) {
            my $numv = values %{ $$celdata{$cel}{$tg}{RA} };
	    $sorthash{$tg}[$numv-1] = undef;  ## set array size to max num tags
	}
    }
    my $count = 0;
    foreach my $cel (keys %$celdata) {
	print STDERR "  $cel " if $$p{V};
	foreach my $tg (keys %{ $$celdata{$cel} }) {
            my $ra_hash = $$celdata{$cel}{$tg}{RA};
            my @sortlist = reverse sort {$a<=>$b} values %$ra_hash;
	    $ranks{$cel}{$tg} = rank_data(\@sortlist, $sorthash{$tg});
	    print STDERR "\t$tg ", scalar(@sortlist) if $$p{V};
	}
	print STDERR "\n" if $$p{V};
	$count++;
    }
    print STDERR "  $count cel files ranked and sorted\n" if $$p{V};
    (\%ranks, \%sorthash);
}

sub rank_data {
    my $l = shift;
    my $s = shift;
    my %r;
    for(my $i=0; $i<@$l; $i++) {
        my $v = $$l[$i];
        push(@{ $$s[$i] }, $v);
        push(@{ $r{$v} }, $i);
    }
    (\%r);
}

sub get_index_values {
    my $size = shift;
    my $numvals = shift;

    my $num = ($numvals > $size/2) ? $size-$numvals : $numvals;
    my @numlist = ();
    my $val = sprintf "%5.2f", ($size/(1+$num));
    for(my $i=0; $i<$num; $i++) {
        push(@numlist, sprintf "%1.0f", ($val*($i+1)));
    }
    my @usenums;
    if ($numvals > $size/2) {
        my ($skipnum) = @numlist ? (shift @numlist) : ($size+1);
        for(my $i=0; $i<$size; $i++) {
            if ($i==$skipnum) {
                ($skipnum) = @numlist ? (shift @numlist) : ($size+1);
            } else {
                push(@usenums, $i);
            }
        }
    }
    return @usenums ? \@usenums : \@numlist;
}

##----------------------------------------------------------------------------

sub get_precalculated_distribution {
    my $p = shift;

    my $d;
    if ($$p{USE_DISTR}) {  $d = get_distr_from_file($$p{USE_DISTR}, $p); } 
    else { return }
    return ($d);
}

sub get_distr_from_file {
    my $file = shift;
    my $p = shift;

    print STDERR "\nGetting distribution from $file.\n";
    my $fh = new FileHandle $file or die "  $file: $!";
    my @lines = $fh->getlines; foreach (@lines) { s/[\r\n]+// }
    $fh->close;
    if (my ($flist) = grep(/^#Files:/, @lines)) {
        $flist =~ s/^#Files:\s*//;
	$$p{DISTR_EXPS} = [ split(/, /, $flist) ];
    }
    @lines = grep(!/^#/, @lines);
    my $head = shift @lines;
    my @f = split(/ *\t */, $head);
    my $idx = find_fields(\@f, [qw/tag_group rank intensity/]);
    my %distr;
    foreach (@lines) {
	my @d = split(/ *\t */, $_);
 	my $tg = $d[$$idx{tag_group}];
 	my $rk = $d[$$idx{rank}] - 1;
 	my $v = $d[$$idx{intensity}];
	$distr{$tg}[$rk] = $v;
    } 
    $fh->close;
    (\%distr);
}

sub calculate_ideal_distribution {
    my $sorted_data = shift;
    my $exps = shift;
    my $p = shift;
    
    $$p{DISTR_EXPS} = $exps;
    my ($csub, $ctype) = $$p{USE_MEDIAN} ? (\&median, 'medians') : 
	                                      (\&mean, 'means');
    my %distr;
    print STDERR "\nCalculating ideal distribution using $ctype.\n";
    foreach my $tg (keys %$sorted_data) {
	foreach my $list (@{ $$sorted_data{$tg} }) {
	    my $num = &$csub($list);
	    my $v = defined($num) ? sprintf "%10.4f", $num : '';
	    $v =~ s/^ *//; $v =~ s/\.?0+$// if $v =~ /\d\.\d/;
	    push(@{ $distr{$tg} }, $v);
        }
    }
    (\%distr);
}

sub print_distribution {
    my $distr = shift;
    my $p = shift;

    unless (ref($distr) eq 'HASH' && %$distr) { 
	warn "  No distribution to print.\n";
    }
    my $outfile = $$p{DISTR_FILE};
    print STDERR "  Printing distribution to $outfile\n    ";
    my $ofh = new FileHandle ">$outfile" or die ">$outfile: $!";
    if (ref($$p{DISTR_EXPS}) eq 'ARRAY') {
        print $ofh "#Files: ", join(", ", @{ $$p{DISTR_EXPS} }), "\n";
    }
    print $ofh "tag_group\trank\tintensity\n";
    foreach my $tg (sort keys %$distr) {
	my $list = $$distr{$tg};
	for(my $i=0; $i<@$list; $i++) {
	    printf $ofh "%s\t%d\t%s\n", $tg, $i+1, $$list[$i];
	}
    }
    $ofh->close;
}

##----------------------------------------------------------------------------

sub reassign_values {
    my $distr = shift;
    my $celdat = shift;
    my $ranks = shift;
    my $p = shift;

    my %norm;
    foreach my $tg (sort keys %$celdat) {
        my $scale;
        if ($$p{MNORM}) {
            my @rawv = values %{ $$celdat{$tg}{RA} };
            my $mn = mean(\@rawv);
            $scale = $mn/$SCALE_MEAN;
#            print STDERR "    $tg: $mn\tscale $scale\n";
        }
	foreach my $tag (sort keys %{ $$celdat{$tg}{RA} }) {
	    my $rawv = $$celdat{$tg}{RA}{$tag};
            next unless defined($rawv);
	    $norm{zygosity}{$tag} = $$celdat{$tg}{ZY}{$tag};
	    $norm{above_thresh}{$tag} = $$celdat{$tg}{AT}{$tag};
            if ($$p{QNORM}) {
	        my $rks = $$ranks{$tg}{$rawv};
	        my $normv = norm_val($rks, $$distr{$tg});
	        $norm{norm_quantile}{$tag} = $normv;
            }
            if ($scale) {
	        $norm{norm_mean}{$tag} = sprintf "%7.2f", $rawv/$scale;
            }
	}
    }
    (\%norm);
}

sub norm_val {
    my $rks = shift;
    my $dis = shift;

    ## Alternative: break ties by using average value
    # my @v; foreach my $r (@$rks) { push(@v, $$dis[$r]); }
    # my $normval = mean(\@v);

    my $r = shift @$rks;
    return $$dis[$r];
}

sub print_norm_data {
    my $outf = shift;
    my $n = shift;
    my $p = shift;

    my @cols;
    push(@cols, 'norm_quantile') if $$p{QNORM};
    push(@cols, 'norm_mean') if $$p{MNORM};
    unless (@cols) { print STDERR "  No data to print\n"; return; }
    push(@cols, 'above_thresh', 'zygosity');
    print STDERR "  Printing normalized data to $outf\n";
    my $ofh = new FileHandle ">$outf" or die ">$outf: $!";
    print $ofh join("\t", 'strain_tag', @cols), "\n";
    foreach my $tag (sort keys %{ $$n{$cols[0]} }) {
        my @v = ($tag);
        foreach my $f (@cols) { 
            push(@v, defined($$n{$f}{$tag}) ? $$n{$f}{$tag} : '');
        }
	printf $ofh join("\t", @v)."\n";
    }
    $ofh->close;
}


##----------------------------------------------------------------------------

sub create_distr_image {
    my $celtagdat = shift;
    my $disdat = shift;
    my $pngf = shift;

    print STDERR "  Creating image $pngf\n";
    my %celdat;
    foreach my $tg (keys %$celtagdat) {
	my @cel = reverse sort {$a<=>$b} values(%{ $$celtagdat{$tg}{RA} });
	$celdat{$tg} = \@cel;
    }
    my ($w, $h, $xs, $ys) = get_dimensions($disdat, \%celdat, 800);
    my $im = new GD::Image($w, $h);
    my $c = allocate_colors($im);
    my $font = $$p{FONT};
    generate_axes($im, $font, $c, $xs, $ys);
    foreach my $tg (keys %$disdat) {
	my $numx = scalar(@{$$disdat{$tg}});
        for(my $i=0; $i<$numx; $i++) {
	    my $ycel = yval($celdat{$tg}[$i], $$ys{$tg});
	    my $ydis = yval($$disdat{$tg}[$i], $$ys{$tg});
	    my $x = xval($numx-$i, $$xs{$tg});
            $im->rectangle($x-2, $ydis-3, $x+2, $ydis+1, $$c{blu}) 
		if defined($ydis);
            $im->ellipse($x-1, $ycel-1, 5, 5, $$c{blk}) if defined($ycel);
	    unless ($i) {
		my $l1 = int(length(int($celdat{$tg}[$i])+1)*4.8) + 9;
		my $l2 = int(length(int($$disdat{$tg}[$i])+1)*4.8) + 9;
        	$im->stringFT($$c{blu}, $font, 8, 0, $x-$l2, $ydis+2,
		    int($$disdat{$tg}[$i]));
        	$im->stringFT($$c{blk}, $font, 8, 0, $x-$l1, $ycel+2,
		    int($celdat{$tg}[$i]));
	    }
        }
    }
    my $ofh = new FileHandle ">$pngf" or die ">$pngf: $!";
    binmode($ofh);
    print $ofh $im->png;
    $ofh->close;
}

sub get_dimensions {
    my $dat1 = shift;
    my $dat2 = shift;
    my $dim = shift;
    my $dim2 = int($dim/2);
    my $margin = 50;

    my @tg = sort keys %$dat1;
    unless (@tg == 2 || @tg == 4) { 
	warn "Strange number of tag groups (@tg)\n"; return undef }
    my ($w, $h, $plotw, %xs, %ys);
    if (@tg >= 2) {
	$w = $dim;
	$h = $dim2;
	$plotw = $dim2 - 2*$margin;
	$xs{$tg[0]}{s} = $margin;
	$xs{$tg[0]}{e} = $dim2-$margin;
	$ys{$tg[0]}{s} = $margin;
	$ys{$tg[0]}{e} = $dim2-$margin;
	$ys{$tg[0]}{top} = $dim2;
	$xs{$tg[1]}{s} = $dim2 + $margin;
	$xs{$tg[1]}{e} = $dim-$margin;
	$ys{$tg[1]}{s} = $margin;
	$ys{$tg[1]}{e} = $dim2-$margin;
	$ys{$tg[1]}{top} = $dim2;
    } 
    if (@tg == 4) {
	$h = $dim;
	$plotw = $dim2 - 2*$margin;
	$xs{$tg[2]}{s} = $margin;
	$xs{$tg[2]}{e} = $dim2-$margin;
	$ys{$tg[2]}{s} = $dim2+$margin;
	$ys{$tg[2]}{e} = $dim-$margin;
	$xs{$tg[3]}{s} = $dim2 + $margin;
	$xs{$tg[3]}{e} = $dim-$margin;
	$ys{$tg[3]}{s} = $dim2+$margin;
	$ys{$tg[3]}{e} = $dim-$margin;
	$ys{$tg[0]}{top} = $dim;
	$ys{$tg[1]}{top} = $dim;
	$ys{$tg[2]}{top} = $dim;
	$ys{$tg[3]}{top} = $dim;
    } 
    foreach my $tg (@tg) { 
	$xs{$tg}{max} = scalar(@{ $$dat1{$tg} }); 
	my $m2 = $$dat2{$tg} ? max($$dat2{$tg}) : 0;
	$ys{$tg}{max} = max([max($$dat1{$tg}), $m2]);
	$xs{$tg}{inc} = $plotw/$xs{$tg}{max};
	$ys{$tg}{inc} = $plotw/$ys{$tg}{max};
    }
    ($w, $h, \%xs, \%ys);
}

sub xval {
    my $x = shift;
    my $xs = shift;
    my $xval = $$xs{s} + int($x*$$xs{inc});
    ($xval);
}

sub yval {
    my $y = shift;
    my $ys = shift;
    return unless defined($y);
    my $yval = $$ys{top} - ($$ys{s} + int($y*$$ys{inc}));
    ($yval);
}

sub generate_axes {
    my $im = shift;
    my $font = shift;
    my $c = shift;
    my $xs = shift;
    my $ys = shift;

    foreach my $tg (keys %$xs) {
	my ($sx, $ex) = ($$xs{$tg}{s}, $$xs{$tg}{e});
	my $sy = $$ys{$tg}{top} - $$ys{$tg}{s};
	my $ey = $$ys{$tg}{top} - $$ys{$tg}{e};
        $im->line($sx, $sy, $ex, $sy, $$c{blk}); # x axis
        $im->line($sx, $sy, $sx, $ey, $$c{blk}); # y axis
        # x tick marks
	for(my $i=2000; $i<$$xs{$tg}{max}-600; $i+=2000) {
	    my $xcoord = xval($i, $$xs{$tg});
            $im->line($xcoord, $sy-5, $xcoord, $sy+5, $$c{blk});
            $im->stringFT($$c{blk}, $font, 9, 0, $xcoord-15, $sy+15, $i);
	}
        $im->line($ex, $sy-5, $ex, $sy+5, $$c{blk});
        $im->stringFT($$c{blk}, $font, 9, 0, $ex-15, $sy+15, $$xs{$tg}{max});
        # y tick marks
	my $yinc = max([int($$ys{$tg}{max}/5000)*1000, 2000]);
	for(my $i=$yinc; $i<$$ys{$tg}{max}; $i+=$yinc) {
	    my $ycoord = yval($i, $$ys{$tg});
            $im->line($sx-5, $ycoord, $sx+5, $ycoord, $$c{blk});
            $im->stringFT($$c{blk}, $font, 9, 1.57, $sx-8, $ycoord+15, $i);
	}
        # x label
        $im->stringFT($$c{blk}, $font, 10, 0, $sx+100, $sy+30,
		"Data point");
        # y label
        $im->stringFT($$c{blk}, $font, 10, 1.57, $sx-30, $sy-130, 
		"Intensity");
        # legend
	my $len = max([length($tg), 10])*4+38;
	my ($lx, $ly) = ($sx+20, $ey+10);
        $im->rectangle($lx, $ly, $lx+$len, $ly+52, $$c{gry});
	my $ly_tt = $ly + 14;
        $im->stringFT($$c{blk}, $font, 9, 0, $lx+5, $ly_tt, $tg);
	my $ly_cir = $ly_tt + 15;
        $im->ellipse($lx+8, $ly_cir-1, 5, 5, $$c{blk});
        $im->stringFT($$c{gry}, $font, 8, 0, $lx+20, $ly_cir+2, "raw");
	my $ly_rct = $ly_cir + 11;
        $im->rectangle($lx+6, $ly_rct, $lx+10, $ly_rct+4, $$c{blu});
        $im->stringFT($$c{gry}, $font, 8, 0, $lx+20, $ly_rct+6, "normalized");
    }
}

sub allocate_colors {
    my $im = shift;
    my %c;
    $c{bg} = $im->colorAllocate(255,255,255);
    $c{blk} = $im->colorAllocate(0,0,0);
    $c{gry} = $im->colorAllocate(80,80,80);
    $c{blu} = $im->colorAllocate(0,0,255);
    (\%c);
}

##----------------------------------------------------------------------------

sub max {
    my $list = shift;
    return unless $list && @$list;
    return ($$list[0]) if @$list == 1;
    my ($x) = reverse sort {$a<=>$b} @$list;
    return $x;
}

sub mean {
    my $list = shift;
    return unless $list && @$list;
    return $$list[0] if @$list == 1;
    my $sum;
    foreach (@$list) { $sum += $_; }
    return $sum/scalar(@$list);
}

sub median {
    my $l = shift;
    return unless $l && @$l;
    return $$l[0] unless @$l > 1;
    my @l = sort{$a<=>$b}@$l;
    return $l[$#l/2] if @l&1;
    my $mid= @l/2;
    return ($l[$mid-1]+$l[$mid])/2;
}

##----------------------------------------------------------------------------

sub find_fields {
    my $fields = shift;
    my $find = shift;
    my %hash;
    for(my $i=0; $i<@$fields; $i++) {
	foreach my $f (@$find) {
	    $hash{$f} = $i if $$fields[$i] eq $f;
	}
    }
    (\%hash);
}

sub outfile_name {
    my $outfile = shift;
    my $file = shift;
    my $outdir = shift;

    my $path;
    if ($outdir) { $path = $outdir; } 
    else {
        # Return output file name in same directory as given file.
        my ($b, $p) = fileparse($file);
	$path = $p;
    }
    return File::Spec->catfile($path, $outfile);
}

sub print_note {
    my $fh = shift || *STDERR;
    my $count = shift;
    my $opts = shift || {};

    my $divbg = $$opts{divbg} || 10;
    my $divsm = $$opts{divsm} || 1;
    my $char = $$opts{char} || '.';

    if (!($count % $divbg)) {
        print $fh "$count";
    } elsif (!($count % $divsm)) {
        print $fh $char;
    }
}

##############################################################################

=head1 COPYRIGHT

Copyright 2007 The Board of Trustees of the Leland Stanford
Junior University. All Rights Reserved.

=cut

__END__