text:
//          Copyright Gavin Band 2008 - 2012.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

/*
 * This program, qctool, selects rows from a GEN file according to certain criteria.
 * - rows where some genotype data is missing, or none is.
 * - rows where genotype data is, or is not in hardy-weinberg equilibrium.
 * - rows with SNP IDs in a given list.
 *
 * Program arguments:
 *
 */

#include <string>
#include <map>
#include <set>
#include <fstream>
#include <memory>
#include <numeric>
#include <algorithm>
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>
#include <boost/iterator/counting_iterator.hpp>
#include <boost/filesystem.hpp>
#include <boost/uuid/uuid.hpp>
#include <boost/uuid/random_generator.hpp>
#include <boost/uuid/uuid_io.hpp>

#include "appcontext/appcontext.hpp"

#include "genfile/db/Error.hpp"
#include "genfile/SNPDataSource.hpp"
#include "genfile/SNPDataSourceChain.hpp"
#include "genfile/SNPDataSourceRack.hpp"
#include "genfile/MergingSNPDataSource.hpp"
#include "genfile/SampleMappingSNPDataSource.hpp"
#include "genfile/SNPDataSinkChain.hpp"
#include "genfile/GenFileSNPDataSink.hpp"
#include "genfile/VCFFormatSNPDataSink.hpp"
#include "genfile/SortingBGenFileSNPDataSink.hpp"
#include "genfile/TrivialSNPDataSink.hpp"
#include "genfile/CategoricalCohortIndividualSource.hpp"
#include "genfile/CountingCohortIndividualSource.hpp"
#include "genfile/SampleFilteringCohortIndividualSource.hpp"
#include "genfile/CohortIndividualSourceChain.hpp"
#include "genfile/DataLinkingCohortIndividualSource.hpp"
#include "genfile/SampleFilteringSNPDataSource.hpp"
#include "genfile/SNPTranslatingSNPDataSource.hpp"
#include "genfile/StrandAligningSNPDataSource.hpp"
#include "genfile/ThreshholdingSNPDataSource.hpp"
#include "genfile/VCFFormatSNPDataSource.hpp"
#include "genfile/PloidyConvertingSNPDataSource.hpp"
#include "genfile/CommonSNPFilter.hpp"
#include "genfile/SNPFilteringSNPDataSource.hpp"
#include "genfile/ReorderingSNPDataSource.hpp"
#include "genfile/ReorderingCohortIndividualSource.hpp"
#include "genfile/VariantIdentifyingDataFilteringSNPDataSource.hpp"
#include "genfile/utility.hpp"
#include "genfile/QuantileNormalisingCrossCohortCovariateValueMapping.hpp"
#include "genfile/ValueMappingCohortIndividualSource.hpp"
#include "genfile/vcf/StrictMetadataParser.hpp"
#include "genfile/WithSNPDosagesCohortIndividualSource.hpp"
#include "genfile/get_set_eigen.hpp"
#include "genfile/SampleFilter.hpp"
#include "genfile/CompoundSampleFilter.hpp"
#include "genfile/SampleFilterNegation.hpp"
#include "genfile/VariableInSetSampleFilter.hpp"
#include "genfile/VariableInRangeSampleFilter.hpp"
#include "genfile/get_ploidy_from_sex.hpp"

#include "statfile/BuiltInTypeStatSource.hpp"
#include "statfile/from_string.hpp"
#include "statfile/read_values.hpp"
#include "statfile/SNPDataSourceAdapter.hpp"

#include "qcdb/FlatFileOutputter.hpp"
#include "qcdb/FlatTableDBOutputter.hpp"

#include "worker/QueuedMultiThreadedWorker.hpp"
#include "worker/SynchronousWorker.hpp"

#include "InputToOutputFilenameMapper.hpp"
#include "SNPIDSink.hpp"
#include "statfile/BuiltInTypeStatSink.hpp"
#include "DataReadTest.hpp"
#include "components/RelatednessComponent/RelatednessComponent.hpp"
#include "components/HaplotypeFrequencyComponent/HaplotypeFrequencyComponent.hpp"
#include "components/SNPSummaryComponent/SNPSummaryComponent.hpp"
#include "components/SNPOutputComponent/SNPOutputComponent.hpp"
#include "components/SNPOutputComponent/SQLiteHaplotypesSNPDataSink.hpp"
#include "components/SNPOutputComponent/SQLiteGenotypesSNPDataSink.hpp"
#include "components/SampleSummaryComponent/SampleSummaryComponent.hpp"
#include "components/SampleSummaryComponent/FlatTableDBOutputter.hpp"
#include "components/SampleSummaryComponent/FlatFileOutputter.hpp"

#include "config/package_revision_autogenerated.hpp"

namespace globals {
	std::string const program_name = "qctool" ;
	std::string const program_version = package_version ;
	std::string const program_revision =  std::string( package_revision ).substr( 0, 7 ) ;
}

namespace impl {
	void set_vector_entry( std::vector< genfile::VariantEntry >* vector, std::size_t i, std::string const& value ) {
		assert( i < vector->size() ) ;
		(*vector)[i] = value ;
	}

	void set_vector_entry_string( std::vector< std::string >* vector, std::size_t i, std::string const& value ) {
		assert( i < vector->size() ) ;
		(*vector)[i] = value ;
	}

	void set_vector_entry_variant( std::vector< std::string >* vector, std::size_t i, genfile::VariantEntry const& value ) {
		assert( i < vector->size() ) ;
		(*vector)[i] = value.as< std::string >() ;
	}
}

struct QCToolOptionProcessor: public appcontext::CmdLineOptionProcessor
{
public:
	std::string get_program_name() const { return globals::program_name ; }

	void declare_options( appcontext::OptionProcessor& options ) {
		// Meta-options
		options.set_help_option( "-help" ) ;
		options.set_spec_option( "-spec" ) ;

		// File options
		options.declare_group( "Input file options" ) ;
	    options[ "-g" ]
	        .set_description( 	"Path of gen file(s) to input.  "
								"The given filename may contain the wildcard character '#', which expands to match a"
								"one- or two-character chromosome identifier.  (For example, \"qctool -g myfile_#.gen\" will match "
								"\"myfile_1.gen\", \"myfile_2.gen\", etc., or \"myfile_01.gen\", \"myfile_02.gen\", etc.)  Only human"
								" autosomes and sex chromosomes are matched this way.\n"
								"This option may also be repeated, in which case each invocation is treated as a seperate cohort and cohorts"
								" are joined together to create one big dataset." )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 1 )
			.set_maximum_multiplicity( 1000 ) ;

	    options[ "-filetype" ]
			.set_description(
				"Specify the filetype of the genotype files specified by -g. "
				"By default, qctool will guess the file type.  Use this option to override that guess. "
				"Possible types are: \"" + genfile::string_utils::join( genfile::SNPDataSource::get_file_types(), "\",\"" ) + "\"." )
			.set_takes_single_value()
			.set_default_value( "guess" ) ;

		options[ "-permissive" ]
			.set_description( "Specify that qctool should be permissive about problems when parsing input."
				"Currently this applies to VCF file input.  Entries that cannot be parsed are interpreted as missing values."
				"If this option is not supplied, qctool will fail on malformed files."
			) ;

	    options[ "-ofiletype" ]
			.set_description(
				"Specify the filetype of the output genotype files specified by -og. "
				"By default, qctool will guess the file type.  Use this option to override that guess. "
				"Possible types are: \"" + genfile::string_utils::join( genfile::SNPDataSink::get_file_types(), "\",\"" ) + "\"." )
			.set_takes_single_value()
			.set_default_value( "guess" ) ;

		options[ "-merge-in" ]
			.set_description( "Specify an additional set of genotypes that should be merged in (in position order) to the dataset. "
				"This must have the same number of samples as the data set for -g. "
				"Note that filtering, strand alignment, allele matching, or other transformations are not applied to the data "
				"specified by -merge-in.")
			.set_takes_values( 2 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 1 ) ;

		options[ "-merge-strategy" ]
			.set_description( "Specify a strategy to use when encountering SNPs with the same position in a merge. "
				"Options are \"" + genfile::string_utils::join( genfile::MergingSNPDataSource::get_merge_strategies(), "\",\"" ) + "\"." )
			.set_takes_single_value()
			.set_default_value( "keep-all" ) ;
		options[ "-merge-prefix" ]
			.set_description( "Specify a string to add as a prefix to ID fields of merged-in variants" )
			.set_takes_single_value()
			.set_default_value( "" ) ;
		options[ "-match-sample-ids" ]
			.set_description( "Specify the columns in the main and comparison dataset sample files that will be used to match samples. "
					"The value should be of the form <main dataset column>~<comparison dataset column>." )
			.set_takes_single_value()
			.set_default_value( "ID_1~ID_1" )
		;

	    options[ "-s" ]
	        .set_description( "Path of sample file to input.  If specified, this option must occur as often as the -g option"
							" to specify one sample file per cohort." )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 1000 ) ;

	    options[ "-sample-data" ]
	        .set_description( "Path of additional sample file(s) (in the same format as accepted by -s) to read." )
			.set_takes_values_until_next_option()
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;

		options[ "-sex-column" ]
			.set_description( "Specify the name of a column (of type D) in the sample files(s) containing information on sex." )
			.set_takes_values( 1 )
			.set_default_value( "sex" )
		;

		options[ "-infer-ploidy-from" ]
			.set_description(
				"Specify that QCTOOL should infer the true ploidy of genotype calls "
				"(and genotype call probabilities) "
				"based on the sex of samples and the chromosome. "
				"Specifically, this is appropriate for some human datasets; "
				"it replaces diploid genotype calls with haploid calls for males on the X and Y chromosomes. "
				"Genotypes are set to missing if the original diploid call is heterozygote. "
				"For genotype call probabilities, the probabilities of homozygote calls are used to infer the haploid "
				"genotype probabilities, while probabilities are set to missing if there is nonzero probability of a heterozygote. "
				"Females calls on the Y chromosome are also set to missing."
			)
			.set_takes_single_value() ;

		options[ "-haploid-genotype-coding" ]
			.set_description( "Specify whether haploid samples (e.g. males on the non-pseudoautosomal region of the X chromosome)"
			" are coded as heterozygotes (\"het\") or homozygotes (\"hom\")." )
			.set_takes_single_value()
			.set_default_value( "hom" )
		;

		options.option_implies_option( "-infer-ploidy-from", "-s" ) ;

