#!/usr/bin/perl
# USGS computer program source code: qwmixbq.pl
# 
# Use of this program is described in USGS Scientific Investigations Report:
# Software for Analysis of Chemical Mixture Occurrence, Distribution, and
# Potential Toxicity by Jonathon Scott, Kenneth Skach, and Patricia Toccalino
# 
# The program is written in the Perl programming language (version 5.10.1)
# and was developed and tested on a Sun V440 with 8-GB of memory running
# SunOS 5.10.  The program has also been successfully using a Lenovo with
# 64-bit Intel G630 processors and 6-GB of memory running Microsoft Windows 7.
#
# Use of trade names is for identification purposes and does not imply
# endorsement by the U.S. Geological Survey.  Perl and applicable Perl
# packages are free software, made available under an Artistic License.
# http://perldoc.perl.org/perlartistic.html
#
# Although this program has been used by the USGS, no warranty, expressed or
# implied, is made by the USGS or the United States Government as to the
# accuracy and functioning of the program and related program material nor
# shall the fact of distribution constitute any such warranty, and no
# responsibility is assumed by the USGS in connection therewith.
#
# Jon Scott, Feb 2011 : Initial coding
#            Sep 2012 : Improved internal documentation
#            Oct 2012 : Version 1.0
#
# The software produces a report that indicates :
#  the mixture
#  how many samples had all compounds in the mixture analyzed
#  how many samples had all compounds detected
#  and the mixture detection frequency.
#  Then for each exceedance threshold:  (defaults: 0.1 and 1.0)
#   how many samples had the sum of normalized concentrations exceed the threshold
#   how frequently that happened among samples where all compounds were analyzed
#   how many samples had the exceedance threshold surpassed when all the compound ratios were below the threshold
#
=head1 Usage

qwmixbq.pl -d directory [arguments]

qwmixbq.pl computes potential mixture toxicity using an additive model.
A previous run of the qwmixtures.pl program must have completed using a
benchmark reference file.  All samples and mixtures from the earlier run
are modelled.  The output report lists each mixture, how often it was
present among samples where all components were analyzed, and for each
threshold: the number and percent of samples where the normalized
concentration exceeded the threshold, and number of samples where the
threshold was exceeded by the normalized sum and not by any of the
normalized component concentrations. 

=head1 Required arguments
 -d directoryname  the directory name for reading and writing files
 
=head1 Optional arguments 
 -T directoryname  the directory for storing temporary files
 -t threshold      a normalized concentration threshold to summarize.

The "-t" threshold may be repeated to specify multiple thresholds.
If no "-t" option is supplied, then the default threshold of 0.1 and 1
are used.

=cut 
# ==================================================================
sub ErrorExit {
    print "Errors:$/  ", join(".$/  ",@error),"$/" if $#error >= 0;
    print "Alerts:$/  ", join(".$/  ",@alert),"$/" if $#alert >= 0;
    exit if $#error >= 0;
}
# ==================================================================
# Perl packages employed by this program.  See:  http://cpan.org
# ==================================================================
use IO::Zlib;
use IO::Pipe;
use File::Basename;
use File::Spec::Functions;
use File::Temp qw/ tempfile tempdir /;
$DEBUG = 0;
# ==================================================================
# Parsing of the command line
# ==================================================================
while ( ($_ = $ARGV[0]) =~ /^-\S/ ) {
    print "Processing command argument: $_$/" if $DEBUG;
    $_ = shift;
    if ( /^-d.*/ ) { $Directory = shift; }
    if ( /^-t.*/ ) { $t = shift; push (@Thresholds, $t); }
    if ( /^-T.*/ ) { $ENV{TMPDIR} = shift; }
}
if ($#ARGV >= 0) {
    push (@error,"Unexpected arguments: " . join(" ",@ARGV));
}
# ==================================================================
# -- Checking the input directory argument
# ==================================================================
if ($Directory eq "") {
    push (@error,"Missing required argument: -d directory");
}
else {
    if (! -e $Directory) {
	push (@error,"No such directory: $Directory");
    }
    else {
	if (! -d "$Directory") {
	    push(@error,"Not a directory: $Directory");
	}
	if (! -w "$Directory") {
	    push(@error,"Directory is not writeable: $Directory");
	}
    }
}
$InFile = catfile($Directory, "xlate.rdb"); 
open (IN,"<", $InFile) or
    push (@error, "Unable to find input-file: $InFile");
