#!/usr/bin/perl
# USGS computer program source code: qwmixtures.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.
#
# Program to determine mixtures from analytical results.
#
# Jon Scott, Oct 2008 : Initial coding
#            Sep 2012 : Improved internal documentation
#            Oct 2012 : Version 1.0
#
=head1 Usage

qwmixtures -f input_file -d output_directory [arguments]

qwmixtures determines mixtures of constituents in a group of samples.
Each row of the input file represents a sample.  The constituents to
be analyzed must be in a contiguous block of adjacent columns.  Other
columns may be present.  Column names for constituents should be
short, because mixture names will be formed by joining the detected-
constituent names with a plus sign (+).  Constituent columns are
expected to contain remark codes that qualify the measurement outcome.
No column in the input file may be named "mixture".

Output is written to the specified directory.  Successful runs produce
two output files.  A file named "raw.rdb" contains a row for every
mixture found in each sample, and reproduces all the non-constituent
columns.  This file can be useful for determining which samples and
sites produced particular mixtures.  Another file named
"mixtures.rdb.gz" contains a row for each mixture that was found at
least the minimum number of times (see "-m" argument).  The file is
compressed and sorted to show the most frequently occurring mixtures
first.

The program may run for a very long time and consume large amounts of
disk space.  Slow performance is caused by a large number of detected
constituents in each sample, because the number of mixtures grows
rapidly.

If the censoring option is selected, then unquantified detections
(remarks M and N) will be censored to non-detections, as will any
quantified detections below the specified censoring level.  The
censoring option will not censor a result that exceeded the method
upper range (remark ">").  Without the censoring option, unquantified
detections are treated as detections.

=head1 Required arguments:

 -f filename       the file name of the input data
 -d directoryname  the directory name for storing output

=head1 Optional arguments (defaults in parentheses):

 -s column  start-constituent column number or name (2)
 -e column  end-constituent column number or name (last column)
 -m n       minimum number of mixture occurrences (2)
 -c value   global censoring level (no censoring)
 -T dir     directory to use for creating temporary files (output directory)
 -v         switch to treat V-remarked data as detections (missing)

 -r file    filename of the compound-reference file    
 -p column  column name of the parameter-code column in reference file
 -b column  column name of the benchmark-concentration column in reference file
 -n column  column name of the compound-name column in reference file