/*	    options[ "-i" ]
	        .set_description( 	"Path of intensity file(s) to input.  "
								"If specified, this option must occur the same number of times as the -g option, to specify"
								" one intensity file per cohort."
								"As with -g the '#' wildcard character may be used to match chromosomes." )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 1 )
			.set_maximum_multiplicity( 100 ) ;
*/
		options.declare_group( "Sample exclusion options" ) ;
		options[ "-incl-samples"]
			.set_description( "Filter out samples whose sample ID does not lie in the given file(s).")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-samples"]
			.set_description( "Filter out samples whose sample ID lies in the given file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-samples-where"]
			.set_description( "Include samples by specifying conditions on the values of columns in the sample file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-samples-where"]
			.set_description( "Exclude samples by specifying conditions on the values of columns in the sample file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;

		// SNP exclusion options
		options.declare_group( "SNP exclusion options" ) ;
		options[ "-excl-variants" ]
			.set_description( "Exclude all SNPs in the given file(s) from the analysis."
				" The format of this file is the same as that output by the -write-snp-excl-list option. "
				" It must have six columns interpreted as SNPID, rsid, chromosome, position, first and second alleles." )
			.set_takes_values( 1 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-variants" ]
			.set_description( "Exclude all SNPs not in the given file(s) from the analysis."
				" The format of this file is the same as that output by the -write-snp-excl-list option. "
				" It must have six columns interpreted as SNPID, rsid, chromosome, position, first and second alleles." )
			.set_takes_values( 1 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-snpids" ]
			.set_description( "Exclude all SNPs whose SNPID is in the given file(s) from the analysis.")
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-rsids" ]
			.set_description( "Exclude all SNPs whose RSID is in the given file(s) from the analysis.")
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-snpids" ]
			.set_description( "Exclude all SNPs whose SNPID is not in the given file(s) from the analysis.")
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-rsids" ]
			.set_description( "Exclude all SNPs whose RSID is not in the given file(s) from the analysis.")
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-positions" ]
			.set_description( "Exclude all SNPs whose position is in the given file(s) from the analysis. "
				"Positions should be in the form [chromosome]:[position] and separated by whitespace." )
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-positions" ]
			.set_description( "Exclude all SNPs whose position is not in the given file(s) from the analysis. "
				"Positions should be in the form [chromosome]:[position] and separated by whitespace." )
			.set_takes_values_until_next_option()
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-variants-matching" ]
			.set_description( "Filter out snps whose rsid or SNPID matches the given value. "
				"The value should be a string which can contain a % wildcard character (which matches any substring). "
				"Optionally, prefix the argument with snpid~ or rsid~ to only match against the SNPID or rsid fields." )
			.set_takes_values_until_next_option() ;
		options[ "-incl-variants-matching" ]
			.set_description( "Filter out snps whose rsid or SNPID does not match the given value. "
				"The value should be a string which can contain a % wildcard character (which matches any substring). "
				"Optionally, prefix the argument with snpid~ or rsid~ to only match against the SNPID or rsid fields." )
			.set_takes_values_until_next_option() ;
		options[ "-incl-range" ]
			.set_description( "Specify a range of SNPs (or comma-separated list of ranges of SNPs) to operate on. "
				"Each range should be in the format CC:xxxx-yyyy where CC is the chromosome and xxxx and yyyy are the "
				"start and end coordinates, or just xxxx-yyyy which matches that range from all chromosomes. "
				"You can also omit either of xxxx or yyyy to get all SNPs from the start or to the end of a chromosome." )
			.set_takes_values_until_next_option() ;
		options[ "-excl-range" ]
			.set_description( "Specify a range of SNPs (or comma-separated list of ranges of SNPs) to exclude from operation. "
				"Each range should be in the format CC:xxxx-yyyy where CC is the chromosome and xxxx and yyyy are the "
				"start and end coordinates, or just xxxx-yyyy which matches that range from all chromosomes. "
				"You can also omit either of xxxx or yyyy to get all SNPs from the start or to the end of a chromosome." )
			.set_takes_values_until_next_option() ;
		options[ "-incl-ranges" ]
			.set_description( "Like -incl-range, but read the ranges from the specified file."
				"The file must contain a whitespace-separated list of genomic ranges to include in the analysis." )
			.set_takes_values_until_next_option() ;
		options[ "-excl-ranges" ]
			.set_description( "Like -excl-range, but read the ranges from the specified file."
				"The file must contain a whitespace-separated list of genomic ranges to exclude from the analysis." )
			.set_takes_values_until_next_option() ;

		options.declare_group( "Options for adjusting SNPs" ) ;
	    options[ "-strand" ]
	        .set_description(
				"Path of strand file(s) to input.  "
				"If specified, this option must occur the same number of times as the -g option, to specify"
				" one intensity file per cohort. "
				"This file can have SNPID, rsid, chromosome, position, alleleA, alleleB, strand columns. "
				"Strand should be specified as +, - or NA.  Alleles for variants with - strand will be complemented. "
				"Variants with strand=NA, or with missing strand information, will be skipped entirely."
			)
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
	    options[ "-flip-to-match-allele" ]
	        .set_description( "Specify that we should flip data if necessary so that the first allele is the given allele given in the "
				"specified column of the strand file (which must be supplied using the -strand option).  First the strand is applied, "
				"then the two alleles are compared to the given allele and data flipped if necessary."
			)
			.set_takes_values( 1 )
		;
		options.option_implies_option( "-flip-to-match-allele", "-strand" ) ;
		options[ "-map-id-data" ]
			.set_description( "Update the chromosome, position, IDs and/or alleles of a set of SNPs with new values."
				" The argument must be a file with six named columns giving the original SNPID, rsid, chromosome, position and alleles,"
				" followed by another six columns containing the values to replace with."
				" SNPs not in this file will be passed to the output file unchanged. "
				" This option only affects the identifying data, not genotypes themselves."
			)
			.set_takes_single_value() ;
		options[ "-flip-to-match-cohort1" ]
			.set_description( "Specify that alleles (and corresponding genotypes) in all cohorts should be switched, if necessary,"
				" so as to match the alleles of the first cohort."
			) ;
		options[ "-compare-variants-by" ]
			.set_description( "By default, matching SNPs between cohorts uses all the available fields"
				" (position, rsid, snpid, and alleles.)"
				" Use this option to specify a comma-separated subset of those fields to use instead."
				" The first entry must be \"position\"."
				" This option can be used, for example, when cohorts are typed on different platforms so have different SNPID fields." )
			.set_takes_single_value()
			.set_default_value( "position,alleles,ids" ) ;
		options[ "-assume-chromosome" ]
			.set_description( "Treat each SNP whose chromosome cannot be determined"
				" as though it lies on the specified chromosome." )
			.set_takes_single_value() ;

		options.declare_group( "Output file options" ) ;
	    options[ "-og" ]
	        .set_description(
				"Specify that qctool should write an output genotype file with the specified filename. "
				"The type of this file will be determined from the filename extension. "
				"If the first occurence of -g uses a '#' wildcard character, the '#' character can "
				"also be used here to specify numbered output files corresponding to the input files."
				"If \"-\" is specifed here, genotypes are written to standard output in GEN format."
			)
	        .set_takes_values( 1 )
			.set_maximum_multiplicity( 1 ) ;
		options[ "-sort" ]
			.set_description( "Sort the genotypes in the output file.  Currently this is only supported if BGEN, unzipped GEN, unzipped VCF format is output." ) ;
		options[ "-os" ]
	        .set_description( "Output sample information to the file specified.  " )
	        .set_takes_single_value() ;
		options[ "-output-sample-format" ]
	        .set_description( "Specify the format of sample files generated by QCTOOL.  Options are \"traditional\" for the"
	 		" traditional format with column types specified on the second line, or \"new\" for a new format with column types"
			" specified in metadata encoded as JSON in the header." )
	        .set_takes_single_value()
			.set_default_value( "traditional" ) ;

		options[ "-osnp" ]
			.set_description( "Set the name of the file used to output results of per-SNP computations." )
			.set_takes_single_value() ;

		options[ "-osample" ]
			.set_description( "Set the name of the file used to output results of per-sample computations." )
			.set_takes_single_value() ;

		options[ "-omit-chromosome" ]
			.set_description( "(This option is specific to output files in the GEN format.) Do not output a chromosome column." ) ;

		options[ "-no-clobber" ]
			.set_description( "Do not run QCTOOL if it would overwrite an existing file." )
		;

		options.declare_group( "GEN-specific options" ) ;
		options[ "-precision" ]
			.set_description( "Currently used when outputting GEN files only.  Use this many bits of precision in output." )
			.set_default_value( 5 )
			.set_takes_single_value()
		;
		options.declare_group( "BGEN-specific options" ) ;
		options[ "-bgen-bits" ]
			.set_description( "For use when outputting BGEN files only.  Tell QCTOOL to use this number"
				" of bits to store each probability." )
			.set_default_value( 16 )
			.set_takes_single_value()
		;
		options[ "-bgen-free-data" ]
			.set_description( "For use when outputting BGEN files only.  Tell QCTOOL to write the specified free data segment"
				" in the BGEN header block." )
			.set_takes_single_value()
		;
		options[ "-bgen-compression" ]
			.set_description( "Specify what compression to use when outputting BGEN files only."
				" This can be \"none\", \"zlib\", or \"zstd\"." )
			.set_default_value( "zlib" )
			.set_takes_single_value()
		;
		options[ "-bgen-permitted-input-rounding-error" ]
			.set_description(
				"Specify the maximum error that will be tolerated in input probability values when writing a BGEN file. "
				"E.g. if data is accurate to 3 decimal places, a value of 0.0005 should work here, while if "
				"your data is less accurate you may need to relax this.  BGEN v1.2 format rescales genotype "
				"probabilities so that they sum to 1; this implies "
				"a scaling factor of up to (1 ± (n × error)), where n is the number of "
				"probability values for the sample. "
				"(BGEN encoding also incurs additional error of up to 1/(2^bits-1) per probability)."
			)
			.set_takes_single_value()
			.set_default_value( 0.001 )
		;
		options[ "-bgen-omit-sample-identifier-block" ]
			.set_description( "For use when outputting BGEN files only.  Tell QCTOOL to omit the sample identifier block.  By default"
				" this is written whenever -s is specified." )
		;

		options.declare_group( "sqlite-specific options" ) ;
		options[ "-sqlite-dont-link-variants" ]
			.set_description( "Do not attempt to link to existing variants when outputting to a sqlite file. " )
		;
		options[ "-sqlite-store-alt-identifiers" ]
			.set_description( "Store variant alternate identifiers. This adds information but can be slower." )
		;

		options.option_implies_option( "-sort", "-og" ) ;
		options.option_implies_option( "-omit-chromosome", "-og" ) ;
		options.option_implies_option( "-output-sample-format", "-os" ) ;

		// VCF file options
		options.declare_group( "VCF file options" ) ;
		options[ "-vcf-genotype-field" ]
			.set_description(
				"Specify the name of the field in a VCF file to read genotypes from.  This must match "
				"the name of a FORMAT field in the VCF file."
			)
			.set_takes_single_value()
			.set_default_value( "GT" ) ;
		options[ "-vcf-intensity-field" ]
			.set_description(
				"Specify the name of the field in a VCF file to read intensities from.  This must match "
				"the name of a FORMAT field in the VCF file."
			)
			.set_takes_single_value()
			.set_default_value( "XY" ) ;
		options[ "-vcf-output-fields" ]
			.set_description(
				"Specify a subset of fields to appear in output vcf files"
			)
			.set_takes_values_until_next_option()
		;

		options[ "-metadata" ]
			.set_description(
				"Specify the name of a file containing VCF metadata to be used to parse "
				"a VCF file.  Keys in this file must either not be present in the VCF file, or must have "
				"identical values."
			)
			.set_takes_single_value() ;

		// SNP filtering options
		options.declare_group( "SNP filtering options" ) ;

		// Sample filtering options
		options.declare_group( "Sample filtering options" ) ;

		// Inclusion / exclusion list options
		options.declare_group( "Inclusion / exclusion list options" ) ;
		options[ "-write-sample-excl-list" ]
			.set_description( "Write a file containing a list of the ids of individuals which are filtered out by the sample filter." )
			.set_takes_single_value() ;
		options [ "-write-snp-excl-list" ]
			.set_description( "Write a file (or files) containing the SNPs that are filtered out." )
			.set_takes_single_value() ;
		;

		// Other options
		options.declare_group( "Other options" ) ;
		options [ "-threshold" ]
			.set_description( "Threshhold genotype call probabilities.  This means replace call probabilities "
			">= threshold with 1 and those < threshhold with 0." )
			.set_takes_single_value() ;
		options [ "-force" ]
			.set_description( "Ignore warnings and proceed with requested action." ) ;
		options [ "-log" ]
			.set_description( "Specify that " + globals::program_name + " should write a log file to the given file." )
			.set_takes_single_value() ;
		options [ "-threads" ]
			.set_description( "Specify the number of worker threads to use in computationally intensive tasks." )
			.set_takes_single_value()
			.set_default_value( 0 ) ;
		options[ "-analysis-name" ]
			.set_description( "Specify a name to label results from this analysis with." )
			.set_takes_single_value()
			.set_default_value( "qctool analysis" ) ;
		options[ "-analysis-chunk" ]
			.set_description( "Specify a name denoting the current genomic region or chunk on which this is run.  This is intended for use in parallel environments." )
			.set_takes_single_value()
			.set_default_value( genfile::MissingValue() ) ;

		options.declare_group( "Options for adjusting sample data" ) ;
		options[ "-quantile-normalise" ]
			.set_description( "Quantile normalise each specified continuous phenotype or covariate "
			 	"by ranking its values and mapping to quantiles of the standard normal distribution N(0,1). "
				"Ties are handled by sending tied values to the average of the corresponding quantiles."
				"The argument should be a comma-separated list of column names from the sample file." )
			.set_takes_values_until_next_option() ;
		options[ "-missing-code" ]
			.set_description( "Specify a comma-separated list of strings to be treated as missing values "
				"when encountered in the sample file(s)." )
			.set_takes_single_value()
			.set_default_value( "NA" ) ;
		options[ "-condition-on" ]
			.set_description( "Condition on the dosages of a given SNP or set of SNPs."
				" The argument should be a comma-separated list of values of the form:\n"
				"   [field]~[value]:[dose][,dose...]\n"
				"where field is \"rsid\", \"snpid\", or \"pos\", value is the value to match, and "
				"each dose is one of add,dom,het,rec,or gen."
				" If the field is omitted, it is assumed to be rsid;"
				" if the dose is omitted it is assumed to be add."
			)
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;

		options [ "-reorder" ]
			.set_description( "Specify that qctool should write output files with individuals re-ordered. "
				"The argument can be one of three things:\n"
				" - the name of a readable file - in this case samples IDs are read from the file and output in the order specified.\n"
				" - or the word \"randomly\" - in this case samples are randomly reordered.\n"
				" - or the word \"backwards\" - in this case the order of samples is reveresed."
			)
			.set_takes_single_value() ;

		options.option_implies_option( "-quantile-normalise", "-s" ) ;
		options.option_implies_option( "-quantile-normalise", "-os" ) ;
		options.option_implies_option( "-condition-on", "-s" ) ;
		options.option_implies_option( "-condition-on", "-os" ) ;

		//Relatotron::declare_options( options ) ;
		RelatednessComponent::declare_options( options ) ;
		DataReadTest::declare_options( options ) ;
//		ClusterFitter::declare_options( options ) ;
		SNPSummaryComponent::declare_options( options ) ;
		SampleSummaryComponent::declare_options( options ) ;
		HaplotypeFrequencyComponent::declare_options( options ) ;

		options.option_excludes_option( "-write-sample-excl-list", "-os" ) ;

		options.option_excludes_option( "-write-snp-excl-list", "-og" ) ;

		options.option_implies_option( "-excl-samples-where", "-s" ) ;
		options.option_implies_option( "-incl-samples-where", "-s" ) ;
	}
} ;

struct QCToolOptionMangler {
	QCToolOptionMangler( appcontext::OptionProcessor const& options ):
		m_options( options )
	{
		process_filenames() ;
	}
	std::vector< std::string > const& input_sample_filenames() const { return m_input_sample_filenames ; }
	std::string const& output_sample_filename() const { return m_output_sample_filename ; }
	std::string const& output_sample_excl_list_filename() const { return m_output_sample_excl_list_filename ; }
	std::string const log_filename() const { return m_options.get_value< std::string > ( "-log" ) ; }
	InputToOutputFilenameMapper const& gen_filename_mapper() const { return m_gen_file_mapper ; }
	std::vector< std::vector< genfile::wildcard::FilenameMatch > > const& gen_filenames() const { return m_gen_filenames ; }
	InputToOutputFilenameMapper const& snp_excl_list_filename_mapper() const { return m_output_snp_excl_file_mapper ; }
private:

	void process_filenames() {
		get_snp_related_filenames() ;
		get_sample_related_filenames() ;
	}

	void get_snp_related_filenames() {
		assert( m_options.check_if_option_was_supplied( "-g" )) ;
		std::vector< std::string > input_gen_filenames_supplied = m_options.get_values< std::string >( "-g" ) ;
		assert( input_gen_filenames_supplied.size() >= 1 ) ;
		std::string
			output_gen_filename = construct_output_gen_filename( input_gen_filenames_supplied ),
			output_snp_excl_filename = construct_output_snp_excl_list_filename( input_gen_filenames_supplied ) ;

		m_gen_file_mapper.add_filename_pair( input_gen_filenames_supplied[0], output_gen_filename ) ;
		m_output_snp_excl_file_mapper.add_filename_pair( input_gen_filenames_supplied[0], output_snp_excl_filename ) ;

		m_gen_filenames.resize( input_gen_filenames_supplied.size() ) ;
		for( std::size_t i = 0; i < input_gen_filenames_supplied.size(); ++i ) {
			m_gen_filenames[i] = genfile::wildcard::find_files_by_chromosome(
				input_gen_filenames_supplied[i],
				genfile::wildcard::eNON_SEX_CHROMOSOMES
			) ;
		}
	}

	void get_sample_related_filenames() {
		if( m_options.check_if_option_was_supplied( "-s" ) ) {
			m_input_sample_filenames = m_options.get_values< std::string >( "-s" ) ;
			if( m_input_sample_filenames.size() != m_gen_filenames.size() ) {
				throw genfile::BadArgumentError(
					"QCToolOptionMangler::get_sample_related_filenames()",
					"-s \"" + genfile::string_utils::join( m_input_sample_filenames, " " ) + "\"",
					"Number of sample files must match number of genotype files (" + genfile::string_utils::to_string( m_options.get_values< std::string >( "-s" ) .size() ) + ") specified using -g."
				) ;
			}
		}

		// Otherwise, we need to write a sample exclusion list file if -write-sample-excl-list was given.
		else if( m_options.check_if_option_was_supplied( "-write-sample-excl-list" )) {
			m_output_sample_excl_list_filename = m_options.get_value< std::string > ( "-write-sample-excl-list" ) ;
		}

		// We need to write a sample file if -os is given:
		// -os is given
		if( m_options.check_if_option_was_supplied( "-os" )) {
			m_output_sample_filename = m_options.get_value< std::string >( "-os" ) ;
		}
	}

	std::string strip_sample_file_extension_if_present( std::string filename ) {
		if( filename.size() >= 7 && filename.substr( filename.size() - 7, 7 ) == ".sample" ) {
			filename.resize( filename.size() - 7 ) ;
		}
		return filename ;
	}

	std::string construct_output_gen_filename( std::vector< std::string > const& input_gen_filenames_supplied ) {
		std::string result ;
		// We need to produce output gen files if
		//  * -og is given
		//  OR
		//  ** -write-snp-excl-list is not set
		//  * AND EITHER
		//  ** some sample filters are given (but not -write-sample-excl-list)
		//  ** OR some SNP filters are given.
		if( m_options.check_if_option_was_supplied( "-og" )
			|| (
				!( m_options.check_if_option_was_supplied( "-write-snp-excl-list" ) || m_options.check_if_option_was_supplied( "-write-snp-incl-list" ) )
				&& (
					(m_options.check_if_option_was_supplied_in_group( "Sample filtering options" ) && !m_options.check_if_option_was_supplied( "-write-sample-excl-list" ))
					||
					m_options.check_if_option_was_supplied_in_group( "SNP filtering options" )
				)
			)
		) {
			if( m_options.check_if_option_was_supplied( "-og" ) ) {
				result = m_options.get_value< std::string >( "-og" ) ;
			} else {
				assert( input_gen_filenames_supplied.size() > 0 ) ;
				std::string stub ;
				if( input_gen_filenames_supplied.size() == 1 )  {
					stub = input_gen_filenames_supplied[0] ;
				}
				else {
					stub = "qctool_cohorts_1-" + genfile::string_utils::to_string( input_gen_filenames_supplied.size() ) ;
				}
				result
					= genfile::strip_gen_file_extension_if_present( stub )
					+ ".fltrd"
					+ genfile::get_gen_file_extension_if_present( input_gen_filenames_supplied[0] ) ;
			}
		}
		return result ;
	}

	std::string construct_output_snp_excl_list_filename( std::vector< std::string > const& input_gen_filenames_supplied ) {
		std::string result ;
		if( m_options.check_if_option_was_supplied( "-write-snp-excl-list" ) ) {
				result = m_options.get_value< std::string >( "-write-snp-excl-list" ) ;
		}
		return result ;
	}

private:
	appcontext::OptionProcessor const& m_options ;
	std::vector< std::string > m_input_sample_filenames ;
	std::string m_output_sample_filename ;
	std::string m_output_sample_excl_list_filename ;
	std::string m_log_filename ;
	InputToOutputFilenameMapper m_gen_file_mapper ;
	InputToOutputFilenameMapper m_output_snp_excl_file_mapper ;
	std::vector< std::vector< genfile::wildcard::FilenameMatch > > m_gen_filenames ;
} ;

struct QCToolCmdLineContext
{
	typedef genfile::SNPDataSource SNPDataSource ;
	typedef genfile::SNPDataSink SNPDataSink ;

	QCToolCmdLineContext( appcontext::OptionProcessor const& options, appcontext::UIContext& ui_context ):
		m_options( options ),
		m_mangled_options( options ),
		m_ui_context( ui_context )
	{
		try {
			setup() ;
		}
		catch( genfile::MalformedInputError const& e ) {
			m_ui_context.logger() << "\nError (" << e.what() <<"): " << e.format_message() << ".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::FileContainsSNPsOfDifferentSizes const& ) {
			m_ui_context.logger() << "\nError: The GEN files specified did not all have the same sample size.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::FileNotFoundError const& e ) {
			m_ui_context.logger() << "\nError: No file matching \"" << e.filespec() << "\" could be found.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::ResourceNotOpenedError const& e ) {
			m_ui_context.logger() << "\nError: (" << e.what() << "): resource \"" << e.source() << "\" could not be opened.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch ( genfile::KeyNotFoundError const& e ) {
			m_ui_context.logger() << "\nKey \"" << e.key() << "\" was not found in source \"" << e.source() << "\".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch ( genfile::DuplicateKeyError const& e ) {
			m_ui_context.logger() << "\nError (" << e.what() << "): Duplicate key \"" << e.key() << "\" was found in source \"" << e.source() << "\".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( statfile::FileNotOpenedError const& e ) {
			m_ui_context.logger() << "\nError: No file matching \"" << e.filename() << "\" could be found.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::OperationFailedError const& e ) {
			m_ui_context.logger() << "\nError (" << e.what() << "): " << e.caller() << ": error attempting to perform the following operation on " << e.object()<< ":\n"
				<< e.operation() << ".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::OperationUnsupportedError const& e ) {
			m_ui_context.logger() << "\nError (" << e.what() << "): " << ": the following operation on " << e.object()<< " is not supported:\n"
				<< e.operation() << ".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
	}

	~QCToolCmdLineContext() {
		write_postamble() ;
		// Ensure a fixed deconstruction order so that samples and pedigree outlive
		// the SNP Data Sinks.
		m_fltrd_out_snp_data_sink.reset() ;
		m_fltrd_in_snp_data_sink.reset() ;
		m_samples.reset() ;
	}

	SNPDataSource& snp_data_source() const {
		return *m_snp_data_source ;
	}

	SNPDataSink& fltrd_in_snp_data_sink() const {
		return *m_fltrd_in_snp_data_sink ;
	}

	SNPDataSink& fltrd_out_snp_data_sink() const {
		return *m_fltrd_out_snp_data_sink ;
	}

	std::vector< std::size_t > const& indices_of_filtered_out_samples() const { return m_indices_of_filtered_out_samples ; }

	genfile::CohortIndividualSource const& samples() const { return *m_samples ; }

	bool ignore_warnings() const {
		return m_ignore_warnings ;
	}

	void write_preamble() {
		m_ui_context.logger() << std::string( 72, '=' ) << "\n\n" ;

		m_ui_context.logger() << std::setw(30) << "Input SAMPLE file(s):" ;
		for( std::size_t i = 0; i < m_mangled_options.input_sample_filenames().size(); ++i ) {
			if( i > 0 ) {
				m_ui_context.logger() << std::setw(30) << "" ;
			}
			m_ui_context.logger() << "  \"" << format_filename( m_mangled_options.input_sample_filenames()[i]) << "\"\n" ;
		}

		if(
			genfile::ValueMappingCohortIndividualSource const* source
				= dynamic_cast< genfile::ValueMappingCohortIndividualSource const* >( m_samples.get() )
		) {
			m_ui_context.logger() << std::setw(30) << "Sample file summary:" << "  " ;
			m_ui_context.logger() << source->get_summary( std::string( 32, ' ' )) << "\n" ;
		}

		m_ui_context.logger() << std::setw(30) << "Output SAMPLE file:"
			<< "  \"" << format_filename( m_mangled_options.output_sample_filename()) << "\".\n" ;
		if( m_options.check( "-osample" )) {
			m_ui_context.logger() << std::setw(30) << "Sample statistic output file:"
				<< "  \"" << format_filename( m_options.get< std::string >( "-osample" )) << "\".\n" ;
		}
		m_ui_context.logger() << std::setw(30) << "Sample exclusion output file:"
			<< "  \"" << format_filename( m_mangled_options.output_sample_excl_list_filename()) << "\".\n" ;
		m_ui_context.logger() << "\n" ;

		m_ui_context.logger() << std::setw(30) << "Input GEN file(s):\n" ;
		m_ui_context.logger()<< m_snp_data_source->get_summary( "", 40 ) ;
		if( m_mangled_options.gen_filename_mapper().input_files().size() > 1 ) {
			m_ui_context.logger() << "\n" ;
		}

		m_ui_context.logger() << std::setw(30) << "Output GEN file(s):" ;
		if( m_mangled_options.gen_filename_mapper().output_filenames().empty() ) {
			m_ui_context.logger() << "  (n/a)\n" ;
		}
		else {
			for( std::size_t i = 0; i < m_mangled_options.gen_filename_mapper().output_filenames().size(); ++i ) {
				if( i > 0 ) {
					m_ui_context.logger() << "\n" << std::string( 30, ' ' ) ;
				}
				m_ui_context.logger() << "  \"" << m_mangled_options.gen_filename_mapper().output_filenames()[i] << "\"" ;
			}
			m_ui_context.logger() << "\n" ;
		}

		m_ui_context.logger() << std::setw(30) << "Output SNP position file(s):" ;
		if( m_mangled_options.snp_excl_list_filename_mapper().output_filenames().empty() ) {
			m_ui_context.logger() << "  (n/a)\n" ;
		}
		else {
			for( std::size_t i = 0; i < m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size(); ++i ) {
				if( i > 0 ) {
					m_ui_context.logger() << "\n" << std::string( 30, ' ' ) ;
				}
				m_ui_context.logger() << "  \"" << m_mangled_options.snp_excl_list_filename_mapper().output_filenames()[i] << "\"" ;
			}
			m_ui_context.logger() << "\n" ;
		}

		m_ui_context.logger() << std::resetiosflags( std::ios::floatfield ) << std::setprecision( 10 ) ;
		{
			if( m_sample_filter.get() ) {
				m_ui_context.logger() << std::setw(30) << "Sample filter:"
					<< "  " << *m_sample_filter << ".\n" ;
			}
		}
		{
			genfile::CommonSNPFilter* snp_filter = get_snp_filter() ;
			if( snp_filter ) {
				m_ui_context.logger() << std::setw(30) << "SNP filter:"
					<< "  " << *snp_filter << ".\n" ;
				m_ui_context.logger() << "\n" ;
			}
		}

		m_ui_context.logger() << std::setw(30) << "# of samples in input files:"
			<< "  " << m_snp_data_source->get_base_source().number_of_samples() << ".\n" ;
		m_ui_context.logger() << std::setw(30) << "# of samples after filtering:"
			<< "  " << m_snp_data_source->number_of_samples()
			<< " (" << m_indices_of_filtered_out_samples.size()
			<< " filtered out).\n" ;

		m_ui_context.logger() << "\n" << std::string( 72, '=' ) << "\n\n" ;

		if( !m_errors.empty() ) {
			for( std::size_t i = 0; i < m_errors.size(); ++i ) {
				m_ui_context.logger() << "!! ERROR: " << m_errors[i] << "\n\n" ;
			}
			m_ui_context.logger() << "!! Please correct the above errors and re-run qctool.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}

		if( !m_warnings.empty() ) {
			for( std::size_t i = 0; i < m_warnings.size(); ++i ) {
				m_ui_context.logger() << "!! WARNING: " << m_warnings[i] << "\n\n" ;
			}
			if( m_ignore_warnings ) {
				m_ui_context.logger() << "!! Warnings were encountered, but proceeding anyway as -force was supplied.\n" ;
				m_ui_context.logger() << "\n" << std::string( 72, '=' ) << "\n\n" ;
			}
			else {
				m_ui_context.logger() << "!! Warnings were encountered.  To proceed anyway, please run again with the -force option.\n" ;
				throw appcontext::HaltProgramWithReturnCode( -1 ) ;
			}
		}

	}

	std::string format_filename( std::string const& filename ) {
		if( filename == "" ) {
			return "(n/a)" ;
		}
		else {
			return filename ;
		}
	}

	void write_postamble() {
		m_ui_context.logger() << std::string( 72, '=' ) << "\n\n" ;
#if 0
		if( m_backup_creator.backed_up_files().size() > 0 ) {
			m_ui_context.logger() << std::setw(36) << "I took backups of the following files before overwriting:\n" ;
			std::size_t max_length = 0u ;
			for(
				std::map< std::string, std::string >::const_iterator i = m_backup_creator.backed_up_files().begin() ;
				i != m_backup_creator.backed_up_files().end() ;
				++i
			) {
				max_length = std::max( max_length, i->first.size() ) ;
			}

			for(
				std::map< std::string, std::string >::const_iterator i = m_backup_creator.backed_up_files().begin() ;
				i != m_backup_creator.backed_up_files().end() ;
				++i
			) {
				m_ui_context.logger() << "  " << std::setw( max_length + 2 ) << std::left << ("\"" + i->first + "\"") << " to \"" << i->second << "\"\n" ;
			}

			m_ui_context.logger() << "\n" ;
			m_ui_context.logger() << std::string( 72, '=' ) << "\n\n" ;
		}
#endif
		m_ui_context.logger() << std::setw(36) << "Number of SNPs:\n"
			<< std::setw(36) << " -- in input file(s):" ;
		if( m_snp_data_source->total_number_of_snps() ) {
			m_ui_context.logger() << "  " << *m_snp_data_source->total_number_of_snps() << ".\n" ;
		}
		else {
			m_ui_context.logger() << "  " << "(not computed).\n" ;
		}

		m_ui_context.logger() << std::setw(36) << " -- in output file(s):"
			<< "  " << m_fltrd_in_snp_data_sink->number_of_snps_written() << "\n" ;

		m_ui_context.logger() << "\n" ;

		m_ui_context.logger() << std::setw(36) << "Number of samples in input file(s):"
			<< "  " << m_snp_data_source->get_base_source().number_of_samples() << ".\n" ;
		if( m_sample_filter->number_of_clauses() > 0 ) {
			for( std::size_t i = 0 ; i < m_sample_filter->number_of_clauses(); ++ i ) {
				m_ui_context.logger() << std::setw(36) << ("...which failed \"" + genfile::string_utils::to_string( m_sample_filter->clause( i )) + "\":")
					<< "  " << ( m_sample_filter_diagnostic_matrix.rows() - m_sample_filter_diagnostic_matrix.col(i).sum() ) << ".\n" ;
			}
			m_ui_context.logger() << std::setw(36) << "(total failures:" << "  " << m_indices_of_filtered_out_samples.size() << ").\n" ;
		}

		m_ui_context.logger() << "\n" ;

		if( m_mangled_options.output_sample_excl_list_filename() != "" ) {
			m_ui_context.logger() << std::setw(36) << "Output sample exclusion list:"
				<< " " << m_mangled_options.output_sample_excl_list_filename() << " ("
				<< m_indices_of_filtered_out_samples.size() << " samples).\n" ;
		}

		if( m_mangled_options.gen_filename_mapper().output_filenames().size() > 0 ) {
			m_ui_context.logger() << std::setw(36) << "Output GEN files:" ;
			for( std::size_t i = 0; i < m_mangled_options.gen_filename_mapper().output_filenames().size(); ++i ) {
				if( i > 0 ) {
					m_ui_context.logger() << std::string( 36, ' ' ) ;
				}
				if( m_fltrd_in_snp_data_sink.get() ) {
					m_ui_context.logger() << "  (" << std::setw(6) << m_fltrd_in_snp_data_sink->sink(i).number_of_snps_written() << " snps)  " ;
				}
				m_ui_context.logger() << "\"" << m_mangled_options.gen_filename_mapper().output_filenames()[i] << "\"\n" ;
			}
			m_ui_context.logger() << std::string( 36, ' ' ) << "  (total " << m_fltrd_in_snp_data_sink->number_of_snps_written() << " snps).\n" ;
		}

		if( m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size() > 0 ) {
			m_ui_context.logger() << std::setw(36) << "Output SNP exclusion list(s):" ;
			for( std::size_t i = 0; i < m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size(); ++i ) {
				if( i > 0 ) {
					m_ui_context.logger() << std::string( 36, ' ' ) ;
				}
				if( m_fltrd_out_snp_data_sink.get() ) {
					m_ui_context.logger() << "  (" << std::setw(6) << m_fltrd_out_snp_data_sink->sink(i).number_of_snps_written() << " snps)  " ;
				}
				m_ui_context.logger() << "  \"" << m_mangled_options.snp_excl_list_filename_mapper().output_filenames()[i] << "\"\n" ;
			}
			m_ui_context.logger() << std::string( 36, ' ' ) << "  (total " << m_fltrd_out_snp_data_sink->number_of_snps_written() << " snps).\n" ;
		}

		if( m_mangled_options.output_sample_filename() != "" ) {
			m_ui_context.logger() << std::setw(36) << "Output SAMPLE files:"
				<< "  \"" << m_mangled_options.output_sample_filename() << "\""
				<< "  (" << m_samples->get_number_of_individuals() << " samples)\n" ;
		}
		if( m_options.check( "-osample" )) {
			m_ui_context.logger() << std::setw(36) << "Sample statistic output file:"
				<< "  \"" << m_options.get< std::string >( "-osample" ) << "\".\n" ;
		}

		if( m_options.check( "-log" )) {
			m_ui_context.logger()[ "screen" ] << std::setw( 36 ) << "\nMore details are in the log file:"
				<< "  \"" << m_mangled_options.log_filename() << "\".\n" ;
		}
		m_ui_context.logger() << std::string( 72, '=' ) << "\n\n" ;
	}

	genfile::CohortIndividualSource const& get_cohort_individual_source() const { assert( m_samples.get() ) ; return *m_samples ; }

private:
	appcontext::OptionProcessor const& m_options ;
	QCToolOptionMangler const m_mangled_options ;
	appcontext::UIContext& m_ui_context ;

	typedef std::map< genfile::VariantIdentifyingData, genfile::VariantIdentifyingData, genfile::VariantIdentifyingData::CompareFields > SNPDictionary ;
	std::auto_ptr< SNPDictionary > m_snp_dictionary ;

	typedef genfile::StrandAligningSNPDataSource::StrandFlipSpec StrandFlipSpec ;
	typedef std::map< genfile::VariantIdentifyingData, StrandFlipSpec, genfile::VariantIdentifyingData::CompareFields > StrandSpec ;
	typedef std::vector< StrandSpec > StrandSpecs ;
	std::auto_ptr< StrandSpecs > m_strand_specs ;
	genfile::CohortIndividualSource::UniquePtr m_samples ; // this must go before the snp_data_sinks.
	std::auto_ptr< genfile::SNPDataSource > m_snp_data_source ;
	std::auto_ptr< genfile::SNPDataSinkChain > m_fltrd_in_snp_data_sink ;
	std::auto_ptr< genfile::SNPDataSinkChain > m_fltrd_out_snp_data_sink ;

	bool m_ignore_warnings ;

	genfile::SampleFilterConjunction::UniquePtr m_sample_filter ;
	mutable genfile::CommonSNPFilter::UniquePtr m_snp_filter ;
	Eigen::MatrixXd m_sample_filter_diagnostic_matrix ;

	std::vector< std::size_t > m_indices_of_filtered_out_samples ;

	std::vector< std::string > m_warnings ;
	std::vector< std::string > m_errors ;

private:

	void setup() {
		m_sample_filter = get_sample_filter() ;
		process_other_options() ;

		if( m_options.check_if_option_has_value( "-map-id-data" )) {
			m_snp_dictionary = load_snp_dictionary( m_options.get_value< std::string >( "-map-id-data" ) ) ;
		}

		if( m_options.check_if_option_has_value( "-strand" )) {
			m_strand_specs = get_strand_specs(
				m_options.get_values< std::string >( "-strand" ),
				m_options.get< std::string >( "-compare-variants-by" )
			) ;
		}

		open_data_sources(
			m_mangled_options.gen_filenames(),
			m_options.check( "-s" ) ? m_options.get_values< std::string >( "-s" ) : boost::optional< std::vector< std::string > >(),
			m_options.check( "-merge-in" ) ? m_options.get_values< std::string >( "-merge-in" ) : boost::optional< std::vector< std::string > >(),
			&m_snp_data_source,
			&m_samples
		) ;

		if( m_options.check( "-sample-data" )) {
			m_samples = open_sample_data( m_samples, m_options.get_values< std::string > ( "-sample-data" )) ;
		}

		{
			if( m_options.check( "-condition-on" )) {
				m_samples = condition_on(
					m_samples,
					*m_snp_data_source,
					genfile::string_utils::join( m_options.get_values< std::string >( "-condition-on" ), "," )
				) ;
				m_snp_data_source->reset_to_start() ;
			}
		}
		m_indices_of_filtered_out_samples = compute_excluded_samples( m_samples ) ;

		if( m_indices_of_filtered_out_samples.size() > 0 ) {
			m_snp_data_source.reset(
				new genfile::SampleFilteringSNPDataSource(
					m_snp_data_source,
					std::set< std::size_t >( m_indices_of_filtered_out_samples.begin(), m_indices_of_filtered_out_samples.end() )
				)
			) ;

			if( m_samples.get() ) {
				m_samples.reset(
					genfile::SampleFilteringCohortIndividualSource::create(
						m_samples,
						std::set< std::size_t >( m_indices_of_filtered_out_samples.begin(), m_indices_of_filtered_out_samples.end() )
					).release()
				) ;
			}
		}

		if( m_options.check_if_option_was_supplied( "-quantile-normalise" )) {
			assert( m_samples.get() ) ;
			m_samples = quantile_normalise_columns(
				m_samples,
				m_options.get_values< std::string >( "-quantile-normalise" )
			) ;
		}

		if( m_options.check( "-reorder" )) {
			std::string const& order_type = m_options.get< std::string > ( "-reorder" ) ;
			using boost::counting_iterator ;
			// generate a random order
			// If order[i] = j then this indicates
			// that sample i in result dataset is sample j in the original dataset.
			std::vector< std::size_t > order(
				counting_iterator< std::size_t >( 0 ),
				counting_iterator< std::size_t >( m_samples->size() )
			) ;
			std::vector< std::string > ids ;
			{
				ids.resize( m_samples->size() ) ;
				for( std::size_t i = 0; i < m_samples->size(); ++i ) {
					ids[i] = m_samples->get_entry( i, "ID_1" ).as< std::string >() ;
				}
			}
			if( boost::filesystem::exists( order_type ) ) {
				if( !m_options.check( "-s" ) && !m_snp_data_source->has_sample_ids() ) {
					throw genfile::BadArgumentError(
						"QCToolCmdLineContext::setup()",
						"-reorder " + order_type,
						"Cannot reorder samples because the input files do not contain sample IDs.  (Do you need to specify -s?)"
					) ;
				}

				// order type is a file specifying the order.
				std::auto_ptr< std::istream > file = genfile::open_text_file_for_input( order_type ) ;
				std::vector< std::string > id_order ;
				std::copy( std::istream_iterator< std::string >( *file ), std::istream_iterator< std::string >(), std::back_inserter< std::vector< std::string > >( id_order )) ;
				if( id_order.size() != ids.size() ) {
					throw genfile::BadArgumentError(
						"QCToolCmdLineContext::setup()",
						"-reorder " + order_type,
						( boost::format( "File has the wrong number of entries (%d, expected %d)" ) % id_order.size() % ids.size() ).str()
					) ;
				}
				for( std::size_t i = 0; i < ids.size(); ++i ) {
					std::vector< std::string >::const_iterator where = std::find( ids.begin(), ids.end(), id_order[i] ) ;
					if( where == ids.end() ) {
						throw genfile::BadArgumentError(
							"QCToolCmdLineContext::setup()",
							"-reorder " + order_type,
							"File contains id (\"" + id_order[i] + "\") that is not found in the dataset."
						) ;
					}
					order[i] = ( where - ids.begin() ) ;
				}
			} else if( order_type == "backwards" ) {
				std::reverse( order.begin(), order.end() ) ;
				std::reverse( ids.begin(), ids.end() ) ;
			} else if( order_type == "randomly" ) {
				// Give the samples a random order, and give them random IDs too.
				{
					unsigned seed = 0 ;
					std::ifstream urandom( "/dev/urandom", std::ios::in | std::ios::binary ) ;
					urandom.read( reinterpret_cast< char* >( &seed ), sizeof( unsigned ) ) ;
					urandom.close() ;
					m_ui_context.logger() << boost::format( "Set random seed to %x.\n" ) % seed ;
					std::srand( seed ) ;
				}
				std::random_shuffle( order.begin(), order.end() ) ;
#if 0
				// This code was to assign UUIDs to each sample
				// Probably easier to do that externally.
				ids = std::vector< std::string >( m_samples->get_number_of_individuals(), "" ) ;
				boost::uuids::random_generator uuidGenerator;
				for( std::size_t i = 0; i < ids.size(); ++i ) {
					ids[i] = boost::uuids::to_string( uuidGenerator() ) ;
				}
				std::sort( ids.begin(), ids.end() ) ;
#endif
			} else {
				throw genfile::BadArgumentError(
					"QCToolCmdLineContext::setup()",
					"-reorder " + order_type,
					"Argument \"" + order_type + "\" should be \"backwards\", \"randomly\", or a file containing ids in the desired order."
				) ;
			}
			m_snp_data_source.reset(
				new genfile::ReorderingSNPDataSource( m_snp_data_source, order )
			) ;
			m_samples.reset(
				new genfile::ReorderingCohortIndividualSource( m_samples, order )
			) ;
		}

		check_for_errors_and_warnings() ;

		write_preamble() ;
		open_snp_data_sinks() ;
	}

	std::auto_ptr< StrandSpecs > get_strand_specs(
		std::vector< std::string > const& filenames,
		genfile::VariantIdentifyingData::CompareFields const& comparator
	) const {
		std::auto_ptr< StrandSpecs > result( new StrandSpecs( filenames.size(), StrandSpec( comparator )) ) ;

		appcontext::UIContext::ProgressContext progress_context = m_ui_context.get_progress_context( "Loading strand files" ) ;
		progress_context.notify_progress( 0, filenames.size() ) ;

		std::vector< std::pair< statfile::BuiltInTypeStatSource::OptionalCount, int > > summaries( filenames.size(), std::make_pair( 0, 0 ) ) ; // first = row count, second = failure count.

		for( std::size_t i = 0; i < filenames.size(); ++i ) {
			statfile::BuiltInTypeStatSource::UniquePtr source( statfile::BuiltInTypeStatSource::open( filenames[i] )) ;
			if( source->number_of_columns() < 6 ) {
				throw genfile::MalformedInputError( filenames[i], 1 ) ;
			}

			std::string platform_column_names ;
			if( source->has_column( "Affy SNP ID" )) {
				platform_column_names = "Chromosome|Physical Position|Probe Set ID|dbSNP RS ID|Allele A|Allele B|Strand" ;
			}
			else if( source->has_column( "IlmnID" )) {
				platform_column_names = "Chromosome|position|IlmnID|rsid|alleleA|alleleB|strand" ;
			}
			else if(
				source->has_column( "chromosome" ) && source->index_of_column( "chromosome" ) == 0
				&& source->has_column( "position" ) && source->index_of_column( "position" ) == 1
				&& source->has_column( "SNPID" ) && source->index_of_column( "SNPID" ) == 2
				&& source->has_column( "rsid" ) && source->index_of_column( "rsid" ) == 3
			) {
				platform_column_names = "chromosome|position|SNPID|rsid|alleleA|alleleB|strand" ;
			}
			else if(
				source->has_column( "chromosome" ) && source->index_of_column( "chromosome" ) == 2
				&& source->has_column( "position" ) && source->index_of_column( "position" ) == 3
				&& source->has_column( "SNPID" ) && source->index_of_column( "SNPID" ) == 0
				&& source->has_column( "rsid" ) && source->index_of_column( "rsid" ) == 1
			) {
				platform_column_names = "chromosome|position|SNPID|rsid|alleleA|alleleB|strand" ;
			}
			else {
				throw genfile::MalformedInputError( source->get_source_spec(), 0 ) ;
			}

			boost::optional< std::size_t > flipColumnIndex ;
			if( m_options.check( "-flip-to-match-allele" )) {
				std::string const columnName = m_options.get< std::string >( "-flip-to-match-allele" ) ;
				if( !source->has_column( columnName ) ) {
					throw genfile::MalformedInputError(
						source->get_source_spec(),
						"Expected column (\"" + columnName + "\") not found in source.", 0
					) ;
				}
				flipColumnIndex = source->index_of_column( columnName ) ;
			}

			summaries[i].first = source->number_of_rows() ; // row count.

			while( *source ) {
				genfile::VariantIdentifyingData snp ;
				std::string SNPID, rsid, allele1, allele2, strand ;
				genfile::GenomePosition position ;
				try {
					if(
						statfile::read_values(
							*source,
							platform_column_names,
							boost::tie(
								position.chromosome(), position.position(),
								SNPID, rsid, allele1, allele2, strand
							)
						)
					) {
						snp = genfile::VariantIdentifyingData(
							SNPID, rsid, position, allele1, allele2
						) ;

						// Support illumina "FWD" and "REV" strands.
						if( strand == "REV" || strand == "-" ) {
							strand = genfile::StrandAligningSNPDataSource::eReverseStrand ;
						}
						else if( strand == "FWD" || strand == "+" ) {
							strand = genfile::StrandAligningSNPDataSource::eForwardStrand ;
						}
						else if( strand == "?" || strand == "NA" ) {
							// std::cerr << "strand = ?: SNP is " << snp << ".\n" ;
							strand = genfile::StrandAligningSNPDataSource::eUnknownStrand ;
						}
						else {
							m_ui_context.logger() << "!! Error: found unrecognised strand \"" << strand << "\" in strand alignment file.\n" ;
							throw genfile::MalformedInputError( source->get_source_spec(), source->number_of_rows_read() + 1 ) ;
						}

						char flip = genfile::StrandAligningSNPDataSource::eNoFlip ;

						if( flipColumnIndex ) {
							flip = genfile::StrandAligningSNPDataSource::eUnknownFlip ;
							(*source) >> statfile::ignore( flipColumnIndex.get() - source->current_column() ) ;
							std::string referenceAllele ;
							(*source) >> referenceAllele ;
							if( strand[0] != genfile::StrandAligningSNPDataSource::eUnknownStrand ) {
								if( referenceAllele == genfile::StrandAligningSNPDataSource::apply_strand( snp.get_allele(0), strand[0] )) {
									flip = genfile::StrandAligningSNPDataSource::eNoFlip ;
								} else if( referenceAllele == genfile::StrandAligningSNPDataSource::apply_strand( snp.get_allele(1), strand[0] )) {
									flip = genfile::StrandAligningSNPDataSource::eFlip ;
								} else {
									throw genfile::MalformedInputError(
										source->get_source_spec(),
										"For SNP ( " + genfile::string_utils::to_string( snp ) + ") with strand " + strand
											+ ": neither alleleA (\"" + snp.get_allele(0) + "\") or alleleB (\"" + snp.get_allele(1)
												+ "\") matches the reference allele (\"" + referenceAllele + "\").",
										source->number_of_rows_read() + 1
									) ;
								}
							}
						}

						StrandSpec::iterator where = result->at(i).find( snp ) ;
						if( where != result->at(i).end() ) {
							throw genfile::DuplicateKeyError( source->get_source_spec(), genfile::string_utils::to_string( snp ) ) ;
						}
						result->at(i)[ snp ] = StrandFlipSpec( strand[0], flip ) ;
						(*source) >> statfile::ignore_all() ;
					}
				}
				catch( genfile::InputError ) {
					++summaries[i].second ; // failure count
					(*source) >> statfile::ignore_all() ;
					//throw ;
				}
			}
			assert( (!source->number_of_rows()) || (*source->number_of_rows() == source->number_of_rows_read() ) ) ;

			progress_context.notify_progress( i+1, filenames.size() ) ;
		}

		m_ui_context.logger() << "\nStrand file summary:\n" ;
		for( std::size_t i = 0; i < summaries.size(); ++i ) {
			if( summaries[i].first ) {
				m_ui_context.logger() << std::setw( 30) << ( "  \"" + filenames[i] + "\"" ) << ": " << std::setw( 7 )
					<< *summaries[i].first << " rows" ;
			} else {
				m_ui_context.logger() << std::setw( 30) << ( "  \"" + filenames[i] + "\"" ) << ": unknown number of rows" ;
			}
			if( summaries[i].second > 0 ) {
				m_ui_context.logger() << " of which " << summaries[i].second << " may be malformed, and will be ignored" ;
			}
			m_ui_context.logger() << ".\n" ;
		}

		return result ;
	}

public:
	genfile::SNPDataSource::UniquePtr open_snp_data_sources(
		std::vector< genfile::wildcard::FilenameMatch > const& filenames
	) const {
		return open_snp_data_sources(
			std::vector< std::vector< genfile::wildcard::FilenameMatch > >( 1, filenames ),
			false
		) ;
	}

public:
	void open_data_sources(
		std::vector< std::vector< genfile::wildcard::FilenameMatch > > const& gen_filenames,
		boost::optional< std::vector< std::string > > const& sample_filenames,
		boost::optional< std::vector< std::string > > const& merge_in_sources,
		genfile::SNPDataSource::UniquePtr* result_snp_data_source,
		genfile::CohortIndividualSource::UniquePtr* result_samples
	) {
		genfile::CohortIndividualSource::UniquePtr samples ;
		genfile::SNPDataSource::UniquePtr snp_data_source ;

		// Load sample file(s) if they are provided
		if( sample_filenames > 0 ) {
			// This is currently supposed to be checked elsewhere
			assert( sample_filenames->size() == gen_filenames.size() ) ;
			samples = open_samples( *sample_filenames ) ;
		}
		snp_data_source = open_snp_data_sources( gen_filenames, m_options.check( "-flip-to-match-cohort1" ) ) ;

		// From this point on we need a sample source.
		// Either one was provided (in which case we check its IDs against the genotype file, if there are any)
		// or none was provided (in which case we create an in-memory one.)
		if( samples.get() && snp_data_source->has_sample_ids() ) {
			check_sample_ids( *samples, *snp_data_source ) ;
		} else if( !samples.get() ) {
			// Create a dummy sample file if none was provided
			samples = create_dummy_samples( *snp_data_source ) ;
		}

		// Check Sample and snp data sources match.
		if( samples->get_number_of_individuals() != snp_data_source->number_of_samples() ) {
			throw genfile::MismatchError(
				"QCToolCmdLineContext::open_data_sources()",
				samples->get_source_spec(),
				"number of samples = " + genfile::string_utils::to_string( samples->get_number_of_individuals() ),
				"expected number of samples = " + genfile::string_utils::to_string( snp_data_source->number_of_samples() )
			) ;
		}

		if( merge_in_sources ) {
			if( merge_in_sources->size() != 2 ) {
				throw genfile::BadArgumentError(
					"open_data_sources()",
					"merge_in_source",
					"Expected two entries, a SNP data source and a sample file"
				) ;
			}
			snp_data_source = open_merged_data_sources(
				snp_data_source, *samples,
				(*merge_in_sources)[0], (*merge_in_sources)[1]
			) ;
		}

		if( m_options.check( "-infer-ploidy-from" )) {
			std::string const sex_column = m_options.get< std::string >( "-infer-ploidy-from" ) ;
			std::vector< char > sexes( samples->size(), '.' ) ;
			for( std::size_t i = 0; i < samples->size(); ++i ) {
				genfile::VariantEntry const& entry = samples->get_entry( i, sex_column ) ;
				if( !entry.is_missing() ) {
					std::string const sex = genfile::string_utils::to_lower( entry.as< std::string >() ) ;
					if( sex == "1" || sex == "m" || sex == "male" ) {
						sexes[i] = 'm' ;
					} else if( sex == "2" || sex == "f" || sex == "female" ) {
						sexes[i] = 'f' ;
					} else {
						throw genfile::MalformedInputError(
							samples->get_source_spec(),
							"Malformed sex value \"" + sex + "\"",
							i,
							samples->get_column_spec().find_column( sex_column )
						) ;
					}
				}
			}

			snp_data_source.reset(
				genfile::PloidyConvertingSNPDataSource::create(
					snp_data_source,
					boost::bind( &genfile::get_ploidy_from_sex, sexes, _1, _2 )
				).release()
			) ;
		}

		if( m_options.check( "-threshold" )) {
			snp_data_source.reset(
				new genfile::ThreshholdingSNPDataSource( snp_data_source, m_options.get< double >( "-threshold" ))
			) ;
		}

		// Put results in output variables
		*result_samples = samples ;
		*result_snp_data_source = snp_data_source ;
	}

private:

	void check_sample_ids(
		genfile::CohortIndividualSource const& samples,
		genfile::SNPDataSource const& snp_data_source
	) const {
		// check sample identifiers agree
		// if not, abort unless -force is specified
		std::vector< std::string > sample_file_ids( samples.size() ) ;
		std::vector< std::string > genotype_file_ids( snp_data_source.number_of_samples() ) ;
		samples.get_column_values( "ID_1", boost::bind( &impl::set_vector_entry_variant, &sample_file_ids, _1, _2 )) ;
		snp_data_source.get_sample_ids( boost::bind( &impl::set_vector_entry_string, &genotype_file_ids, _1, _2 )) ;
		assert( sample_file_ids.size() == genotype_file_ids.size() ) ;
		if( sample_file_ids != genotype_file_ids ) {
			// find first mismatch for a useful error message
			std::size_t k = 0 ;
			for( ; k < sample_file_ids.size(); ++k ) {
				if( sample_file_ids[k] != genotype_file_ids[k] ) {
					break ;
				}
			}
			m_ui_context.logger() << "!! Error: sample file and genotype file sample IDs do not match!\n" ;
			m_ui_context.logger() << "!! First mismatch is at sample " << (k+1) << ": IDs are \"" << sample_file_ids[k] << "\" and \"" << genotype_file_ids[k] << "\".\n" ;
			if( m_options.check( "-force" ) ) {
				m_ui_context.logger() << "!! ...continuing anyway as -force was specified.\n" ;
			} else {
				throw genfile::MismatchError(
					"QCToolCmdLineContext::check_sample_ids()",
					genfile::string_utils::join( m_options.get_values< std::string >( "-s" ), "," ),
					sample_file_ids[k],
					genotype_file_ids[k]
				) ;
			}
		}
	}

private:
	genfile::SNPDataSource::UniquePtr open_snp_data_sources(
		std::vector< std::vector< genfile::wildcard::FilenameMatch > > const& filenames,
		bool match_alleles_between_cohorts
	) const {
		genfile::SNPDataSourceRack::UniquePtr rack ;
		if( filenames.size() > 1 ) {
			rack.reset(
				new genfile::SNPDataSourceRack(
					genfile::VariantIdentifyingData::CompareFields(
						m_options.get_value< std::string >( "-compare-variants-by" ),
						match_alleles_between_cohorts
					)
				)
			) ;
		}

		// count files.
		std::size_t file_count = 0 ;
		for( std::size_t i = 0; i < filenames.size(); ++i ) {
			file_count += filenames[i].size() ;
		}

		appcontext::UIContext::ProgressContext progress_context = m_ui_context.get_progress_context( "Opening genotype files" ) ;
		progress_context.notify_progress( 0, file_count ) ;

		std::size_t progress_count = 0 ;

		genfile::SNPDataSource::UniquePtr source ;

		std::vector< genfile::VariantIdentifyingData > cohort1_snps ;

		for( std::size_t i = 0; i < filenames.size(); ++i ) {
			genfile::SNPDataSourceChain::UniquePtr chain( new genfile::SNPDataSourceChain() ) ;

			for( std::size_t j = 0; j < filenames[i].size(); ++j ) {
				source = open_snp_data_source(
					filenames[i][j].filename(),
					filenames[i][j].match()
				) ;

				// Add the source to the chain.
				chain->add_source( genfile::SNPDataSource::UniquePtr( source ) ) ;

				progress_context.notify_progress( ++progress_count, file_count ) ;
			}
			source.reset( chain.release() ) ;

			// If we have strand alignment information, implement it now
			if( m_strand_specs.get() ) {
				if( m_strand_specs->size() != filenames.size() ) {
					throw genfile::BadArgumentError(
						"QCToolCmdLineContext::open_snp_data_sources()",
						"-strand",
						"Number of strand files (" + genfile::string_utils::to_string( m_strand_specs->size() ) + ") should match number of genotype files (" + genfile::string_utils::to_string( filenames.size() ) + ")."
					) ;
				}
				source.reset(
					genfile::StrandAligningSNPDataSource::create(
						source,
						m_strand_specs->at(i)
					)
					.release()
				) ;
			}

			if( rack.get() ) {
				rack->add_source( source ) ;
			} else {
				break ;
			}
		}

		if( rack.get() ) {
			source.reset( rack.release() ) ;
		}

		return source ;
	}

	// Find mapping of needle samples in haystack samples
	std::vector< std::size_t > get_sample_id_mapping(
		genfile::CohortIndividualSource const& needle_samples,
		std::string const& needle_column,
		genfile::CohortIndividualSource const& haystack_samples,
		std::string const& haystack_column
	) {
		std::size_t const eNotIncluded = genfile::ReorderingSNPDataSource::eNotIncluded ;
		std::vector< std::size_t > result( needle_samples.size() ) ;

		typedef std::vector< genfile::VariantEntry > IdList ;
		IdList haystack_ids ;
		IdList needle_ids ;
		void(IdList::*push_back)(genfile::VariantEntry const&) = &IdList::push_back;
		needle_samples.get_column_values( needle_column, boost::bind( push_back, &needle_ids, _2 ) ) ;
		haystack_samples.get_column_values( haystack_column, boost::bind( push_back, &haystack_ids, _2 ) ) ;

		assert( needle_ids.size() == needle_samples.size() ) ;
		assert( haystack_ids.size() == haystack_samples.size() ) ;

		// typedef boost::unordered_map< VariantEntry, std::size_t > HaystackMap ;
		typedef std::map< genfile::VariantEntry, std::size_t > IdMap ;
		IdMap haystack_map ;
		for( std::size_t i = 0; i < haystack_ids.size(); ++i ) {
			if( !haystack_map.insert( std::make_pair( haystack_ids[i], i )).second ) {
				throw genfile::BadArgumentError(
					"genfile::impl::SampleMapping::setup()",
					"haystack_column=\"" + haystack_column + "\"",
					"Haystack sample " + haystack_column + "=\"" + haystack_ids[i].as< std::string >() + "\" occurs more than once."
				) ;
			} ;
		}
		for( std::size_t i = 0; i < needle_ids.size(); ++i ) {
			IdMap::const_iterator where = haystack_map.find( needle_ids[i] ) ;
			result[i] = ( where == haystack_map.end() ) ? eNotIncluded : where->second ;
		}
		return result ;
	}

	genfile::SNPDataSource::UniquePtr open_merged_data_sources(
		genfile::SNPDataSource::UniquePtr snp_data_source,
		genfile::CohortIndividualSource const& samples,
		std::string const& merge_in_gen_file,
		std::string const& merge_in_sample_file
	) {
		genfile::MergingSNPDataSource::UniquePtr merged_source = genfile::MergingSNPDataSource::create(
			m_options.get< std::string >( "-merge-strategy" ),
			genfile::VariantIdentifyingData::CompareFields( m_options.get_value< std::string >( "-compare-variants-by" ) )
		) ;

		merged_source->add_source( snp_data_source ) ;
		std::string id_prefix = "" ;
		if( m_options.check( "-merge-prefix" )) {
			id_prefix = m_options.get< std::string >( "-merge-prefix" ) ;
		}

		{
			genfile::SNPDataSource::UniquePtr merge_in_source ;
			{
				genfile::SNPDataSourceChain::UniquePtr source( new genfile::SNPDataSourceChain ) ;
				std::vector< genfile::wildcard::FilenameMatch > filenames = genfile::wildcard::find_files_by_chromosome( merge_in_gen_file ) ;
				for( std::size_t i = 0; i < filenames.size(); ++i ) {
					source->add_source( open_snp_data_source( filenames[i].filename(), filenames[i].match())) ;
				}
				merge_in_source.reset( source.release() ) ;
			}
			genfile::CohortIndividualSource::UniquePtr merge_in_samples(
				new genfile::CategoricalCohortIndividualSource(
					merge_in_sample_file,
					m_options.get_value< std::string >( "-missing-code" )
				)
			) ;

			std::vector< std::string > id_columns = genfile::string_utils::split_and_strip( m_options.get_value< std::string >( "-match-sample-ids" ), "~" ) ;
			if( id_columns.size() != 2 ) {
				throw genfile::BadArgumentError(
					"QCToolCmdLineContext::open_merged_data_sources()",
					"-match-sample-ids=\"" + m_options.get_value< std::string >( "-match-sample-ids" ) + "\"",
					"Value should be two strings separated by ~."
				) ;
			}

			merge_in_source = genfile::ReorderingSNPDataSource::create(
				merge_in_source,
				get_sample_id_mapping( samples, id_columns[0], *merge_in_samples, id_columns[1] )
			) ;

			// Make the merged-in source respect the filter.
			genfile::CommonSNPFilter* snp_filter = get_snp_filter() ;
			if( snp_filter ) {
				merge_in_source = genfile::VariantIdentifyingDataFilteringSNPDataSource::create(
					merge_in_source,
					*snp_filter
				) ;
			}

			merged_source->add_source(
				merge_in_source,
				id_prefix
			) ;
		}

		return genfile::SNPDataSource::UniquePtr( merged_source.release() ) ;
	}

	void write_excluded_SNP( genfile::VariantIdentifyingData const& snp ) const {
		if( m_fltrd_out_snp_data_sink.get() ) {
			Eigen::MatrixXd matrix( m_samples->get_number_of_individuals(), 3 ) ;
			matrix.setZero() ;
			m_fltrd_out_snp_data_sink->write_snp(
				m_samples->get_number_of_individuals(),
				snp,
				genfile::GenotypeGetter< Eigen::MatrixXd >( matrix, 0 ),
				genfile::GenotypeGetter< Eigen::MatrixXd >( matrix, 1 ),
				genfile::GenotypeGetter< Eigen::MatrixXd >( matrix, 2 )
			) ;
		}
	}

	genfile::SNPDataSource::UniquePtr
	open_snp_data_source( std::string const& filename, std::string chromosome_indicator ) const {
		genfile::Chromosome chromosome_hint ;
		if( chromosome_indicator == "" && m_options.check_if_option_was_supplied( "-assume-chromosome" )) {
			chromosome_indicator = m_options.get_value< std::string >( "-assume-chromosome" ) ;
		}

		if( chromosome_indicator != "" ) {
			chromosome_hint = chromosome_indicator ;
		}
		genfile::SNPDataSource::UniquePtr source ;

		std::pair< std::string, std::string > uf = genfile::uniformise( filename ) ;

		{
			boost::optional< genfile::vcf::MetadataParser::Metadata > metadata ;
			if( m_options.check( "-metadata" )) {
				metadata = genfile::vcf::StrictMetadataParser(
					m_options.get_value< std::string >( "-metadata" )
				).get_metadata() ;
			}
			source = genfile::SNPDataSource::create(
				filename,
				chromosome_hint,
				metadata,
				m_options.get< std::string >( "-filetype" )
			) ;
		}

		// vcf-specific options
		if( genfile::VCFFormatSNPDataSource* vcf_source = dynamic_cast< genfile::VCFFormatSNPDataSource* >( source.get() ) ) {
			std::string genotype_field = m_options.get< std::string >( "-vcf-genotype-field" ) ;
			vcf_source->set_field_mapping( ":genotypes:", genotype_field ) ;
			std::string intensity_field = m_options.get< std::string >( "-vcf-intensity-field" ) ;
			vcf_source->set_field_mapping( ":intensities:", intensity_field ) ;
			vcf_source->set_strict_mode( !m_options.check( "-permissive" )) ;
		}

		genfile::CommonSNPFilter* snp_filter = get_snp_filter() ;
		// Filter SNPs if necessary
		if( snp_filter ) {
			genfile::VariantIdentifyingDataFilteringSNPDataSource::UniquePtr snp_filtering_source
				= genfile::VariantIdentifyingDataFilteringSNPDataSource::create(
					source,
					*snp_filter
				) ;

			if( m_options.check( "-write-snp-excl-list" )) {
				snp_filtering_source->send_filtered_out_SNPs_to(
					boost::bind(
						&QCToolCmdLineContext::write_excluded_SNP,
						this,
						_1
					)
				) ;
			}

			source.reset(
				snp_filtering_source.release()
			) ;
		}

		// Translate SNP identifying data if necessary
		if( m_snp_dictionary.get() ) {
			source.reset(
				genfile::SNPTranslatingSNPDataSource::create(
					source,
					*m_snp_dictionary
				).release()
			) ;
		}

		return source ;
	}

	genfile::SNPDataSource::UniquePtr
	open_vcf_format_snp_data_source( std::pair< std::string, std::string > const& uf ) const {
		genfile::VCFFormatSNPDataSource::UniquePtr source ;
		genfile::vcf::MetadataParser::Metadata metadata ;
		if( m_options.check( "-metadata" )) {
			metadata = genfile::vcf::StrictMetadataParser(
				m_options.get_value< std::string >( "-metadata" )
			).get_metadata() ;
		}

		if( uf.second == "-" ) {
			std::cerr << "BOOOOOOD\n" ;
			std::auto_ptr< std::istream > str( new std::istream( std::cin.rdbuf() ) ) ;
			source.reset(
				new genfile::VCFFormatSNPDataSource(
					str,
					metadata
				)
			) ;
		} else {
			source.reset(
				new genfile::VCFFormatSNPDataSource( uf.second )
			) ;
		}

		std::string genotype_field = m_options.get< std::string >( "-vcf-genotype-field" ) ;
		source->set_field_mapping( ":genotypes:", genotype_field ) ;
		std::string intensity_field = m_options.get< std::string >( "-vcf-intensity-field" ) ;
		source->set_field_mapping( ":intensities:", intensity_field ) ;
		source->set_strict_mode( !m_options.check( "-permissive" )) ;

		return genfile::SNPDataSource::UniquePtr( source.release() ) ;
	}

	std::auto_ptr< SNPDictionary > load_snp_dictionary( std::string const& filename ) const {
		std::auto_ptr< SNPDictionary > result( new SNPDictionary( genfile::VariantIdentifyingData::CompareFields( m_options.get_value< std::string >( "-compare-variants-by" ) ) )) ;
		statfile::BuiltInTypeStatSource::UniquePtr source(
			statfile::BuiltInTypeStatSource::open(
				genfile::wildcard::find_files_by_chromosome( filename )
			).release()
		) ;

		appcontext::UIContext::ProgressContext progress_context = m_ui_context.get_progress_context( "Opening position translation dictionary" ) ;

		if( source->number_of_columns() != 12 ) {
			throw genfile::MalformedInputError( filename, 1 ) ;
		}

		bool found_duplicates = false ;
		std::string SNPID1, SNPID2, rsid1, rsid2, allele11, allele12, allele21, allele22 ;
		genfile::GenomePosition position1, position2 ;
		while(
			(*source)
				>> SNPID1 >> rsid1
				>> position1.chromosome() >> position1.position()
				>> allele11 >> allele12
				>> SNPID2 >> rsid2
				>> position2.chromosome() >> position2.position()
				>> allele21 >> allele22
		) {
			genfile::VariantIdentifyingData
				data1( SNPID1, rsid1, position1, allele11, allele12 ),
				data2( SNPID2, rsid2, position2, allele21, allele22 ) ;

			std::map< genfile::VariantIdentifyingData, genfile::VariantIdentifyingData >::const_iterator where = result->find( data1 ) ;
			if( where != result->end() ) {
				if( where->second == data2 && !found_duplicates ) {
					m_ui_context.logger() << "qctool::load_snp_dictionary(): found duplicate rows in source "
						<< source->get_source_spec()
						<< ".\nFirst duplicated variant: " << data1 << ".\n" ;
					found_duplicates = true ;
				} else if( where->second != data2 ) {
					m_ui_context.logger() << "qctool::load_snp_dictionary(): found mismatching duplicate rows in source "
						<< source->get_source_spec()
						<< " for variant: " << data1 << ".\n" ;
					throw genfile::DuplicateSNPError( filename, genfile::string_utils::to_string( data1 ) ) ;
				}
			}
			(*result)[ data1 ] = data2 ;
			(*source) >> statfile::end_row() ;
			progress_context( source->number_of_rows_read(), source->number_of_rows() ) ;
		}
		return result ;
	}

	genfile::CommonSNPFilter* get_snp_filter() const {
		if( !m_snp_filter.get() ) {
			m_snp_filter = construct_snp_filter() ;
		}
		return m_snp_filter.get() ;
	}

	genfile::CommonSNPFilter::UniquePtr construct_snp_filter() const {
		genfile::CommonSNPFilter::UniquePtr snp_filter ;

		if( m_options.check_if_option_was_supplied_in_group( "SNP exclusion options" )) {
			snp_filter.reset( new genfile::CommonSNPFilter ) ;

			if( m_options.check_if_option_was_supplied( "-excl-variants" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-excl-variants" ) ;

				BOOST_FOREACH( std::string const& filename, files ) {
					genfile::SNPDataSource::UniquePtr source ;
					try {
						source.reset(
							new statfile::SNPDataSourceAdapter(
								statfile::BuiltInTypeStatSource::open(
									genfile::wildcard::find_files_by_chromosome( filename )
								)
							)
						) ;
					} catch( genfile::InputError const& e ) {
						source.reset(
							genfile::SNPDataSourceChain::create(
								genfile::wildcard::find_files_by_chromosome( filename )
							).release()
						) ;
					}
					snp_filter->exclude_snps(
						source->list_snps(),
						genfile::VariantIdentifyingData::CompareFields( m_options.get_value< std::string >( "-compare-variants-by" ) )
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-variants" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-incl-variants" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					genfile::SNPDataSource::UniquePtr source ;
					try {
						source.reset(
							new statfile::SNPDataSourceAdapter(
								statfile::BuiltInTypeStatSource::open(
									genfile::wildcard::find_files_by_chromosome( filename )
								)
							)
						) ;
					} catch( genfile::InputError const& e ) {
						std::cerr << "Uh-oh...\n" ;
						source.reset(
							genfile::SNPDataSourceChain::create(
								genfile::wildcard::find_files_by_chromosome( filename )
							).release()
						) ;
					}
					snp_filter->include_snps(
						source->list_snps(),
						genfile::VariantIdentifyingData::CompareFields( m_options.get_value< std::string >( "-compare-variants-by" ) )
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-excl-snpids" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-excl-snpids" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->exclude_snps_in_file(
						filename,
						genfile::CommonSNPFilter::SNPIDs
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-snpids" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-incl-snpids" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->include_snps_in_file(
						filename,
						genfile::CommonSNPFilter::SNPIDs
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-excl-rsids" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-excl-rsids" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->exclude_snps_in_file(
						filename,
						genfile::CommonSNPFilter::RSIDs
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-rsids" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-incl-rsids" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->include_snps_in_file(
						filename,
						genfile::CommonSNPFilter::RSIDs
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-excl-positions" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-excl-positions" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->exclude_snps_in_file(
						filename,
						genfile::CommonSNPFilter::Positions
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-positions" )) {
				std::vector< std::string > files = m_options.get_values< std::string > ( "-incl-positions" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					snp_filter->include_snps_in_file(
						filename,
						genfile::CommonSNPFilter::Positions
					) ;
				}
			}
			if( m_options.check_if_option_was_supplied( "-excl-variants-matching" )) {
				std::vector< std::string > specs = m_options.get_values< std::string >( "-excl-variants-matching" ) ;
				BOOST_FOREACH( std::string const& spec, specs ) {
					snp_filter->exclude_snps_matching(
						spec
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-variants-matching" )) {
				std::vector< std::string > specs = m_options.get_values< std::string >( "-incl-variants-matching" ) ;
				BOOST_FOREACH( std::string const& spec, specs ) {
					snp_filter->include_snps_matching(
						spec
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-range" )) {
				std::vector< std::string > specs = m_options.get_values< std::string >( "-incl-range" ) ;
				for ( std::size_t i = 0; i < specs.size(); ++i ) {
					snp_filter->include_snps_in_range(
						genfile::GenomePositionRange::parse( specs[i] )
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-excl-range" )) {
				std::vector< std::string > specs = m_options.get_values< std::string >( "-excl-range" ) ;
				for ( std::size_t i = 0; i < specs.size(); ++i ) {
					snp_filter->exclude_snps_in_range(
						genfile::GenomePositionRange::parse( specs[i] )
					) ;
				}
			}

			if( m_options.check_if_option_was_supplied( "-incl-ranges" )) {
				std::vector< std::string > files = m_options.get_values< std::string >( "-incl-ranges" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					std::auto_ptr< std::istream > in = genfile::open_text_file_for_input( filename ) ;
					std::string range ;
					while( (*in) >> range ) {
						snp_filter->include_snps_in_range(
							genfile::GenomePositionRange::parse( range )
						) ;
					}
				}
			}

			if( m_options.check_if_option_was_supplied( "-excl-ranges" )) {
				std::vector< std::string > files = m_options.get_values< std::string >( "-excl-ranges" ) ;
				BOOST_FOREACH( std::string const& filename, files ) {
					std::auto_ptr< std::istream > in = genfile::open_text_file_for_input( filename ) ;
					std::string range ;
					while( (*in) >> range ) {
						snp_filter->exclude_snps_in_range(
							genfile::GenomePositionRange::parse( range )
						) ;
					}
				}
			}
		}
		return snp_filter ;
	}

	void move_to_next_output_file( std::size_t index ) {
		// Kludge: don't move through outputs until they are actually created.
		// Actually, we only want to move through outputs once during the main
		// processing, but for example id filtering takes place first and has nothing
		// to do with the outputs. On the other hand, filtering should probably
		// happen before the
		if( m_fltrd_in_snp_data_sink.get() ) {
			if( index < m_mangled_options.gen_filename_mapper().input_files().size() ) {
				if( m_mangled_options.gen_filename_mapper().output_filenames().size() > 0 ) {
					if( m_mangled_options.gen_filename_mapper().filename_corresponding_to( index ) != m_fltrd_in_snp_data_sink->index_of_current_sink() ) {
						m_fltrd_in_snp_data_sink->move_to_next_sink() ;
					}
				}

				if( m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size() > 0 ) {
					if( m_mangled_options.snp_excl_list_filename_mapper().filename_corresponding_to( index ) != m_fltrd_out_snp_data_sink->index_of_current_sink() ) {
						m_fltrd_out_snp_data_sink->move_to_next_sink() ;
					}
				}
			}
		}
	}

	void open_snp_data_sinks() {
		open_filtered_in_snp_data_sink() ;
		open_filtered_out_snp_data_sink() ;
	}

	void reset_filtered_in_snp_data_sink() {
		m_fltrd_in_snp_data_sink = std::auto_ptr< genfile::SNPDataSinkChain >( new genfile::SNPDataSinkChain() ) ;
	}

	void open_filtered_in_snp_data_sink() {
		reset_filtered_in_snp_data_sink() ;
		if( m_mangled_options.gen_filename_mapper().output_filenames().size() == 0 ) {
			m_fltrd_in_snp_data_sink->add_sink( genfile::SNPDataSink::UniquePtr( new genfile::TrivialSNPDataSink() )) ;
		}
		else {
			for( std::size_t i = 0; i < m_mangled_options.gen_filename_mapper().output_filenames().size(); ++i ) {
				std::string const& filename = m_mangled_options.gen_filename_mapper().output_filenames()[i] ;
				genfile::SNPDataSink::UniquePtr sink ;
				if( m_options.get< std::string >( "-ofiletype" ) == "sqlite_haplotypes" ) {
					sink.reset(
						new SQLiteHaplotypesSNPDataSink(
							qcdb::DBOutputter::create(
								filename,
								m_options.get< std::string >( "-analysis-name" ),
								m_options.get< std::string >( "-analysis-chunk" ),
								m_options.get_values_as_map(),
								boost::optional< genfile::db::Connection::RowId >(),
								m_options.get< std::string >( "-compare-variants-by" )
							),
							*m_samples
						)
					) ;
				} else if( m_options.get< std::string >( "-ofiletype" ) == "sqlite_genotypes" ) {
					sink.reset(
						new SQLiteGenotypesSNPDataSink(
							qcdb::DBOutputter::create(
								filename,
								m_options.get< std::string >( "-analysis-name" ),
								m_options.get< std::string >( "-analysis-chunk" ),
								m_options.get_values_as_map(),
								boost::optional< genfile::db::Connection::RowId >(),
								m_options.get< std::string >( "-compare-variants-by" )
							),
							*m_samples
						)
					) ;
			 	}
				else {
					sink = genfile::SNPDataSink::create(
						filename,
						genfile::SNPDataSink::Metadata(),
						m_options.get< std::string >( "-ofiletype" )
					) ;
					if( m_options.check_if_option_was_supplied( "-omit-chromosome" )) {
						genfile::GenLikeSNPDataSink* gen_sink = dynamic_cast< genfile::GenLikeSNPDataSink* >( sink.get() ) ;
						if( gen_sink ) {
							gen_sink->omit_chromosome() ;
						}
					}
					if( m_options.check_if_option_was_supplied( "-vcf-output-fields" )) {
						genfile::VCFFormatSNPDataSink* vcf_sink = dynamic_cast< genfile::VCFFormatSNPDataSink* >( sink.get() ) ;
						if( vcf_sink ) {
							std::vector< std::string > const values = m_options.get_values< std::string >( "-vcf-output-fields" ) ;
							vcf_sink->set_output_fields( std::set< std::string >( values.begin(), values.end() ) ) ;
						}
					}
					// bgen-specific options
					{
						genfile::BGenFileSNPDataSink* bgen_sink = dynamic_cast< genfile::BGenFileSNPDataSink* >( sink.get() ) ;
						if( bgen_sink ) {
							bgen_sink->set_number_of_bits( m_options.get< std::size_t >( "-bgen-bits" )) ;
							bgen_sink->set_compression_type( m_options.get< std::string >( "-bgen-compression" ) ) ;
							bgen_sink->set_permitted_input_rounding_error( m_options.get< double >( "-bgen-permitted-input-rounding-error" ) ) ;

							if( m_options.check( "-bgen-free-data" )) {
								bgen_sink->set_free_data( m_options.get< std::string >( "-bgen-free-data" ) ) ;
							}

							bgen_sink->set_write_sample_identifier_block(
								(m_snp_data_source->has_sample_ids() || m_options.check( "-s" ))
								&& ! m_options.check( "-bgen-omit-sample-identifier-block" )
							) ;
						}
					}
					// gen-specific options
					{
						genfile::GenFileSNPDataSink* gen_sink = dynamic_cast< genfile::GenFileSNPDataSink* >( sink.get() ) ;
						if( gen_sink ) {
							if( m_options.check( "-precision" )) {
								gen_sink->set_precision( m_options.get< std::size_t >( "-precision" )) ;
							}
						}
					}
					if( m_options.check( "-sort" )) {
						sink.reset(
							new genfile::SortingBGenFileSNPDataSink(
								filename,
								sink,
								m_options.get< std::string >( "-compare-variants-by" )
							)
						) ;
					}
				}
				m_fltrd_in_snp_data_sink->add_sink( sink ) ;
			}
		}
	}

	void open_filtered_out_snp_data_sink() {
		reset_filtered_out_snp_data_sink() ;
		if( m_options.check_if_option_was_supplied( "-write-snp-excl-list" ) && m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size() > 0 ) {
			for( std::size_t i = 0; i < m_mangled_options.snp_excl_list_filename_mapper().output_filenames().size(); ++i ) {
				std::string const& filename = m_mangled_options.snp_excl_list_filename_mapper().output_filenames()[i] ;
				m_fltrd_out_snp_data_sink->add_sink( genfile::SNPDataSink::UniquePtr( new SNPIDSink( filename ))) ;
			}
		}
		else {
			m_fltrd_out_snp_data_sink->add_sink( genfile::SNPDataSink::UniquePtr( new genfile::TrivialSNPDataSink() )) ;
		}
	}

	void reset_filtered_out_snp_data_sink() {
		m_fltrd_out_snp_data_sink.reset( new genfile::SNPDataSinkChain() ) ;
	}

	genfile::SampleFilterConjunction::UniquePtr get_sample_filter() const {
		genfile::SampleFilterConjunction::UniquePtr filter( new genfile::SampleFilterConjunction() ) ;

		if( m_options.check_if_option_was_supplied( "-sample-missing-rate" ) ) {
			filter->add_clause(
				genfile::SampleFilter::UniquePtr(
					new genfile::VariableInInclusiveRangeSampleFilter( "missing", -1, m_options.get_value< double >( "-sample-missing-rate" ) )
				)
			) ;
		}

		if( m_options.check_if_option_was_supplied( "-heterozygosity" ) ) {
			filter->add_clause(
				genfile::SampleFilter::UniquePtr(
					new genfile::VariableInInclusiveRangeSampleFilter( "heterozygosity", m_options.get_values< double >( "-heterozygosity" )[0], m_options.get_values< double >( "-heterozygosity" )[1] )
				)
			) ;
		}

		genfile::SampleFilterDisjunction::UniquePtr sample_inclusion_filter( new genfile::SampleFilterDisjunction() ) ;

		if( m_options.check_if_option_was_supplied( "-incl-samples" ) ) {
			genfile::SampleFilter::UniquePtr id1_filter = get_sample_id_filter(
				m_options.get_values< std::string >( "-incl-samples" )
			) ;
			sample_inclusion_filter->add_clause( genfile::SampleFilter::UniquePtr( id1_filter.release() ) ) ;
		}

		if( m_options.check_if_option_was_supplied( "-excl-samples" ) ) {
			genfile::SampleFilter::UniquePtr id1_filter = get_sample_id_filter(
				m_options.get_values< std::string >( "-excl-samples" )
			) ;
			filter->add_clause(
				genfile::SampleFilter::UniquePtr(
					new genfile::SampleFilterNegation(
						id1_filter
					)
				)
			) ;
		}

		if( m_options.check_if_option_was_supplied( "-incl-samples-where" ) ) {
			std::vector< std::string > conditions = m_options.get_values< std::string >( "-incl-samples-where" ) ;
			// we OR together different WHERE clauses.
			for( std::size_t i = 0; i < conditions.size(); ++i ) {
				sample_inclusion_filter->add_clause( genfile::SampleFilter::create( conditions[i] )) ;
			}
		}

		if( m_options.check_if_option_was_supplied( "-excl-samples-where" ) ) {
			std::vector< std::string > conditions = m_options.get_values< std::string >( "-excl-samples-where" ) ;
			// we AND together different WHERE clauses.
			genfile::SampleFilterDisjunction::UniquePtr where( new genfile::SampleFilterDisjunction() ) ;
			for( std::size_t i = 0; i < conditions.size(); ++i ) {
				where->add_clause( genfile::SampleFilter::create( conditions[i] )) ;
			}
			filter->add_clause(
				genfile::SampleFilter::UniquePtr(
					new genfile::SampleFilterNegation(
						genfile::SampleFilter::UniquePtr( where.release() )
					)
				)
			) ;
		}
		if( sample_inclusion_filter->number_of_clauses() > 0 ) {
			filter->add_clause( genfile::SampleFilter::UniquePtr( sample_inclusion_filter.release() ) ) ;
		}

		return filter ;
	}

	genfile::SampleFilter::UniquePtr get_sample_id_filter( std::vector< std::string > const& filenames ) const {
		genfile::VariableInSetSampleFilter::UniquePtr result( new genfile::VariableInSetSampleFilter( "ID_1" )) ;
		for( std::size_t i = 0; i < filenames.size(); ++i ) {
			std::string elt ;
			std::ifstream is( filenames[i].c_str() ) ;
			while( is >> elt ) {
				result->add_level( elt ) ;
			}
		}
		return genfile::SampleFilter::UniquePtr( result.release() ) ;
	}

	void process_other_options() {
		m_ignore_warnings = m_options.check_if_option_was_supplied( "-force" ) ;
	}

	genfile::CohortIndividualSource::UniquePtr open_samples(
		std::vector< std::string > const& filenames
	) {
		try {
			return unsafe_open_samples( filenames ) ;
		}
		catch( genfile::MalformedInputError const& e ) {
			m_ui_context.logger() << "\n\n!! ERROR (" << e.what() << "): the sample file \"" << e.source() << "\" is malformed on line "
				<< e.line() + 1 ;
			if( e.has_column() ) {
				m_ui_context.logger() << ", column " << e.column() + 1 ;
			}
			m_ui_context.logger() << ".  Quitting.\n" ;
			throw ;
		}
	}

	genfile::CohortIndividualSource::UniquePtr unsafe_open_samples(
		std::vector< std::string > const& filenames
	) {
		genfile::CohortIndividualSourceChain::UniquePtr source_chain( new genfile::CohortIndividualSourceChain() ) ;
		for( std::size_t i = 0; i < filenames.size(); ++i ) {
			source_chain->add_source(
				genfile::CohortIndividualSource::UniquePtr(
					new genfile::CategoricalCohortIndividualSource(
						filenames[i],
						m_options.get_value< std::string >( "-missing-code" )
					)
				),
				"cohort_" + genfile::string_utils::to_string( i+1 )
			) ;
		}
		genfile::CohortIndividualSource::UniquePtr result( source_chain.release() ) ;
		return result ;
	}

	genfile::CohortIndividualSource::UniquePtr create_dummy_samples(
		genfile::SNPDataSource const& snp_data_source
	) {
		genfile::CohortIndividualSource::UniquePtr result ;
		std::vector< genfile::VariantEntry > sample_ids( snp_data_source.number_of_samples(), genfile::MissingValue() ) ;
		snp_data_source.get_sample_ids( boost::bind( &impl::set_vector_entry, &sample_ids, _1, _2 )) ;
		result.reset( new genfile::CountingCohortIndividualSource( sample_ids ) ) ;
		return result ;
	}

	genfile::CohortIndividualSource::UniquePtr open_sample_data( genfile::CohortIndividualSource::UniquePtr samples, std::vector< std::string > sample_files ) const {
		try {
			return unsafe_open_sample_data( samples, sample_files ) ;
		}
		catch( genfile::InputError const& e ) {
			m_ui_context.logger() << "\n\n!! ERROR (" << e.what() << "): " << e.format_message() << ".\n" ;
			m_ui_context.logger() << "Quitting.\n" ;
			throw ;
		}
	}

	genfile::CohortIndividualSource::UniquePtr unsafe_open_sample_data( genfile::CohortIndividualSource::UniquePtr samples, std::vector< std::string > sample_files ) const {
		genfile::DataLinkingCohortIndividualSource::UniquePtr result( new genfile::DataLinkingCohortIndividualSource( samples )) ;
		for( std::size_t i = 0; i < sample_files.size(); ++i ) {
			result->add_source(
				genfile::CohortIndividualSource::UniquePtr(
					new genfile::CategoricalCohortIndividualSource(
						sample_files[i],
						m_options.get_value< std::string >( "-missing-code" )
					)
				)
			) ;
		}
		return genfile::CohortIndividualSource::UniquePtr( result.release() )	 ;
	}

	std::vector< std::size_t > compute_excluded_samples( genfile::CohortIndividualSource::UniquePtr const& sample_source ) {
		std::vector< std::size_t > included_samples ;
		void(std::vector<std::size_t>::*push_back)(std::size_t const&) = &std::vector<std::size_t>::push_back ;
		m_sample_filter->test(
			*sample_source,
			boost::bind(
				push_back,
				&included_samples,
				_1
			),
			&m_sample_filter_diagnostic_matrix
		) ;
		std::set< std::size_t > excluded_samples(
			boost::counting_iterator< std::size_t >(0),
			boost::counting_iterator< std::size_t >( sample_source->get_number_of_individuals() )
		) ;
		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
			excluded_samples.erase( included_samples[i] ) ;
		}
		return std::vector< std::size_t >( excluded_samples.begin(), excluded_samples.end() ) ;
	}

	genfile::CohortIndividualSource::UniquePtr condition_on(
		genfile::CohortIndividualSource::UniquePtr samples,
		genfile::SNPDataSource& snps,
		std::string const& conditioning_spec
	) {
		genfile::WithSNPDosagesCohortIndividualSource::SNPDosageSpec spec ;

		using genfile::string_utils::split_and_strip ;
		using genfile::string_utils::split ;
		std::vector< std::string > elts = split_and_strip( conditioning_spec, ",", " \t" ) ;
		for( std::size_t i = 0; i < elts.size(); ++i ) {
			std::vector< std::string > parts = split( elts[i], "(" ) ;
			if( parts.size() == 1 ) {
				parts.push_back( "add" ) ;
			}
			else if( parts.size() != 2 ) {
				throw genfile::BadArgumentError( "QCToolContext::condition_on()", "conditioning_spec=\"..." + elts[i] + "...\"" ) ;
			} else if( parts[1][ parts[1].size() - 1 ] != ')' ) {
				throw genfile::BadArgumentError( "QCToolContext::condition_on()", "conditioning_spec=\"..." + elts[i] + "...\"" ) ;
			}

			if( parts[1][ parts[1].size() - 1 ] == ')' ) {
				parts[1].resize( parts[1].size() - 1 ) ;
			}

			genfile::VariantIdentifyingDataTest::SharedPtr snp_matcher(
				genfile::WithSNPDosagesCohortIndividualSource::create_snp_matcher( parts[0] ).release()
			) ;

			std::vector< std::string > types = split( parts[1], "|" ) ;
			for( std::size_t j = 0; j < types.size(); ++j ) {
				if( types[j] == "gen" ) {
					types[j] = "add" ;
					types.push_back( "het" ) ;
				}
			}

			spec[ snp_matcher ] = std::set< std::string >( types.begin(), types.end() ) ;
		}

		genfile::CohortIndividualSource::ConstUniquePtr const_samples( samples.release() ) ;
		return genfile::CohortIndividualSource::UniquePtr(
			genfile::WithSNPDosagesCohortIndividualSource::create( const_samples, snps, spec )
			.release()
		) ;
	}

	genfile::CohortIndividualSource::UniquePtr quantile_normalise_columns(
		genfile::CohortIndividualSource::UniquePtr source,
		std::vector< std::string > const& column_names
	) const {
		genfile::ValueMappingCohortIndividualSource::UniquePtr value_mapping_source(
			new genfile::ValueMappingCohortIndividualSource( source )
		) ;
		for( std::size_t i = 0; i < column_names.size(); ++i ) {
			add_normalisation( *value_mapping_source, column_names[i] ) ;
		}
		return genfile::CohortIndividualSource::UniquePtr( value_mapping_source.release() ) ;
	}

 	void add_normalisation(
		genfile::ValueMappingCohortIndividualSource& source,
		std::string const& column_name
	) const {
		genfile::CohortIndividualSource::ColumnSpec column_spec = source.get_column_spec() ;
		std::size_t column_i = column_spec.find_column( column_name ) ;
		if( !column_spec[ column_i ].is_continuous() ) {
			throw genfile::BadArgumentError( "qctool::quantile_normalise_columns()", "column_name = \"" + column_name + "\"" ) ;
		}
		genfile::CrossCohortCovariateValueMapping::UniquePtr mapping(
			new genfile::QuantileNormalisingCrossCohortCovariateValueMapping( column_name )
		) ;
		mapping->add_source( source ) ;
		source.add_mapping( column_name, column_name + ":quantile_normalised", mapping ) ;
	}

	void check_for_errors_and_warnings() {
		check_for_errors() ;
		check_for_warnings() ;
	}

	void check_for_errors() {
		if( m_mangled_options.gen_filename_mapper().input_files().size() == 0 ) {
			m_errors.push_back( "At least one GEN input file must be supplied." ) ;
		}

		for( std::size_t i = 0; i < m_mangled_options.gen_filename_mapper().input_files().size(); ++i ) {
			for( std::size_t j = 0; j < m_mangled_options.gen_filename_mapper().output_filenames().size(); ++j ) {
				if(
					strings_are_nonempty_and_equal(
						m_mangled_options.gen_filename_mapper().output_filenames()[j],
						m_mangled_options.gen_filename_mapper().input_file( i ).filename()
					) && m_mangled_options.gen_filename_mapper().output_filenames()[j] != "-"
				) {
					m_errors.push_back( "Output GEN file \"" + m_mangled_options.gen_filename_mapper().output_filenames()[j] +"\" also specified as input GEN file." ) ;
					break ;
				}
			}
			if( strings_are_nonempty_and_equal( m_mangled_options.output_sample_filename(), m_mangled_options.gen_filename_mapper().input_file( i ).filename() )) {
				m_errors.push_back( "Output SAMPLE file \"" + m_mangled_options.output_sample_filename() +"\" also specified as input GEN file." ) ;
				break ;
			}
		}
		for( std::size_t j = 0; j < m_mangled_options.gen_filename_mapper().output_filenames().size(); ++j ) {
			if( strings_are_nonempty_and_equal( m_mangled_options.output_sample_filename(), m_mangled_options.gen_filename_mapper().output_filenames()[j] )) {
				m_errors.push_back( "The GEN and SAMPLE output filenames must differ." ) ;
			}
		}
		//if( m_mangled_options.input_sample_filenames().size() == 0 && m_sample_filter->number_of_clauses() != 0 ) {
		//	m_errors.push_back( "To filter on samples, please supply input sample files." ) ;
		//}
	}

	void check_for_warnings() {
		if( m_options.check( "-osample" ) && m_mangled_options.input_sample_filenames().size() == 0 ) {
			m_warnings.push_back( "You are outputting per-sample stats, but no input sample files have been supplied.\n"
			"   Statistics will be output but the ID fields will be templates.") ;
		}
		if( m_snp_data_source->total_number_of_snps() && *m_snp_data_source->total_number_of_snps() == 0 ) {
			m_warnings.push_back( "There are no SNPs in the source files (after exclusions, translation, aligning and matching between cohorts where relevant).\n" ) ;
		}
		if( m_samples->get_number_of_individuals() == 0 ) {
			m_warnings.push_back( "There are no individuals in the source files (after exclusions, translation, aligning and matching between cohorts where relevant).\n" ) ;
		}
	}

	bool strings_are_nonempty_and_equal( std::string const& left, std::string const& right ) {
		return (!left.empty()) && (!right.empty()) && (left == right) ;
	}

} ;

struct QCToolApplication: public appcontext::ApplicationContext
{
public:
	static void check_options( appcontext::OptionProcessor const& options ) {
		// Bail out early if -no-clobber is given and we might overwrite files
		if( options.check( "-no-clobber" ) ) {
			std::vector< std::string > outputOptions ;
			outputOptions.push_back( "-og" ) ;
			outputOptions.push_back( "-os" ) ;
			outputOptions.push_back( "-osnp" ) ;
			outputOptions.push_back( "-osample" ) ;
			outputOptions.push_back( "-olog" ) ;

			for( std::size_t i = 0; i < outputOptions.size(); ++i ) {
				std::string const& option = outputOptions[i] ;
				if( options.check( option ) && boost::filesystem::exists( options.get< std::string >( option ) ) ) {
					throw genfile::BadArgumentError(
						"QCToolApplication::check_options()",
						option + "=\"" + options.get< std::string >( option ) + "\"",
						"File \"" + options.get< std::string >( option ) + "\" exists, but -no-clobber is specified."
					) ;
				}
			}
		}
	}

public:
	QCToolApplication( int argc, char** argv ):
		appcontext::ApplicationContext(
			globals::program_name,
			"version: " + globals::program_version + ", revision: " + globals::program_revision,
			std::auto_ptr< appcontext::OptionProcessor >( new QCToolOptionProcessor ),
			argc,
			argv,
			"-log",
			QCToolApplication::check_options
		)
	{
		process() ;
	}

private:

	void process() {
		try {
			unsafe_process() ;
		}
		catch( genfile::InputError const& e ) {
			get_ui_context().logger() << "!! Error (" << e.what() << "): " << e.format_message() << ".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::FileNotFoundError const& e ) {
			get_ui_context().logger() << "\nError: No file matching \"" << e.filespec() << "\" could be found.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::db::Error const& e ) {
			get_ui_context().logger() << "!! Error (" << e.what() << ") with the following statement: \""
				<< e.sql()
				<< "\".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
	}

	struct ProgressWrapper: public genfile::SNPDataSourceProcessor::Callback, boost::noncopyable {
		static UniquePtr create( appcontext::UIContext& ui_context, std::string const& name ) {
			return UniquePtr( new ProgressWrapper( ui_context, name )) ;
		}

		ProgressWrapper( appcontext::UIContext& ui_context, std::string const& name ):
			m_progress_context( ui_context.get_progress_context( name ) ),
			m_number_of_snps( 0 )
		{}

		void begin_processing_snps( std::size_t, genfile::SNPDataSource::Metadata const& ) {
			m_progress_context.restart_timer() ;
			m_progress_context.notify_progress( 0 ) ;
		}

		void processed_snp( genfile::VariantIdentifyingData const&, genfile::VariantDataReader::SharedPtr data_reader ) {
			m_progress_context.notify_progress( ++m_number_of_snps ) ;
		}
		void end_processing_snps() {
			m_progress_context.finish() ;
			std::cerr << "Total: " << m_number_of_snps << "SNPs.\n" ;
		}

	private:
		appcontext::UIContext::ProgressContext m_progress_context ;
		std::size_t m_number_of_snps ;
	} ;

	void write_samples( genfile::CohortIndividualSource const& samples, std::string const& filename, std::string format = "traditional" ) {
		using genfile::string_utils::to_string ;
		genfile::CohortIndividualSource::ColumnSpec const spec = samples.get_column_spec() ;

		assert( format == "traditional" || format == "new" ) ;

		std::ofstream str( filename.c_str() ) ;

		if( format == "new" ) {
			str << "# Written by qctool, " + appcontext::get_current_time_as_string() << "\n" ;
			str << "# Number of samples: " << to_string( samples.get_number_of_individuals() ) << "\n" ;
			str << "# metadata: {\n" ;
			for( std::size_t j = 0; j < spec.size(); ++j ) {
				str << "#   \""
					<< spec[j].name()
					<< "\": { \"type\": \""
					<< spec[j].type()
					<< "\" }"
					<< ( ( j < ( spec.size() - 1 ) ) ? ",\n" : "\n" )
				;
			}
			str << "# }\n" ;

			for( std::size_t j = 0; j < spec.size(); ++j ) {
				str << (( j > 0 ) ? " " : "" ) << spec[j].name() ;
			}
			str << "\n" ;

		} else if( format == "traditional" ) {
			for( std::size_t j = 0; j < spec.size(); ++j ) {
				str << (( j > 0 ) ? " " : "" ) << spec[j].name() ;
			}
			str << "\n" ;
			for( std::size_t j = 0; j < spec.size(); ++j ) {
				str << (( j > 0 ) ? " " : "" ) << spec[j].type() ;
			}
			str << "\n" ;
		} else {
			assert(0) ;
		}

		for( std::size_t i = 0; i < samples.get_number_of_individuals(); ++i ) {
			for( std::size_t j = 0; j < spec.size(); ++j ) {
				str << ( j > 0 ? " " : "" ) << samples.get_entry( i, spec[j].name() ) ;
			}
			str << "\n" ;
		}
	}

	std::vector< std::string > parse_filespec( std::string spec ) {
		std::vector< std::string > result(2) ;
		result[0] = "flat" ;

		// First get rid of leading sqlite specifier, if any.
		if( spec.size() >= 9 && spec.substr( 0, 9 ) == "sqlite://" ) {
			result[0] = "sqlite" ;
			spec = spec.substr( 9, spec.size() ) ;
		}

		std::vector< std::string > elts = genfile::string_utils::split( spec, ":" ) ;
		assert( elts.size() > 0 ) ;

		if( elts[0].size() >= 7 && elts[0].substr( elts[0].size() - 7, 7 ) == ".sqlite" ) {
			result[0] = "sqlite" ;
		}
		result[1] = elts[0] ;

		if( elts.size() > 2 ) {
			throw genfile::BadArgumentError(
				"QCToolProcessor::parse_filespec()",
				"spec=\"" + spec + "\"",
				"Expected format for filespec is <filename> or sqlite://<filename>[:<tablename>]."
			) ;
		}
		if( elts.size() == 2 ) {
			if( result[0] != "sqlite" ) {
				throw genfile::BadArgumentError(
					"QCToolProcessor::parse_filespec()",
					"spec=\"" + spec + "\"",
					"Table spec is not expected unless a sqlite file is specified (i.e. sqlite://<filename>:<tablename>)"
				) ;
			}
			result.push_back( elts[1] ) ;
		}

		return( result ) ;
	}

	void unsafe_process() {
		std::size_t const number_of_threads = options().get_value< std::size_t >( "-threads" ) ;
		worker::Worker::UniquePtr worker ;
		if( number_of_threads > 0 ) {
			worker.reset( new worker::QueuedMultiThreadedWorker( number_of_threads )) ;
		} else {
			worker.reset( new worker::SynchronousWorker() ) ;
		}

		QCToolCmdLineContext context(
			options(),
			get_ui_context()
		) ;

		genfile::SimpleSNPDataSourceProcessor processor ;

		qcdb::Storage::SharedPtr per_snp_storage ;
		if( SNPSummaryComponent::is_needed( options() )) {
			if( !options().check( "-osnp" )) {
				throw genfile::BadArgumentError(
					"QCToolApplication::unsafe_process()",
					"-osnp option",
					"-osnp must be supplied to output per-SNP computations."
				) ;
			}
			std::vector< std::string > const file_spec = parse_filespec( options().get< std::string >( "-osnp" ) ) ;
			if( file_spec[0] == "sqlite" ) {
				qcdb::FlatTableDBOutputter::SharedPtr table_storage = qcdb::FlatTableDBOutputter::create_shared(
					file_spec[1],
					options().get< std::string >( "-analysis-name" ),
					options().get< std::string >( "-analysis-chunk" ),
					get_application_metadata(),
					options().get< std::string >( "-compare-variants-by" )
				) ;
				if( file_spec.size() == 3 ) {
					table_storage->set_table_name( file_spec[2] ) ;
				}
				if( options().check( "-sqlite-dont-link-variants" )) {
					table_storage->set_dont_link_variants() ;
				}
				if( !options().check( "-sqlite-store-alt-identifiers" )) {
					table_storage->set_no_alt_identifiers() ;
				}
				per_snp_storage = table_storage ;
			} else {
				assert( file_spec[0] == "flat" ) ;
				per_snp_storage = qcdb::FlatFileOutputter::create_shared(
					file_spec[1],
					options().get< std::string >( "-analysis-name" ),
					options().get_values_as_map()
				) ;
			}

			SNPSummaryComponent component(
				context.get_cohort_individual_source(),
				options(),
				get_ui_context()
			) ;

			component.setup( processor, per_snp_storage ) ;
		}

		sample_stats::SampleStorage::SharedPtr per_sample_storage ;
		if( options().check( "-osample" )) {
			std::vector< std::string > const file_spec = parse_filespec( options().get< std::string >( "-osample" ) ) ;
			if( file_spec[0] == "sqlite" ) {
				// Catch the case where we write both to the same db.
				boost::optional< genfile::db::Connection::RowId > analysis_id ;
				if( options().check( "-osnp" ) && file_spec[1] == options().get< std::string >( "-osnp" ) && per_snp_storage ) {
					analysis_id = per_snp_storage->analysis_id() ;
				}
				sample_stats::FlatTableDBOutputter::SharedPtr
					dbStorage = sample_stats::FlatTableDBOutputter::create_shared(
						file_spec[1],
						options().get< std::string >( "-analysis-name" ),
						options().get< std::string >( "-analysis-chunk" ),
						options().get_values_as_map(),
						context.get_cohort_individual_source(),
						analysis_id
					) ;
				if( file_spec.size() == 3 ) {
					dbStorage->set_table_name( file_spec[2] ) ;
				}
				per_sample_storage = dbStorage ;
			} else {
				assert( file_spec[0] == "flat" ) ;
				per_sample_storage = sample_stats::FlatFileOutputter::create_shared(
					context.get_cohort_individual_source(),
					file_spec[1],
					options().get< std::string >( "-analysis-name" ),
					options().get_values_as_map()
				) ;
			}
		}

		if( SampleSummaryComponent::is_needed( options() ) ) {
			if( !per_sample_storage ) {
				throw genfile::BadArgumentError(
					"QCToolApplication::unsafe_process()",
					"-osample",
					"An output sample data file (-osample) must be supplied to output per-sample computations."
				) ;
			}

			SampleSummaryComponent::UniquePtr sample_summary_component
				= SampleSummaryComponent::create(
					options(),
					context.get_cohort_individual_source(),
					get_ui_context()
				) ;
			sample_summary_component->setup( processor, per_sample_storage ) ;
		}

		std::auto_ptr< DataReadTest > data_read_test ;
		if( options().check_if_option_was_supplied( "-read-test" )) {
			data_read_test.reset( new DataReadTest() ) ;
			processor.add_callback( *data_read_test ) ;
		}

		if( options().check_if_option_was_supplied_in_group( "Kinship options" )) {
			RelatednessComponent::UniquePtr relatedness_component = RelatednessComponent::create(
				options(),
				context.get_cohort_individual_source(),
				worker.get(),
				get_ui_context()
			) ;
			relatedness_component->setup( processor, per_sample_storage ) ;
		}

		HaplotypeFrequencyComponent::UniquePtr haplotype_frequency_component ;
		if( options().check_if_option_was_supplied_in_group( "LD computation options" )) {
			std::vector< std::string > filenames = options().get_values< std::string >( "-compute-ld-with" ) ;

			genfile::CohortIndividualSource::UniquePtr ld_samples ;
			genfile::SNPDataSource::UniquePtr ld_variants ;

			if( filenames.size() != 2 ) {
				throw genfile::BadArgumentError(
					"QCToolApplication::unsafe_process()",
					"-compute-ld-with\"" + genfile::string_utils::join( filenames, " " ) + "\"",
					"Value should be two filenames - a genotype file and a sample file."
				) ;
			}

			context.open_data_sources(
				std::vector< std::vector< genfile::wildcard::FilenameMatch > >( 1,
						genfile::wildcard::find_files_by_chromosome(
						filenames[0],
						genfile::wildcard::eNON_SEX_CHROMOSOMES
					)
				),
				std::vector< std::string >( 1, filenames[1] ),
				boost::optional< std::vector< std::string > >(),
				&ld_variants,
				&ld_samples
			) ;

			std::vector< std::string > id_columns
				= genfile::string_utils::split_and_strip( options().get_value< std::string >( "-match-sample-ids" ), "~" ) ;
			if( id_columns.size() != 2 ) {
				throw genfile::BadArgumentError(
					"QCToolApplication::unsafe_process()",
					"-match-sample-ids=\"" + options().get_value< std::string >( "-match-sample-ids" ) + "\"",
					"Value should be two strings separated by ~."
				) ;
			}

			haplotype_frequency_component = HaplotypeFrequencyComponent::create(
				context.get_cohort_individual_source(),
				id_columns[0],
				ld_samples,
				id_columns[1],
				ld_variants,
				options(),
				get_ui_context()
			) ;

			processor.add_callback( *haplotype_frequency_component ) ;
		}

		if( options().check( "-og" ) || options().check( "-op" ) ) {
			SNPOutputComponent component(
				context.get_cohort_individual_source(),
				options(),
				get_ui_context()
			) ;
			component.setup(
				context.fltrd_in_snp_data_sink(),
				processor
			) ;
		}

		if( options().check( "-os" )) {
			write_samples(
				context.get_cohort_individual_source(),
				options().get< std::string >( "-os" ),
				options().get< std::string >( "-output-sample-format" )
			) ;
		}
		// Process it (but only if there was something to do) !
		if( processor.get_callbacks().size() > 0 ) {
			processor.add_callback( ProgressWrapper::create( get_ui_context(), "Processing SNPs" ) ) ;
			processor.process( context.snp_data_source() ) ;
		} else {
			get_ui_context().logger() << "SNPs do not need to be visited -- skipping.\n" ;
		}
		if( per_sample_storage ) {
			per_sample_storage->finalise() ;
		}
		if( per_snp_storage ) {
			per_snp_storage->finalise() ;
		}
	}
} ;

int main( int argc, char** argv ) {
    // needed for efficient output to standard output streams
    std::ios_base::sync_with_stdio( false ) ;
    try {
		QCToolApplication app( argc, argv ) ;
    }
	catch( appcontext::HaltProgramWithReturnCode const& e ) {
		return e.return_code() ;
	}
	return 0 ;
}

max_stars_repo_name: gavinband/qctool

max_issues_repo_name: gavinband/qctool

max_forks_repo_name: gavinband/qctool