# ===================================================================
# -- Set default thresholds if none supplied, otherwise verify numeric  
# ===================================================================
if ($#Thresholds == -1) { @Thresholds = (0.1,1.0); }
else {
    $t = 0;
    foreach $i (@Thresholds) {
	$j = $i * -1;
	$t ++ if ($j >= 0);
    }
    push (@error,"$t non-positive thresholds supplied in: " . join(", ",@Thresholds))
	if $t > 0;
}
# ===================================================================
# -- Checking and setting the temporary-file directory
# ===================================================================
if ($ENV{TMPDIR} ne "") {
    $TempDir = $ENV{TMPDIR};
    if (! -d $TempDir ) {
	push(@alert,"Temporary-file directory does not exist: $TempDir");
	$TempDir = $Directory;
    }
    elsif (! -w "$TempDir") {
	push(@alert,"Temporary-file directory not writeable: $TempDir");
	$TempDir = $Directory;
    }
}
else {
    $TempDir = $Directory;
}
push(@alert,"Using temporary-file directory: $TempDir");
# ===================================================================
# -- Verifying QWMIXDIR software directory environment variable
# ===================================================================
if ( ! -f catfile($ENV{QWMIXDIR},"qwmixsort.pl") ) {
     if ( -f catfile(curdir,"qwmixsort.pl") ) {
	 $ENV{QWMIXDIR} = curdir;
     }
     else {
	 push(@error,"Environment variable QWMIXDIR is not properly set");
     }
 }
# ==================================================================
# -- Print report milestone
# ==================================================================
@Time = localtime;
$DateTime = sprintf ("%4.2d%2.2d%2.2d.%2.2d%2.2d%2.2d",
             $Time[5]+1900,$Time[4]+1,$Time[3],$Time[2],$Time[1],$Time[0]);
open (REPORT,">>", catfile($Directory, "report.txt")) or
    push (@error,"Unable to open report.txt for writing");
print REPORT "$/Threshold processing began at: $DateTime$/";

# ==================================================================
# -- Verify normalized-concentration input data,
#     by looking for "tea leaves" left by the qwmixtures.pl program.
# ==================================================================
$okay = 0;
head: while (<IN>) {
      if (/^#/) {
	  $okay = 1 if /^# Normalized concentrations/;
	  next head;
      }
      $_ =~ s/[\r\n]+$//;
      @F = split("\t",$_);
      $i=0;
      foreach $col (@F) {$Column{$col} = $i; $i++;}
      last head;
} close IN;
push (@error, "Mixtures were not normalized to benchmarks in xlate.rdb")
    if (! $okay);

# ==================================================================
# -- Initialize temporary, input, and output files
# ==================================================================
($BQTMP, $BQTmpFile) = tempfile( "bqtmpXXXXXX", DIR => $TempDir);
$BQTmpFile = basename($BQTmpFile);
print REPORT "Temp directory is $TempDir$/" if $DEBUG;

tie *MIX, 'IO::Zlib', catfile($Directory, "mixtures.rdb.gz"), "rb" or
    push (@error,"Unable to find input file: mixtures.rdb.gz");
tie *BQ,  'IO::Zlib', catfile($Directory, "bq.rdb.gz"), "wb" or
    push (@error,"Unable to open bq.rdb.gz for writing");
&ErrorExit;

@row = ("mixture","analyzed","present","pct_present");
foreach $col (@Thresholds) {
    push (@row,"count>" . $col,"pct>" . $col, "sum>" . $col );
}    
print BQ join("\t",@row),"\n",
    join("\t","40s",("7n") x (3 + 3 * ($#Thresholds+1))),"\n";

# ==================================================================
# -- Count number of occurrences of each mixture when the summed
#    normalized concentrations exceeds a threshold and compute 
#    exceedence percentage among samples where components were measured
# ==================================================================

Mixture: while (<MIX>) { # next mixture
    next Mixture if /^#/;
    $n++;
    next Mixture if $n < 3; # skipping RDB header of mixture-list file
    $_ =~ s/[\r\n]+$//;
    ($mix,$count,$junk) = split("\t",$_);
    @Cmpds = split(/\+/,$mix);
    open (IN,"<", $InFile);
    $n = 0;
    $checked = 0; # denominator for percentage computation
    @Exceed  = ();   @Below = ();
    for ($i = 0; $i <= $#Thresholds; $i++) {push(@Exceed,0);push(@Below,0);}

  body: while (<IN>) { # next sample
      next body if /^#/;
      $n++;
      next body if $n < 3;
      $_ =~ s/[\r\n]+$//;
      @F = split("\t",$_);
      $analyzed = 0;
      $here = 0;
      $sumbq = 0;
      @ratio = ();
      foreach $cmpd (@Cmpds) { # were mixture compounds analyzed in this sample?
	  if ($F[$Column{$cmpd}] ne "") {
	      $analyzed ++;
	      $here ++ if ($F[$Column{$cmpd}] > 0) ; # detected
	      $sumbq +=    $F[$Column{$cmpd}];
	      push (@ratio,$F[$Column{$cmpd}]);
	  }
      }
      if ($analyzed > $#Cmpds) { # all mixture compounds were analyzed
	  $checked ++;
	  if ($here == $analyzed) { # only count exceedence when all cmpds found 
	      for ($i = 0; $i <= $#Thresholds; $i++) {
		  if ($sumbq > $Thresholds[$i]) {
		      $Exceed[$i] ++;
		      $nbelow = 0;
		      for ($j = 0; $j <= $#ratio; $j++) { # detected below threshold
			  $nbelow ++ if ($ratio[$j] <= $Thresholds[$i]) ;
		      }
		      $Below[$i] ++ if $nbelow == $here;
		  }
	      }
	  }
      }
  } close IN;
    
    $pct = 0;
    $pct = sprintf ("%5.1f",100 * $count / $checked) if $checked > 0;
    @row = ($mix,$checked,$count,$pct);
    for ($i = 0; $i <= $#Thresholds; $i++) {
	$pct = 0;
	$pct = sprintf ("%5.1f",100 * $Exceed[$i] / $checked) if $checked > 0;
	push (@row,$Exceed[$i],$pct,$Below[$i]);
    }
    print $BQTMP join("\t",@row),"\n";
}
close $BQTMP;

# ==================================================================
# -- Print report milestone
# ==================================================================
@Time = localtime;
$DateTime = sprintf ("%4.2d%2.2d%2.2d.%2.2d%2.2d%2.2d",
             $Time[5]+1900,$Time[4]+1,$Time[3],$Time[2],$Time[1],$Time[0]);
print REPORT "Threshold sort beginning at: $DateTime$/";

# ==================================================================
# -- Sort sums by frequency of occurrence of the final threshold
#    (possible performance gain by threading?)
# ==================================================================
$SortKey = 4 + 3 * ($#Thresholds+1) . "nr";
$SortProg = catfile($ENV{QWMIXDIR},"qwmixsort.pl");
$TEMP = new IO::Pipe;
$TEMP -> reader ("$SortProg $TempDir $BQTmpFile $SortKey");
while (<$TEMP>)  {print BQ;}
close $TEMP;
untie *BQ; # close

# ==================================================================
# -- Print report milestone
# ==================================================================
@Time = localtime;
$DateTime = sprintf ("%4.2d%2.2d%2.2d.%2.2d%2.2d%2.2d",
             $Time[5]+1900,$Time[4]+1,$Time[3],$Time[2],$Time[1],$Time[0]);
print REPORT "Threshold processing finished at: $DateTime$/";