=cut
# ==================================================================    
# -- Enumerate mixtures in a section of a list and sort
# ==================================================================    
sub combinate {
    local ($lower, $upper, $carry, $TempDir, @array) = @_;

    my $sortext = Sort::External->new (
				       working_dir => $TempDir
				       );
    # mem_threshold # default 1024 ** 2 * 16 (8 MiB)
    # cache_size    # default disabled

    my $num = 0;
    for (my $setsize=$lower; $setsize<=$upper; $setsize++) {
	my $engine = Math::Combinatorics->new (count => $setsize, data=>[@array]);
	
	while (my @mix = $engine->next_combination) {
	    my $mixture = join("+",sort(@mix));
	    $sortext->feed( $mixture );
            $num ++;
	}
    }
    my ($OUT, $OutFile) = tempfile( "sampleXXXXXX", 
                                    DIR => $TempDir,
                                    UNLINK => 0,
                                    SUFFIX => ".gz") or
       die "Unable to open temp file during enumeration in $TempDir\n";

    tie *OUT, 'IO::Zlib', $OutFile, "wb" or
       die "Unable to open compressed file during enumeration: $OutFile\n";    
    $sortext->finish();
    if ($carry ne "") {
        while ( defined( $_ = $sortext->fetch ) ) {
            print OUT join("\t",$_,$carry),"\n";
        }
    } else {    
        while ( defined( $_ = $sortext->fetch ) ) {print OUT $_,"\n";}
    }   
    untie *OUT;
    close OUT;
    return ($OutFile,$num);
}
# ==================================================================    
# -- Print one row in the final mixture enumeration table
# ==================================================================    
sub reducemix {
    # Prints one row to the mixture table:  mixture name and count of occurrences.
    # The row also includes the expected and actual mixture-detection frequency
    # Mixtures that occur less than the user-specified minimum number of times
    # (mincnt) are omitted.
    if ($count >= $mincnt) {
	@Compounds = split(/\+/,$previous); 
	$exp = 1;
	@PatternArray = @TemplateArray;
	foreach $compound (@Compounds) {
	    $exp *= $Freq{$compound};
	    $PatternArray[$Column{$compound}] = "[01]";
	}
	$Pattern = join("",@PatternArray);
	@measurements = grep (/$Pattern/,@SampleHitArray);
	$actual = sprintf("%12.6f",(100. * $count / ($#measurements + 1)));
	$exp  = sprintf ("%12.6f",(100. * $exp));
	$exp  =~ tr/ //d;
	print $COUNT join("\t",$previous,$count,$exp,$actual),"\n";
    }
}
# ==================================================================
# -- Process each constituent result, based on remarks and arguments.
# ==================================================================    
sub ParseResult {
    local ($Compound,$vremark,$censor,$Result) = @_;
    $Result =~ tr/ //d;
    if ($Result =~ /^([^\d.-]{0,1})(.*)/o) {
	$rmk = $1;
	$rmk =~ tr/[a-z]/[A-Z]/;
	$val  = $2;
	$val *= -1 if $rmk eq "-";
    } else {print STDERR "Parse error for $Result\n";}

    if    ($rmk eq "<") {
	$value = 0;
    }
    elsif ($rmk =~ /[MN]/o) {
	if ($censor) { $value = 0; }
	else         { $value = 0.001; } # quantify non-quantified result
    }
    elsif ($rmk eq "V") {
	if ($vremark) {
	    if ("$censor" && $val < $censor) {$value = 0;}
	    else {$value = $val;}
	}
	else {$value = "";}
    }
    elsif ($rmk eq "U") {$value = 0;}
    elsif ($rmk eq ">") {$value = $val;}
    else {
	if ($val eq "") {$value = "";}
	elsif ("$censor" && $val < $censor) {$value = 0;}
	else {$value = $val;}
    }
    if ($BenchColNo ne "") { $value = $value / 	$CompoundBenchmark{$Compound};}

    return $value;
}
# ==================================================================    
# -- Print warning and error messages, and stop if error.
# ==================================================================    
sub ErrorExit {
    print "Errors:$/  ", join(".$/  ",@error),".$/" if $#error >= 0;
    print "Alerts:$/  ", join(".$/  ",@alert),".$/" if $#alert >= 0;
    exit if $#error >= 0;
}
# ==================================================================    
# -- Calculate a factorial
# ==================================================================    
sub factorial {
    local ($value) = @_;
    my $factored = 1;
    for ($i = $value; $i = 2; $i--) {$factored *= $i;}
    return $factored;
}
# ==================================================================    
# -- Calculate the number of combinations of N items
# ==================================================================    
sub combinations {
    local ($items) = @_;
    my $sum   = Math::BigInt -> bzero();
    my $fact1 = Math::BigInt->new(&factorial($items));
    for ($i = 2; $i <= $items; $i++) {
	my $fact2 = Math::BigInt->new(&factorial($i));
	my $fact3 = Math::BigInt->new(&factorial($items - $i));
	$fact2 -> bmul ($fact3);
	my $factg = $fact1 -> copy() -> bdiv ($fact2);
	$sum -> badd ($factg);
    }
    return $sum;
}
# ==================================================================
# -- Perl packages used, see: http://cpan.org
# ==================================================================
use File::Spec::Functions ;
use File::MergeSort;
use File::Sort qw(sort_file);
use File::Temp qw/ tempfile tempdir /;
use File::Basename;
use Sort::External;
use Math::Combinatorics;
use Math::BigInt;
use IO::Pipe;
use IO::Zlib;
#use threads; # perhaps add in the future?

# ==================================================================
# -- Program defaults
# ==================================================================
$DEBUG = 0;
$Clock = 1;
$MixDef="40s";
$limit = 24; # default for QWMIXLIMIT

# ==================================================================
# -- Process command-line arguments
# ==================================================================
$CommandLine = join(" ",@ARGV);
while ( ($_ = $ARGV[0]) =~ /^-\S/ ) {
    print "Processing command argument: $_\n" if $DEBUG;
    $_ = shift;
    if ( /^-d.*/ ) { $Directory = shift; }
    if ( /^-f.*/ ) { $InFile = shift; }
    if ( /^-s.*/ ) { $startcol = shift; }
    if ( /^-e.*/ ) { $endcol = shift; }
    if ( /^-m.*/ ) { $mincnt = shift; }
    if ( /^-c.*/ ) { $censor = shift; }
    if ( /^-T.*/ ) { $ENV{TMPDIR} = shift; }
    if ( /^-v.*/ ) { $vremark = 1; }

    if ( /^-r.*/ ) { $RefFile  = shift; } # Compound reference file name 
    if ( /^-b.*/ ) { $BenchCol = shift; } # Benchmark column in reference file
    if ( /^-n.*/ ) { $NameCol  = shift; } # Compound name column in ref file
    if ( /^-p.*/ ) { $ParmCol  = shift; } # Pcode column in reference file
}
if ($#ARGV >= 0) {
    push (@error,"Unexpected arguments: " . join(" ",@ARGV));
}
# ==================================================================
# -- Check the input-filename argument
# ==================================================================
if ($InFile eq "") {
    push (@error,"Missing required argument: -f filename");
}
else {
    if (! -f $InFile) {
	push (@error,"Input file not found: $InFile");
    }
    if (! -r $InFile) {
	push (@error,"Input file not readable: $InFile");
    }
}
# ==================================================================
# -- Check the output-directory name argument
# ==================================================================
if ($Directory eq "") {push (@error,"Missing required argument: -d directory")}
else {
    if (! -e $Directory) {
	$i = mkdir $Directory;
	if ($i != 1){
	    push (@error,"Unable to create directory: $Directory");
	}
	else {
	    push (@alert,"Creating non-existing directory: $Directory");
	}
    }
    else {
	if (! -d "$Directory") {
	    push(@error,"Not a directory: $Directory");
	}
	if (! -w "$Directory") {
	    push(@error,"Directory is not writeable: $Directory");
	}
    }
}
# ==================================================================
# -- Force minimum mixture occurrence count to 2 if not specified
# ==================================================================
if ($mincnt eq "") {
    $mincnt = 2 ;
    push (@alert,"Setting minimum mixture-occurrence count to 2");
}
# ===================================================================
# -- Alert user that M/N remarked data treated as non-detections when
#     censoring is in effect
# ===================================================================
if ($censor) {
    push (@alert,"M & N remarked results will be censored to non-detections");
}
# ===================================================================
# -- Check and set 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 "$Directory") {
	push(@alert,"Temporary-file directory not writeable: $TempDir");
	$TempDir = $Directory;
    }
}
else {
    $TempDir = $Directory;
}
push(@alert,"Using temporary-file directory: $TempDir");
# ===================================================================
# -- Verify 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");
     }
 }
# ===================================================================
# -- Verify QWMIXLIMIT environment variable
# ===================================================================
if ( $ENV{QWMIXLIMIT} ne "" && $ENV{QWMIXLIMIT} !~ /[^0-9]/o ) {
    $limit = $ENV{QWMIXLIMIT} ;
}
else {
    push(@alert,"Environment variable QWMIXLIMIT not properly set, defaults to $limit");
}
	
# ===================================================================
# -- Check and read the reference-file 
# ===================================================================
if ( $RefFile ne "" ) {
    if (open (REF,"<$RefFile")) {
	$row = 0;
	RefRow: while (<REF>) {
	    next if /^#/;
	    $row ++;
	    chomp;
	    @F = split("\t",$_);
	    if ($row == 1) {
		for ($col = 0; $col <= $#F; $col++) {
		    if   ($F[$col] eq $BenchCol){$BenchColNo = $col;}
		    elsif($F[$col] eq $NameCol) {$NameColNo = $col;}
		    elsif($F[$col] eq $ParmCol) {$ParmColNo = $col;}
		}
	    }
	    elsif ($row == 2) {
		$i = $#error;
		if ($ParmCol eq "") {push(@error,"-p Parameter column required with -r reference file");}
		elsif ($ParmColNo eq ""){push(@error,"-p $ParmCol column not found in reference file -r $RefFile");}
		else {
		    if ($NameCol eq "" && $BenchCol eq "") {
			# no name nor benchmark column specified for reference file
			push (@error,"No -n name_column nor -b benchmark_column for -r $RefFile specified");
		    }
		    elsif ($NameCol ne "" && $NameColNo eq "") {
			# didn't find the name column in reference file
			push (@error,"In reference file -r $RefFile the -n $NameCol does not exist");
		    }
		    elsif ($BenchCol ne "" && $BenchColNo eq "") {
			# didn't find the benchmark column in reference file
			push (@error,"In reference file -r $RefFile the -b $BenchCol does not exist");
		    }
		}
		if ($#error != $i) { # Fatal flaws with ref file
		    push (@error,"Fatal flaws with reference file prevent further evaluation");
		    last RefRow ;
		}
	    }
	    else {
		if ($NameColNo ne "") {
		    $CompoundName{$F[$ParmColNo]} = $F[$NameColNo];
		    push (@error, "Name for $F[$ParmColNo] is missing")
			if ($F[$NameColNo] eq "");
		}
		if ($BenchColNo ne "") {
		    $CompoundBenchmark{$F[$ParmColNo]} = $F[$BenchColNo];
		    push (@error, "Benchmark for $F[$ParmColNo] is zero")
			if ($F[$BenchColNo] == 0);
		}
	    }
	}
    close REF;
    }
    else {
	push (@error, "Unable to open reference file: " . $RefFile);
    }
}
else {
    if ($NameCol ne "" || $BenchCol ne "" || $ParmCol ne "") {
	push (@alert,"No reference file supplied, arguments: -b, -n, -p are inapplicable");
    }
}
# ==================================================================
# -- Get timestamp, initialize output files
# ==================================================================
@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]);

$liststarttime = times;

open (DATA,"<$InFile");
if ($Clock) {
    open (TIMING,">", catfile(curdir(), $Directory, "timing.rdb")) ;
    print TIMING join("\t","hits","seconds"),"\n",
    join("\t","4n","9n"),"\n";
}

open (REPORT,">", catfile(curdir(), $Directory, "report.txt")) or
    push (@error,"Unable to open report.txt for writing");
print REPORT "Processing began at: $DateTime using command line:$/",
    $CommandLine,$/;

open (XLATE,">", catfile(curdir(), $Directory, "xlate.rdb")) or
    push (@error,"Unable to open xlate.rdb for writing");
print XLATE "# Normalized concentrations\n" if keys(%CompoundBenchmark);
print XLATE "# Censored concentrations to $censor\n" if $censor;

tie *RAW, 'IO::Zlib', catfile(curdir(), $Directory, "raw.rdb.gz"), "wb" or    
    push (@error,"Unable to open raw.rdb.gz for writing");
# ==================================================================
# -- Iterate through the analytical data
# ==================================================================
$row = 0;
input: while (<DATA>) {
    next input if /^#/o;
    $row ++;
    $_ =~ s/[\r\n]+$//;
    if ($row == 1) {
# ==================================================================
# -- Process the column names
# ==================================================================
	@H = split("\t",$_);
	if ($startcol eq "") {$startcol = 1;}
	elsif ($startcol =~ s/^([\d]+)$/$1/) {
#	elsif ($startcol =~ /^[\d]+$/) {
	    $startcol -= 1;
	    print REPORT "Starting column number: $startcol$/" if $DEBUG; 
	}
	else {
	    $fnd = -1;
	    for ($i=0;$i<=$#H;$i++) {$fnd = $i if $H[$i] eq $startcol;}
	    push (@error,"Start column not found: $startcol") if $fnd == -1;
	    $startcol = $fnd;
	    print REPORT "Starting column number: $startcol$/" if $DEBUG; 
	}
	if ($endcol eq "") {$endcol = $#H;}
	elsif ($endcol =~ s/^([\d+])$/$1/) {
	    $endcol -= 1;
	}
	else {
	    $fnd = -1;
	    for ($i=0;$i<=$#H;$i++) {$fnd = $i if $H[$i] eq $endcol;}
	    push (@error,"End column not found: $endcol") if $fnd == -1;
	    $endcol = $fnd;
	}
	foreach $i ((0 .. $startcol-1), ($endcol+1 .. $#H)) {
	    push (@carrycol,$i);
	    push (@carryname,$H[$i]);
	}
	$NCompounds = $endcol - $startcol + 1;
	$pot_mixtures = Math::BigInt->new(&combinations ($NCompounds));
	print "$pot_mixtures possible mixtures of $NCompounds constituents.\n";
	# =========================================================
	# When supplying benchmarks, ensure every constituent has one
	#
	if ($RefFile && $BenchColNo ne "") {
	    for ($i=$startcol; $i<=$endcol; $i++ ) {
		if (! exists $CompoundBenchmark{$H[$i]}){
		    push (@NoBenchmark,$H[$i]);
		}
	    }
	    push (@error,
		  "Reference file: $RefFile does not contain benchmarks for compounds: " .
		  join(", ",@NoBenchmark)) if $#NoBenchmark >= 0;
	}
	# =========================================================
	# Print compound names in translated input file 
	    @F = ();
	    for ($i=0; $i <= $#H; $i++ ) {
		if (! exists ($CompoundName{$H[$i]})) {
		    push (@F,$H[$i]);
		} else {push (@F,$CompoundName{$H[$i]});}
	    }
	print XLATE join("\t",@F),"\n";
	# =========================================================
	# -- Save column names and a template search pattern
	#    for later calculation of mixture-detection frequencies
	for ($i=0; $i < $NCompounds; $i++) {
	    push (@TemplateArray,".");
	    $name = $H[$startcol+$i];
	    $name = $CompoundName{$name} if exists $CompoundName{$name};
	    $Column{$name} = $i;
	}
	# =========================================================
	&ErrorExit if ($#error >= 0 || $#alert >= 0) ; # Last Escape
	next input;
    }
    elsif ($row == 2) {
	print XLATE "$_\n"; # Save column definitions for XLATE
	@F = split("\t",$_);
        @Carry = ();
        @Merges = ();
	foreach $i (@carrycol) {push (@Carry,$F[$i]);}
	print RAW join("\t","mixture",@carryname),"\n",
	join("\t",$MixDef,@Carry),"\n";
    
	next input;
    }
    else {
# ==================================================================
# -- Process each sample
# ==================================================================
	$starttime = times;
	@F = split("\t",$_);
	@Compounds = ();
	@Carry = ();
	$BinList = "";
	foreach $i (@carrycol) {push (@Carry,$F[$i]);}
	for ($i = $startcol; $i <= $endcol; $i++) {
	    $result = &ParseResult ($H[$i],$vremark,$censor,$F[$i]);
	    $F[$i] = $result;	    # replace result onto row for XLATE
		if ($result ne "") {
		    $NMeas{$H[$i]} ++;
		    if ($result != 0) {$NDet{$H[$i]} ++;$bin=1; }
		    else {$bin = 0;}
		}
		else {$bin = "_";}
	    ###########################################
	    push (@Compounds,$CompoundName{$H[$i]} || $H[$i]) if $result > 0;
	    ###########################################
	    $BinList .= $bin;
	}
	print XLATE join("\t",@F),"\n"; # Save translated input
    }
    push (@SampleHitArray, $BinList);
# ==================================================================
# -- Determine mixtures in the sample
# ==================================================================
    $nhit = $#Compounds + 1;
    $carry = join("\t",@Carry);
    @Compounds = sort (@Compounds);
    print $nhit, " detections in sample ",join(" ",@Carry),"\n"
	if $DEBUG;
    if    ($nhit <= 1) {next input;} # no mixtures
    elsif ($nhit <= $limit) {
	# No special hoops necessary, enumerate the mixtures
	($sample,$numix) = &combinate (2, $#Compounds+1, $carry, $TempDir, @Compounds);
	push (@Merges,$sample);
    }
    elsif ($nhit <= $limit*2) {
	# Employ a dual enumeration with temp-files method
	# enumerate mixtures in the first group of detections
        # opportunity for performance improvement through parallelism here?
        # if so, then split into two equal size groups 
	($file1,$numix) = &combinate (2, $limit, "", $TempDir, @Compounds[0..$limit-1]);
        @files = ($file1);
	
	# enumerate mixtures in the second group of detections
	($file2,$numix) = &combinate (2, $limit, "", $TempDir, @Compounds[$limit..$#Compounds]);
        push (@files,$file2) if $numix > 0;

        my $sortext = Sort::External->new (working_dir => $TempDir);
        				       
	tie *ONE, 'IO::Zlib', $file1, "rb" ;
	while ($one = <ONE>) {
	    # combine each single detection in group 2 with group 1 mixtures
	    chop $one;
	    foreach $cmpd (@Compounds[$limit..$#Compounds]) { 
		$sortext->feed(join("+",$one,$cmpd));
	    }
	    # list every combination of mixtures in groups 1 and 2
	    tie *TWO, 'IO::Zlib', $file2, "rb";
	    while ($two = <TWO>) {
                chop $two;
                $sortext->feed(join("+",$one, $two));
            }
	    untie *TWO;
	}

	# swap the next two loops to achieve sorting by mixture name?
	tie *TWO, 'IO::Zlib', $file2, "rb";
	while ($two = <TWO>) {
	    # combine each single detection in group 1 with group 2 mixtures
	    chop $two;
	    foreach $one (@Compounds[0..$limit-1]) { 
		$sortext->feed(join("+",$one,$two));
	    }
	}
	untie *ONE;
	untie *TWO;

	# enumerate all binary mixtures between groups 1 and 2
	foreach $cmpda (@Compounds[0..$limit-1]) { 
	    foreach $cmpdb (@Compounds[$limit..$#Compounds]) { 
		$sortext->feed(join("+",$cmpda,$cmpdb));
	    }
	}
	
        ($MERGE, $file3) = tempfile( "mergeXXXXXX", DIR => $TempDir, 
                                UNLINK => 0, SUFFIX => ".gz");
	tie *THREE, 'IO::Zlib', $file3, "wb";
        close $MERGE;
        $sortext->finish;
        while ( defined( $_ = $sortext->fetch ) ) { print THREE $_,"\n";}
        untie *THREE;

        ($MERGE, $sample) = tempfile( "sampleXXXXXX", DIR => $TempDir,
                                UNLINK => 0, SUFFIX => ".gz");
	tie *SAMPLE, 'IO::Zlib', $sample, "wb";
        close $MERGE;

        push(@files, $file3);
        $reader = sub { $_=shift;chop;return $_;};        
        my $Merger = File::MergeSort->new(\@files, $reader,{skip_empty_lines=>1});
        while ( my $row = $Merger->next_line() ) {
            chop $row;
            print SAMPLE join("\t",$row,$carry),"\n";
        }
        untie *SAMPLE;
        foreach $file (@files) {unlink $file;}
	push (@Merges, $sample);
    }
    else {
	print STDERR "Too many compounds ($nhit) for sample $carry\n";
	exit;
    }
		
    if ($Clock) {
	$endtime = times;
	$elapsed = $endtime - $starttime;
	print TIMING join("\t",$nhit,$elapsed),"\n" if $elapsed > 0;
    }
    print "$sample file created from $nhit detection sample\n" if $DEBUG;  
}

# ==================================================================
# -- Produce the compound-detection frequency report
# ==================================================================
    $samples = $row - 2;
    print REPORT "Samples ",$samples,"$/";
    print REPORT "Compounds ",$endcol - $startcol + 1,"$/$/";
    printf REPORT ("%-30s%12s%13s%12s%12s%s",
		   "Compound","Detected","Not Detected","Missing","Frequency%",$/);
    
    for ($i=$startcol; $i<=$endcol; $i++) {
	$name = $CompoundName{$H[$i]} || $H[$i];
	$Freq{$name} = "--";
	if ($NMeas{$H[$i]} != 0) {
	    $Freq{$name} = $NDet{$H[$i]} / $NMeas{$H[$i]}; }

	printf REPORT ("%-30s%12.0f%13.0f%12.0f%12.6f%s",
		       $name,
		       $NDet {$H[$i]},
		       $NMeas{$H[$i]} - $NDet{$H[$i]},
		       $samples - $NMeas{$H[$i]},
		       ( 100. * $Freq{$name} )),
                       $/;
    }

$elapsed = times - $liststarttime;
$starttime = times;
print REPORT "$/Potential mixtures of $NCompounds constituents: $pot_mixtures$/";
print REPORT "Mixture-listing phase completed in $elapsed seconds.$/";

# ==================================================================
# -- Merge mixtures identified for all the samples (separate files)
#    and, Count the frequency of each unique mixture occurrence
# ==================================================================

($COUNT, $CountFile) = tempfile( "countXXXXXX", DIR => $TempDir);
print REPORT "Established tempfile for counting$/" if $DEBUG;
$previous = "";
$count = 0;
$uniq_mixtures = 0;

$mreader = sub { $_=shift; my @a = split("\t",$_);return $a[0];};        
#print "Files to merge are: ",join(", ",@Merges),"$/" if $DEBUG;
print "Number of Files to merge is ",($#Merges+1),"$/" if $DEBUG;

# If the number of files is large, then the UNIX ulimit may need to increase
 
my $Merger = File::MergeSort->new( \@Merges, $mreader,{skip_empty_lines=>1});

print "Merging...$/" if $DEBUG;
while ( my $row = $Merger->next_line() ) {
    print RAW $row;
    ($mix,$other) = split("\t",$row,2);
    if ($previous ne "" && $previous ne $mix) {
	&reducemix;
	$previous = $mix;
	$count = 1;
	$uniq_mixtures ++;
    } else {$count ++;}
    $previous = $mix;
}
&reducemix;
$uniq_mixtures ++;

untie *RAW;
foreach $file (@Merges) {unlink $file;}

print REPORT "Finished counting unique mixtures$/" if $DEBUG;
print        "Finished counting unique mixtures$/" if $DEBUG;
close $COUNT;

# ==================================================================
# -- Sort the mixture-occurrence frequencies, most frequent first
# ==================================================================
$CountFile = basename($CountFile);
print REPORT "Temp directory is $TempDir$/" if $DEBUG;

$SortProg = catfile($ENV{QWMIXDIR},"qwmixsort.pl");
$TEMP = new IO::Pipe;
$TEMP -> reader ("$SortProg $TempDir $CountFile 2nr");

print REPORT "Preparing to sort unique mixtures$/" if $DEBUG;
# ==================================================================
# -- Create the output mixture-frequency file, with column headings
#    listing unique mixture, occurrence count, expected & actual
#    frequency
# ==================================================================

tie *MIX, 'IO::Zlib', catfile($Directory,"mixtures.rdb.gz"), "wb";

print MIX join("\t","mixture","count","expected%","actual%"),"\n";
print MIX join("\t",$MixDef,"7n","9n","9n"),"\n";
$listedmix = 0;

while (<$TEMP>)  {print MIX;$listedmix ++;}
close $TEMP;
untie *MIX; # close
# ==================================================================
# -- Print summary messages with timestamp
# ==================================================================
$elapsed = times - $starttime;
@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 "Mixture-enumeration phase completed in $elapsed seconds.$/";
print REPORT "$uniq_mixtures unique mixtures identified in data file $InFile.$/";
print REPORT "$listedmix unique mixtures listed that occurred at least $mincnt times.$/";
print REPORT "Processing completed at $DateTime$/";
close REPORT;

# Screen
print "$uniq_mixtures unique mixtures identified.\n";
print "$listedmix unique mixtures listed that occurred at least $mincnt times.\n";
unlink catfile($TempDir, $CountFile);
