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)
  
  #include <string>
  #include <iostream>
  #include <limits>
  #include <typeinfo>
  #include <fstream>
  
  #include <boost/filesystem.hpp>
  #include <boost/format.hpp>
  #include <boost/math/distributions/normal.hpp>
  #include <boost/math/distributions/chi_squared.hpp>
  #include <boost/timer/timer.hpp>
  #include <boost/ptr_container/ptr_vector.hpp>
  #include <boost/bind.hpp>
  #include <boost/foreach.hpp>
  #include <boost/iterator/counting_iterator.hpp>
  
  #include "config/package_revision_autogenerated.hpp"
  
  #include "appcontext/CmdLineOptionProcessor.hpp"
  #include "appcontext/ApplicationContext.hpp"
  
  #include "genfile/SNPDataSource.hpp"
  #include "genfile/SNPDataSourceChain.hpp"
  #include "genfile/ThreshholdingSNPDataSource.hpp"
  #include "genfile/ToGP.hpp"
  #include "genfile/CommonSNPFilter.hpp"
  #include "genfile/VariantIdentifyingDataFilteringSNPDataSource.hpp"
  #include "genfile/CrossCohortCovariateValueMapping.hpp"
  #include "genfile/SampleFilter.hpp"
  #include "genfile/SampleFilteringSNPDataSource.hpp"
  #include "genfile/SampleFilteringCohortIndividualSource.hpp"
  #include "genfile/DataLinkingCohortIndividualSource.hpp"
  #include "genfile/CompoundSampleFilter.hpp"
  #include "genfile/SampleFilterNegation.hpp"
  #include "genfile/VariableInSetSampleFilter.hpp"
  #include "genfile/GPThresholdingGTSetter.hpp"
  #include "genfile/db/Error.hpp"
  #include "genfile/CartesianProductVisitor.hpp"
  
  #include "components/SNPSummaryComponent/InfoComputation.hpp"
  
  #include "metro/constants.hpp"
  #include "metro/regression/Design.hpp"
  #include "metro/SampleRange.hpp"
  #include "metro/intersect_ranges.hpp"
  #include "metro/regression/BinomialLogistic.hpp"
  #include "metro/regression/ThreadedLogLikelihood.hpp"
  #include "metro/regression/IndependentNormalWeightedLogLikelihood.hpp"
  #include "metro/regression/IndependentLogFWeightedLogLikelihood.hpp"
  #include "metro/regression/LogPosteriorDensity.hpp"
  #include "metro/ValueStabilisesStoppingCondition.hpp"
  #include "metro/Snptest25StoppingCondition.hpp"
  #include "metro/maximisation.hpp"
  #include "metro/fit_model.hpp"
  #include "metro/CholeskyStepper.hpp"
  #include "metro/SampleStratification.hpp"
  
  #include "qcdb/MultiVariantStorage.hpp"
  #include "qcdb/FlatTableDBOutputter.hpp"
  
  #include <Eigen/Core>
  #include <Eigen/QR>
  #include "Eigen/Eigenvalues"
  
  #define DEBUG 0
  
  namespace globals {
  	std::string const program_name = "hptest" ;
  	std::string const program_version = package_version ;
  	std::string const program_revision =  std::string( package_revision ).substr( 0, 7 ) ;
  }
  
  namespace {
  	bool in_set( std::string const& arg, std::vector< std::string > const& values, std::set< std::string > const& permitted ) {
  		for( auto v: values ) {
  			if( permitted.find( v ) == permitted.end() ) {
  				throw appcontext::OptionValueInvalidException(
  					arg,
  					values,
  					arg + ": unexpected value (\"" + v + "\")"
  				) ;
  			}
  		}
  		return true ;
  	}
  	
  	std::vector< std::string > collect_unique_ids( std::vector< std::string > const& ids_or_filenames ) {
  		std::vector< std::string > result ;
  		for( auto elt: ids_or_filenames ) {
  			if( boost::filesystem::exists( elt )) {
  				std::ifstream f( elt ) ;
  				std::copy(
  					std::istream_iterator< std::string >( f ),
  					std::istream_iterator< std::string >(),
  					std::back_inserter< std::vector< std::string > >( result )
  				) ;
  			} else {
  				result.push_back( elt ) ;
  			}
  		}
  		// now sort and uniqueify them...
  		std::sort( result.begin(), result.end() ) ;
  		std::vector< std::string >::const_iterator newBack = std::unique( result.begin(), result.end() ) ;
  		result.resize( newBack - result.begin() ) ;
  		return result ;
  	}
  
  	double extract_mapped_continuous_value(
  		std::size_t sample_index,
  		genfile::CohortIndividualSource const& samples,
  		std::size_t column_index,
  		genfile::CrossCohortCovariateValueMapping const& mapping
  	) {
  		double value = metro::NA ;
  		genfile::VariantEntry entry = samples.get_entry( sample_index, column_index ) ;
  		if( !entry.is_missing() ) {
  			value = mapping.get_mapped_value( entry ).as< double >() ;
  		}
  		return value ;
  	}
  
  	int extract_mapped_categorical_value(
  		std::size_t sample_index,
  		genfile::CohortIndividualSource const& samples,
  		std::size_t column_index,
  		genfile::CrossCohortCovariateValueMapping const& mapping
  	) {
  		int value = -1 ;
  		genfile::VariantEntry entry = samples.get_entry( sample_index, column_index ) ;
  		if( !entry.is_missing() ) {
  			value = mapping.get_mapped_value( entry ).as< int >() ;
  		}
  		return value ;
  	}
  
  	std::string get_unmapped_level(
  		genfile::CrossCohortCovariateValueMapping const& mapping,
  		int level
  	) {
  		genfile::VariantEntry entry = mapping.get_unmapped_value( level ) ;
  		return entry.as< std::string >() ;
  	}
  }
  
  struct HPTestOptions: 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[ "-predictor" ]
  	        .set_description( 	"Path to predictor genotype files."
  								"The given filename may contain the wildcard character '#', which expands to match a"
  								"one- or two-character chromosome identifier." )
  			.set_takes_values( 1 )
  			.set_minimum_multiplicity( 1 )
  			.set_maximum_multiplicity( 1 ) ;
  		options[ "-incl-predictor-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[ "-incl-predictor-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[ "-outcome" ]
  			.set_description(
  				"Path to outcome genotype files."
  				"The given filename may contain the wildcard character '#', which expands to match a"
  				"one- or two-character chromosome identifier." )
  			.set_takes_values( 1 )
  			.set_minimum_multiplicity( 1 )
  			.set_maximum_multiplicity( 1 ) ;
  		options[ "-outcome-name" ]
  			.set_description( "Specify a name for the outcome variable" )
  			.set_takes_single_value()
  			.set_default_value( "outcome" ) ;
  		options[ "-incl-outcome-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[ "-incl-outcome-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[ "-treat-outcome-as-haploid" ]
  			.set_description( "Turn outcome homozygous genotype calls into haploid calls. "
  				"Heterozygous calls will be treated as missing" ) ;
  		options[ "-treat-predictor-as-haploid" ]
  			.set_description( "Turn predictor homozygous genotype calls or probabilities into haploid calls or probabilities. "
  				"Samples with nonzero probability of heterozygous genotype will be treated as missing." ) ;
  		options[ "-s" ]
  			.set_description( "Path of sample file" )
  			.set_takes_values( 1 )
  			.set_minimum_multiplicity( 1 )
  			.set_maximum_multiplicity( 1 ) ;
  
  		options[ "-sample-data" ]
  			.set_description( "Path of additional sample file(s) (in the same format(s) as accepted by -s) to read. "
  				"These files are used to look up additional sample information, using a sample ID in the first column as key. "
  				"For each sample and column, the relevant entry will be taken from the superset of all sample files supplied. "
  				"(It is an error to supply two different, non-missing values for columns with the same name for the same sample.)"
  			)
  			.set_takes_values_until_next_option()
  			.set_minimum_multiplicity( 0 )
  			.set_maximum_multiplicity( 100 ) ;
  
  		options.declare_group( "Output file options" ) ;
  		options[ "-o" ]
  			.set_description( "Output file" )
  			.set_takes_values( 1 )
  			.set_minimum_multiplicity( 1 )
  			.set_maximum_multiplicity( 1 ) ;
  		options[ "-output-models" ]
  			.set_description( "Specify which models to output estimates, standard errors, and Wald P-values for. "
  				"This can be either \"alternative\", in which case only non-null models are output, or \"all\" in which "
  				"case all models are output. "
  			)
  			.set_takes_single_value()
  			.set_default_value( "alternative" )
  			.set_check( boost::bind( &in_set, _1, _2, std::set< std::string >({"all", "alternative"}) ) ) ;
  		options[ "-output-parameters" ]
  			.set_description( "Specify which parameters to output estimates, standard errors, and Wald P-values for."
  				"By default only data on the genetic parameters are output; specify \"all\" to output all parameters."
  			)
  			.set_takes_single_value()
  			.set_default_value( "genetic" )
  			.set_check( boost::bind( &in_set, _1, _2, std::set< std::string >({"all", "genetic"}) ) ) ;
  		options[ "-output-counts" ]
  			.set_description( "Specify that hptest should output predictor and outcome genotype counts" ) ;
  		options[ "-output-all-variants" ]
  			.set_description( "Specify that hptest should output all variant combinations, even those where tests are skipped." ) ;
  		
  		options.declare_group( "Model options" ) ;
  		options[ "-model" ]
  			.set_description( "Specify models to fit.  Syntax is <name>:<model>" )
  			.set_takes_values_until_next_option()
  			.set_default_value( "gen:add+het" ) ;
  		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 ) ;
  		options[ "-minimum-predictor-count" ]
  			.set_description( "Skip predictors that have less than this minor allele count" )
  			.set_takes_values(1)
  			.set_default_value( 10 ) ;
  		options[ "-minimum-outcome-count" ]
  			.set_description( "Skip outcomes that have less than this count" )
  			.set_takes_values(1)
  			.set_default_value( 10 ) ;
  		options[ "-prior" ]
  			.set_description( "Specify a prior for bayesian computation."
  				"Prior spec should be in the form <parameter name>~<distribution family>(parameters)."
  				"(It is recommended to put single quotes round the whole argument to protect from shell expansion.)"
  				"Currently, logF and gaussian priors are supported."
  			)
  			.set_takes_values_until_next_option()
  			.set_default_value( "add/[stratum]/[outcome]=1~logf(2,2)" )
  			.set_default_value( "overdominance/[stratum]/[outcome]=1~logf(4,4)" )
  			.set_default_value( "het/[stratum]/[outcome]=1~logf(2,2)" )
  			.set_default_value( "dom/[stratum]/[outcome]=1~logf(2,2)" )
  			.set_default_value( "rec/[stratum]/[outcome]=1~logf(2,2)" )
  		;
  		options[ "-no-prior" ]
  			.set_description( "Specify that no prior should be used.  Only frequentist summaries will be output." ) ;
  		options[ "-no-covariate-priors" ]
  			.set_description( "Specify that no default prior should be placed on covariates. "
  				"(By default a N(0,1000) prior is used)." ) ;
  		options[ "-covariates" ]
  			.set_description( "Specify the name of one or more covariates to include in the model."
  				" These must be columns named in the sample file." )
  			.set_takes_values_until_next_option() ;
  		options[ "-stratify" ]
  			.set_description( "Specify the name of a discrete covariate to stratify across" )
  			.set_takes_single_value() ;
  		options.declare_group( "Model fitting options" ) ;
  		options[ "-tolerance" ]
  			.set_description( "Tolerance" )
  			.set_takes_values(1)
  			.set_default_value( 0.001 ) ;
  		options[ "-max-iterations" ]
  			.set_description( "Maximum fitting iterations" )
  			.set_takes_values(1)
  			.set_default_value( 100 ) ;
  
  		options.option_excludes_option( "-prior", "-no-prior" ) ;
  
  		options.declare_group( "Output options" ) ;
  		options[ "-analysis-name" ]
  			.set_description( "Specify a name to label results from this analysis with." )
  			.set_takes_single_value()
  			.set_default_value( "hptest 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[ "-analysis-id" ]
  			.set_description( "Specify an integer ID for the current analysis." )
  			.set_takes_single_value() ;
  		options[ "-outcome-genotype-call-threshold" ]
  			.set_description( "The threshold to apply to outcome genotype probabilities, if necessary, to make calls." )
  			.set_takes_single_value()
  			.set_default_value( 0.9 ) ;
  		options[ "-compare-variants-by" ]
  			.set_description( "A comma-separated list of fields to compare variants by."
  				" This affects how results for variants"
  				" (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[ "-debug" ]
  			.set_description( "Output debugging information." ) ;
  		options[ "-threads" ]
  			.set_description( "Number of additional threads to use in likelihood computations."
  				" The value 0 indicates that all work will take place in the main thread." )
  			.set_takes_single_value()
  			.set_default_value(0) ;
  	}
  } ;
  
  
  namespace {
  	// counts the number of each allele carried by each sample
  	// Outputs count of 1st allele in column 0 of output genotypes,
  	// 2nd allele in column 1, etc.
  	struct AlleleCounter: public genfile::VariantDataReader::PerSampleSetter {
  		typedef std::vector< metro::SampleRange > SampleRanges ;
  		
  		AlleleCounter(
  			Eigen::MatrixXd* genotypes,
  			Eigen::VectorXd* ploidy,
  			SampleRanges* nonmissing_samples
  		):
  			m_genotypes( genotypes ),
  			m_ploidy( ploidy ),
  			m_nonmissing_samples( nonmissing_samples ),
  			m_sample_i(0),
  			m_coerce_to_haploid( false )
  		{
  			assert( genotypes != 0 && nonmissing_samples != 0 && ploidy != 0 ) ;
  		}
  
  		void set_coerce_to_haploid() { m_coerce_to_haploid = true ; }
  
  		void initialise( std::size_t nSamples, std::size_t nAlleles ) {
  			if( nAlleles != 2 ) {
  				throw genfile::BadArgumentError(
  					"AlleleCounter::initialise()",
  					"nAlleles=" + genfile::string_utils::to_string( nAlleles ),
  					"I only support biallelic variants"
  				) ;
  			}
  			m_genotypes->resize( nSamples, nAlleles ) ;
  			m_genotypes->setZero() ;
  
  			m_ploidy->resize( nSamples ) ;
  			m_ploidy->setZero() ;
  
  			m_nonmissing_samples->clear() ;
  			m_last_nonmissing_sample_i = 0 ;
  		}
  
  		bool set_sample( std::size_t n ) {
  			m_sample_i = n ;
  			return true ;
  		}
  
  		void set_number_of_entries(
  			uint32_t ploidy, std::size_t n,
  			genfile::OrderType const order_type,
  			genfile::ValueType const value_type
  		) {
  			if( value_type != genfile::eAlleleIndex ) {
  				throw genfile::BadArgumentError(
  					"AlleleCounter::set_number_of_entries()",
  					"value_type=" + genfile::string_utils::to_string( value_type ),
  					"Expected a hard-called genotype, consider using threshholded calls."
  				) ;
  			}
  			if( order_type != genfile::ePerOrderedHaplotype && order_type != genfile::ePerUnorderedHaplotype ) {
  				throw genfile::BadArgumentError(
  					"AlleleCounter::set_number_of_entries()",
  					"order_type=" + genfile::string_utils::to_string( order_type ),
  					"Expected a hard-called genotype, consider using threshholded calls."
  				) ;
  			}
  			if( m_coerce_to_haploid ) {
  				(*m_ploidy)(m_sample_i) = 1 ;
  			} else {
  				(*m_ploidy)(m_sample_i) = ploidy ;
  			}
  		}
  
  		void set_value( std::size_t entry_i, genfile::MissingValue const value ) {
  			// This sample has missing data, end our current sample range here.
  			record_missing_sample() ;
  		}
  		
  		void set_value( std::size_t entry_i, Integer const value ) {
  			// Only accumulate if not already set missing
  			if( m_sample_i >= m_last_nonmissing_sample_i ) {
  				if( m_coerce_to_haploid && entry_i > 0 ) {
  					// Set anything that isn't homozygous to missing
  					if( (*m_genotypes)(m_sample_i,value) != 1 ) {
  						(*m_genotypes).row(m_sample_i).setZero() ;
  						record_missing_sample() ;
  					}
  				} else {
  					(*m_genotypes)(m_sample_i,value) += 1 ;
  				}
  			}
  		}
  
  		void set_value( std::size_t entry_i, double const value ) {
  			assert(0) ; // expecting GT field
  		}
  
  		void record_missing_sample() {
  			if( m_sample_i > m_last_nonmissing_sample_i ) {
  				m_nonmissing_samples->push_back(
  					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
  				) ;
  			}
  			m_last_nonmissing_sample_i = m_sample_i + 1 ;
  		}
  
  		void finalise() {
  			// add the final sample range if needed.
  			if( m_sample_i+1 > m_last_nonmissing_sample_i ) {
  				m_nonmissing_samples->push_back(
  					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i+1 )
  				) ;
  			}
  		}
  		
  	private:
  		Eigen::MatrixXd* m_genotypes ;
  		Eigen::VectorXd* m_ploidy ;
  		SampleRanges* m_nonmissing_samples ;
  		std::size_t m_last_nonmissing_sample_i ;
  		std::size_t m_sample_i ;
  		bool m_coerce_to_haploid ;
  	} ;
  	
  	struct ProbSetter: public genfile::VariantDataReader::PerSampleSetter {
  		typedef std::vector< metro::SampleRange > SampleRanges ;
  
  		ProbSetter(
  			Eigen::MatrixXd* genotypes,
  			Eigen::VectorXi* ploidy,
  			SampleRanges* nonmissing_samples
  		):
  			m_genotypes( genotypes ),
  			m_ploidy( ploidy ),
  			m_coerce_to_haploid( false ),
  			m_nonmissing_samples( nonmissing_samples ),
  			m_sample_i(0),
  			m_number_of_entries(0)
  		{
  			assert( genotypes != 0 && nonmissing_samples != 0 && ploidy != 0 ) ;
  		}
  
  		void set_coerce_to_haploid() {
  			m_coerce_to_haploid = true ;
  		}
  
  		void initialise( std::size_t nSamples, std::size_t nAlleles ) {
  			if( nAlleles != 2 ) {
  				throw genfile::BadArgumentError(
  					"ProbSetter::initialise()",
  					"nAlleles=" + genfile::string_utils::to_string( nAlleles ),
  					"I only support biallelic variants"
  				) ;
  			}
  			// This currently supports constant ploidy only, up to ploidy = 2
  			// Allocate enough space for this with biallelic variants
  			m_genotypes->resize( nSamples, 3 ) ;
  			m_genotypes->setZero() ;
  			m_ploidy->resize( nSamples ) ;
  			m_ploidy->setZero() ;
  			
  			m_nonmissing_samples->clear() ;
  			m_last_nonmissing_sample_i = 0 ;
  		}
  
  		bool set_sample( std::size_t n ) {
  			m_sample_i = n ;
  			return true ;
  		}
  
  		void set_number_of_entries(
  			uint32_t ploidy, std::size_t number_of_entries,
  			genfile::OrderType const order_type,
  			genfile::ValueType const value_type
  		) {
  			if( value_type != genfile::eProbability ) {
  				throw genfile::BadArgumentError(
  					"Probsetter::set_number_of_entries()",
  					"value_type=" + genfile::string_utils::to_string( value_type ),
  					"Expected genotype call probabilities (e.g. GP field)."
  				) ;
  			}
  			if( order_type != genfile::ePerUnorderedGenotype ) {
  				throw genfile::BadArgumentError(
  					"ProbSetter::set_number_of_entries()",
  					"order_type=" + genfile::string_utils::to_string( order_type ),
  					"Expected genotype call probabilities (e.g. GP field)."
  				) ;
  			}
  			if( ploidy > 2 ) {
  				throw genfile::BadArgumentError(
  					"ProbSetter::set_number_of_entries()",
  					"order_type=" + genfile::string_utils::to_string( order_type ),
  					"Expected an at-most diploid sample."
  				) ;
  			}
  			(*m_ploidy)(m_sample_i) = m_coerce_to_haploid ? 1 : ploidy ;
  			m_number_of_entries = number_of_entries ;
  		}
  
  		void set_value( std::size_t entry_i, genfile::MissingValue const value ) {
  			// This sample has missing data, end our current sample range here.
  			record_missing_sample() ;
  		}
  
  		void set_value( std::size_t entry_i, Integer const value ) {
  			assert(0) ; // expected a probabilty
  		}
  
  		void set_value( std::size_t entry_i, double const value ) {
  			if( m_sample_i >= m_last_nonmissing_sample_i && value != 0.0 ) {
  				if( m_coerce_to_haploid ) {
  					int const index = (entry_i / (m_number_of_entries-1)) ;
  					if( ( entry_i == 0 || (entry_i+1) == m_number_of_entries )) {
  						(*m_genotypes)(m_sample_i, index) = value ;
  					} else {
  						record_missing_sample() ;
  					}
  				} else {
  					(*m_genotypes)(m_sample_i, entry_i) = value ;
  				}
  			}
  		}
  
  		void finalise() {
  			// Add the final sample range if needed.
  			if( (m_sample_i+1) > m_last_nonmissing_sample_i ) {
  				m_nonmissing_samples->push_back(
  					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i+1 )
  				) ;
  			}
  			
  			// For probability data, we need to handle the case of missingness encoded as zero probabilities.
  			// Do that now by inspecting the genotypes.
  			{
  				std::vector< metro::SampleRange > ranges ;
  				std::size_t last_nonmissing_sample_i = 0 ;
  				std::size_t i = 0 ;
  				for( ; i < m_genotypes->rows(); ++i ) {
  					if( m_genotypes->row(i).sum() == 0.0 ) {
  						if( i > last_nonmissing_sample_i ) {
  							ranges.push_back( metro::SampleRange( last_nonmissing_sample_i, i )) ;
  						}
  						last_nonmissing_sample_i = i+1 ;
  					}
  				}
  				if( i > last_nonmissing_sample_i ) {
  					ranges.push_back( metro::SampleRange( last_nonmissing_sample_i, i )) ;
  				}
  				
  				*m_nonmissing_samples = metro::impl::intersect_ranges( *m_nonmissing_samples, ranges ) ;
  			}
  		}
  		
  	private:
  		Eigen::MatrixXd* m_genotypes ;
  		Eigen::VectorXi* m_ploidy ;
  		bool m_coerce_to_haploid ;
  		SampleRanges* m_nonmissing_samples ;
  		std::size_t m_last_nonmissing_sample_i ;
  		std::size_t m_sample_i ;
  		std::size_t m_number_of_entries ;
  		
  		void record_missing_sample() {
  			if( m_sample_i > m_last_nonmissing_sample_i ) {
  				// remember the current range of nonmissing samples
  				m_nonmissing_samples->push_back(
  					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
  				) ;
  			}
  			// Record this sample as missing
  			m_last_nonmissing_sample_i = m_sample_i + 1 ;
  		}
  	} ;
  }
  
  struct HPTestApplication: public appcontext::ApplicationContext
  {
  public:
  	
  	typedef Eigen::MatrixXd Matrix ;
  	typedef Eigen::VectorXd Vector ;
  	typedef Eigen::VectorXi IntegerVector ;
  	typedef Eigen::RowVectorXd RowVector ;
  	
  public:
  	HPTestApplication( int argc, char** argv ):
  		appcontext::ApplicationContext(
  			globals::program_name,
  			"version: " + globals::program_version + ", revision: " + globals::program_revision,
  			std::auto_ptr< appcontext::OptionProcessor >( new HPTestOptions ),
  			argc,
  			argv,
  			"-log"
  		)
  	{
  		process() ;
  	}
  	
  private:
  
  	metro::concurrency::threadpool::UniquePtr m_pool ;
  	metro::SampleStratification m_stratification ;
  	std::vector< std::string > m_prior_specs ;
  	Eigen::MatrixXd m_regression_predictor_probs ;
  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 ) ;
  		}
  	}
  	
  	void unsafe_process() {
  		using genfile::string_utils::to_string ;
  		
  		genfile::SampleFilterConjunction::UniquePtr sample_filter = get_sample_filter() ;
  
  		genfile::CohortIndividualSource::UniquePtr
  			samples = genfile::CohortIndividualSource::create( options().get< std::string >( "-s" ) ) ;
  
  		if( options().check( "-sample-data" )) {
  			samples = open_sample_data( samples, options().get_values< std::string > ( "-sample-data" )) ;
  		}
  
  		std::set< std::size_t > const excluded_samples = compute_excluded_samples( *sample_filter, *samples ) ;
  		if( excluded_samples.size() > 0 ) {
  			samples.reset(
  				genfile::SampleFilteringCohortIndividualSource::create( samples, excluded_samples ).release()
  			) ;
  		}
  
  		if( options().check( "-stratify" )) {
  			m_stratification = compute_stratification( *samples, options().get< std::string >( "-stratify" )) ;
  			std::cerr << "Stratification:\n"
  				<< m_stratification << ".\n" ;
  		}
  
  		genfile::SNPDataSource::UniquePtr host = open_genotype_data_sources(
  			options().get< std::string >( "-predictor" ),
  			get_variant_filter(
  				"-incl-predictor-range",
  				"-incl-predictor-rsids"
  			),
  			excluded_samples
  		) ;
  
  		genfile::SNPDataSource::UniquePtr para = open_genotype_data_sources(
  			options().get< std::string >( "-outcome" ),
  			get_variant_filter(
  				"-incl-outcome-range",
  				"-incl-outcome-rsids"
  			),
  			excluded_samples
  		) ;
  
  		std::size_t const N = samples->size() ;
  		if( host->number_of_samples() != N ) {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::unsafe_process()",
  				"-predictor \"" + options().get< std::string >( "-predictor" ) + "\"",
  				"Wrong number of samples (" + to_string(host->number_of_samples()) + ", expected " + to_string(N) + ")"
  			) ;
  		}
  		if( para->number_of_samples() != N ) {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::unsafe_process()",
  				"-outcome \"" + options().get< std::string >( "-outcome" ) + "\"",
  				"Wrong number of samples (" + to_string(para->number_of_samples()) + ", expected " + to_string(N) + ")"
  			) ;
  		}
  		
  		if( !options().check( "-no-prior" )) {
  			m_prior_specs = expand_prior_specs( options().get_values< std::string >( "-prior" )) ;
  		}
  	
  		write_preamble( *host, *para, *samples ) ;
  		
  		{
  			boost::optional< genfile::db::Connection::RowId > analysis_id ;
  			if( options().check( "-analysis-id" ) ) {
  				analysis_id = options().get< genfile::db::Connection::RowId >( "-analysis-id" ) ;
  			}
  			qcdb::MultiVariantStorage::UniquePtr storage = qcdb::MultiVariantStorage::create(
  				options().get< std::string > ( "-o" ),
  				2,
  				options().get< std::string > ( "-analysis-name" ),
  				options().get< std::string > ( "-analysis-chunk" ),
  				get_application_metadata(),
  				options().get< std::string > ( "-compare-variants-by" ),
  				analysis_id
  			) ;
  			storage->set_variant_names( std::vector< std::string >({ "predictor", "outcome" })) ;
  			
  			genfile::CartesianProductVisitor visitor ;
  			visitor.add_source( "predictor", host.get() ) ;
  			visitor.add_source( "outcome", para.get() ) ;
  			
  			test( visitor, *samples, *storage ) ;
  			storage->finalise() ;
  		}
  	}
  	
  	genfile::CohortIndividualSource::UniquePtr 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::create( sample_files[i] )) ;
  		}
  		return genfile::CohortIndividualSource::UniquePtr( result.release() ) ;
  	}
  	
  	genfile::CommonSNPFilter::UniquePtr get_variant_filter(
  		std::string const& incl_range_option,
  		std::string const& incl_ids_option
  	) const {
  		genfile::CommonSNPFilter::UniquePtr result ;
  		if( options().check( incl_range_option ) || options().check( incl_ids_option ) ) {
  			result.reset( new genfile::CommonSNPFilter() ) ;
  			if( options().check( incl_range_option )) {
  				std::vector< std::string > specs = collect_unique_ids( options().get_values< std::string >( incl_range_option ) ) ;
  				for ( std::size_t i = 0; i < specs.size(); ++i ) {
  					result->include_snps_in_range(
  						genfile::GenomePositionRange::parse( specs[i] )
  					) ;
  				}
  			}
  			if( options().check( incl_ids_option )) {
  				std::vector< std::string > files = collect_unique_ids( options().get_values< std::string > ( incl_ids_option ) ) ;
  				result->include_snps_in_set(
  					std::set< std::string >( files.begin(), files.end() ),
  					genfile::CommonSNPFilter::RSIDs
  				) ;
  			}
  		}
  		return result ;
  	}
  	
  	genfile::SampleFilterConjunction::UniquePtr get_sample_filter() const {
  		genfile::SampleFilterConjunction::UniquePtr filter( new genfile::SampleFilterConjunction() ) ;
  		genfile::SampleFilterDisjunction::UniquePtr sample_inclusion_filter( new genfile::SampleFilterDisjunction() ) ;
  
  		if( options().check_if_option_was_supplied( "-incl-samples" ) ) {
  			genfile::SampleFilter::UniquePtr id1_filter = get_sample_id_filter(
  				options().get_values< std::string >( "-incl-samples" )
  			) ;
  			sample_inclusion_filter->add_clause( genfile::SampleFilter::UniquePtr( id1_filter.release() ) ) ;
  		}
  
  		if( options().check_if_option_was_supplied( "-excl-samples" ) ) {
  			genfile::SampleFilter::UniquePtr id1_filter = get_sample_id_filter(
  				options().get_values< std::string >( "-excl-samples" )
  			) ;
  			filter->add_clause(
  				genfile::SampleFilter::UniquePtr(
  					new genfile::SampleFilterNegation(
  						id1_filter
  					)
  				)
  			) ;
  		}
  
  		if( options().check_if_option_was_supplied( "-incl-samples-where" ) ) {
  			std::vector< std::string > conditions = 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( options().check_if_option_was_supplied( "-excl-samples-where" ) ) {
  			std::vector< std::string > conditions = 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 ;
  	}
  
  	std::set< std::size_t > compute_excluded_samples( genfile::SampleFilter const& filter, genfile::CohortIndividualSource const& samples ) const {
  		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 ;
  		filter.test(
  			samples,
  			boost::bind(
  				push_back,
  				&included_samples,
  				_1
  			)
  		) ;
  		std::set< std::size_t > excluded_samples(
  			boost::counting_iterator< std::size_t >(0),
  			boost::counting_iterator< std::size_t >( samples.get_number_of_individuals() )
  		) ;
  		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
  			excluded_samples.erase( included_samples[i] ) ;
  		}
  		return excluded_samples ;
  	}
  
  	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() ) ;
  	}
  	
  	metro::SampleStratification compute_stratification(
  		genfile::CohortIndividualSource const& samples,
  		std::string const& column
  	) const {
  		using genfile::string_utils::to_string ;
  
  		metro::SampleStratification result ;
  		genfile::CohortIndividualSource::ColumnSpec const& column_spec = samples.get_column_spec() ;
  
  		if( !column_spec[ column ].is_discrete() ) {
  			throw genfile::BadArgumentError(
  				"PerVariantComputationManager::get_stratification()",
  				"column=\"" + column + "\"",
  				"argument to -stratify must be a discrete covariate in the sample file."
  			) ;
  		}
  		// Add strata.
  		// We do this via the a value mapping in order to ensure that values are taken in sorted order.
  		genfile::CrossCohortCovariateValueMapping::UniquePtr mapping
  			= genfile::CrossCohortCovariateValueMapping::create( column_spec[ column ] ) ;
  		mapping->add_source( samples ) ;
  		int const highest_level = int( mapping->get_number_of_distinct_mapped_values() ) ;
  		for( int level = 0; level < highest_level; ++level ) {
  			genfile::VariantEntry const unmapped_value = mapping->get_unmapped_value( level ) ;
  			std::vector< std::size_t > const samples_in_strata = samples.find_samples_by_value( column, unmapped_value ) ;
  			for( std::size_t i = 0; i < samples_in_strata.size(); ++i ) {
  				result.add_sample( column + "=" + to_string( unmapped_value ), samples_in_strata[i] ) ;
  			}
  		}
  		return result ;
  	}
  	
  	genfile::SNPDataSource::UniquePtr open_genotype_data_sources(
  		std::string const& filename,
  		genfile::CommonSNPFilter::UniquePtr filter,
  		std::set< std::size_t > const& excluded_samples,
  		bool threshold = false
  	) {
  		std::vector< genfile::wildcard::FilenameMatch > filenames
  			= genfile::wildcard::find_files_by_chromosome(
  				filename,
  				genfile::wildcard::eALL_CHROMOSOMES
  			) ;
  
  		genfile::SNPDataSource::UniquePtr source( 
  			genfile::SNPDataSourceChain::create(
  				filenames
  			).release()
  		) ;
  		if( filter.get() ) {
  			source.reset(
  				genfile::VariantIdentifyingDataFilteringSNPDataSource::create(
  					source,
  					genfile::VariantIdentifyingDataTest::UniquePtr( filter.release() )
  				).release()
  			) ;
  		}
  		if( threshold ) {
  			source.reset(
  				new genfile::ThreshholdingSNPDataSource( source, 0.9 )
  			) ;
  		}
  		if( excluded_samples.size() > 0 ) {
  			source.reset(
  				new genfile::SampleFilteringSNPDataSource( source, excluded_samples )
  			) ;
  		}
  		return source ;
  	}
  
  	void write_preamble(
  		genfile::SNPDataSource const& host,
  		genfile::SNPDataSource const& para,
  		genfile::CohortIndividualSource const& samples
  	) {
  		get_ui_context().logger() << "Loaded data for " << samples.size() << " samples.\n" ;
  		get_ui_context().logger() << "  Predictor data:\n" << host.get_summary() << "\n" ;
  		get_ui_context().logger() << "    Outcome data:\n" << para.get_summary() << "\n" ;
  	}
  
  	void test(
  		genfile::CartesianProductVisitor& visitor,
  		genfile::CohortIndividualSource const& samples,
  		qcdb::MultiVariantStorage& output
  	) {
  		using namespace metro ;
  		typedef std::vector< SampleRange > SampleRanges ;
  
  		std::vector< std::string > model_names ;
  		boost::ptr_vector< regression::Design > designs ;
  		std::vector< Eigen::MatrixXd > predictorCodings ;
  		build_models( samples, &model_names, &designs, &predictorCodings ) ;
  
  		// Reserve some space for our variants
  		std::vector< genfile::VariantIdentifyingData > variants( 2 ) ;
  //		genfile::VariantIdentifyingData& pv = variants[0] ;
  //		genfile::VariantIdentifyingData& ov = variants[1] ;
  
  		// Storage
  		Matrix outcome = Matrix::Zero( samples.size(), 2 ) ;
  		Vector outcome_ploidy = Vector::Zero( samples.size() ) ;
  		SampleRanges nonmissing_outcome ;
  
  		Matrix predictor_probabilities = Matrix::Zero( samples.size(), 3 ) ;
  		IntegerVector predictor_ploidy = IntegerVector::Zero( samples.size() ) ;
  		SampleRanges nonmissing_predictor ;
  
  		stats::impl::InfoComputation info_computation;
  
  		{
  			appcontext::UIContext::ProgressContext progress_context = get_ui_context().get_progress_context( "Testing" ) ;
  			std::size_t count = 0 ;
  			
  			bool have_summarised_models = false ;
  			bool outputAllVariants = options().check( "-output-all-variants" ) ;
  			bool const debug = options().check( "-debug" ) ;
  
  			while(
  				visitor.step(
  					[&] (
  						std::vector< int > const& changed,
  						std::vector< genfile::VariantIdentifyingData > const& variants,
  						std::vector< genfile::VariantDataReader::SharedPtr > const& readers
  					) {
  						bool predictorIncluded = true ;
  						double predictorCount = 0.0 ;
  						bool outcomeIncluded = true ;
  						double outcomeCount = 0.0 ;
  						if( changed[0] == 1 ) {
  //							std::cerr << "Predictor changed!" << ".\n" ;
  							std::pair< bool, double > const result = set_predictor(
  								variants[0],
  								readers[0],
  								&predictor_probabilities,
  								&predictor_ploidy,
  								&nonmissing_predictor,
  								predictorCodings,
  								&designs
  							) ;
  							predictorIncluded = result.first ;
  							predictorCount = result.second ;
  						}
  
  						if( changed[1] == 1 ) {
  //							std::cerr << "Outcome changed!" << ".\n" ;
  							std::pair< bool, double > const result = set_outcome(
  								variants[1],
  								readers[1],
  								&outcome,
  								&outcome_ploidy,
  								&nonmissing_outcome,
  								&designs
  							) ;
  							outcomeIncluded = result.first ;
  							outcomeCount = result.second ;
  						}
  				
  						if( !have_summarised_models ) {
  							summarise_models( model_names, designs ) ;
  							have_summarised_models = true ;
  						}
  					
  						if( outputAllVariants || (predictorIncluded && outcomeIncluded)) {
  							output_design( designs[1], variants, designs[1].nonmissing_samples(), output ) ;
  
  							output.store_data_for_key( variants, "minimum_outcome_count", outcomeCount ) ;
  							output.store_data_for_key( variants, "minimum_predictor_count", predictorCount ) ;
  						
  							std::vector< metro::SampleRange > const included_samples = metro::impl::intersect_ranges(
  								nonmissing_predictor,
  								nonmissing_outcome
  							) ;
  
  							if( debug ) {
  								std::cerr << "INCLUDED SAMPLES: " << included_samples << "\n" ;
  							}
  						
  							output_cross_counts(
  								variants,
  								output,
  								cross_tabulate(
  									predictor_probabilities,
  									outcome,
  									outcome_ploidy,
  									included_samples
  								)
  							) ;
  								
  							info_computation.compute(
  								variants[0],
  								predictor_probabilities,
  								predictor_ploidy,
  								included_samples
  							) ;
  							output.store_data_for_key( variants, "predictor_info", info_computation.info() ) ;
  						}
  						// Set up ll
  						if( predictorIncluded && outcomeIncluded ) {
  							test( model_names, designs, variants, output ) ;
  						} 
  					}
  				)
  			) {
  				progress_context( ++count, visitor.count() ) ;
  			}
  		}
  	}
  
  	std::pair< bool, double > set_predictor(
  		genfile::VariantIdentifyingData const& variant,
  		genfile::VariantDataReader::SharedPtr predictor_source,
  		Matrix* predictor_probabilities,
  		IntegerVector* predictor_ploidy,
  		std::vector< metro::SampleRange >* nonmissing_predictor,
  		std::vector< Eigen::MatrixXd > const& predictorCodings,
  		boost::ptr_vector< metro::regression::Design >* designs
  	) {
  		bool result = true ;
  		bool const debug = options().check( "-debug" ) ;
  		ProbSetter predictor_setter( predictor_probabilities, predictor_ploidy, nonmissing_predictor ) ;
  		if( options().check( "-treat-predictor-as-haploid" )) {
  			predictor_setter.set_coerce_to_haploid() ;
  		}
  		predictor_source->get( ":genotypes:", genfile::to_GP_unphased( predictor_setter )) ;
  
  		Eigen::MatrixXd stratified_predictor_probs ;
  
  		stratify_predictors(
  			*predictor_probabilities,
  			m_stratification,
  			&m_regression_predictor_probs
  		) ;
  		
  		// designs[0] is null model, specify a 0-column matrix.
  		(*designs)[0].set_predictors(
  			Eigen::MatrixXd::Zero( 1, 0 ),
  			m_regression_predictor_probs.rowwise().sum(),
  			*nonmissing_predictor
  		) ;
  		// designs[1..] are alternative models
  		for( std::size_t i = 1; i < designs->size(); ++i ) {
  			assert( predictorCodings[i].rows() == m_regression_predictor_probs.cols() ) ;
  
  			// optimisation: ignore columns with zero probability
  			int lower = 0;
  			int upper = m_regression_predictor_probs.cols() ;
  //			for( ; lower < upper && (*predictor_probabilities).col(lower).array().abs().maxCoeff() == 0; ++lower ) ;
  //			for( ; upper > lower && (*predictor_probabilities).col(upper).array().abs().maxCoeff() == 0; --upper ) ;
  //			assert( upper >= lower ) ;
  
  			(*designs)[i].set_predictors(
  				predictorCodings[i].block( lower, 0, upper - lower, predictorCodings[i].cols() ),
  				m_regression_predictor_probs.block(
  					0, lower,
  					m_regression_predictor_probs.rows(), upper - lower
  				),
  				*nonmissing_predictor
  			) ;
  		}
  
  		double const predictorCount = 
  			std::min(
  				(predictor_probabilities->col(1) + 2.0 * predictor_probabilities->col(2)).array().sum(),
  				(predictor_probabilities->col(1) + 2.0 * predictor_probabilities->col(0)).array().sum()
  			) ;
  		if( predictorCount < options().get< double >( "-minimum-predictor-count" )) {
  			result = false ;
  		}
  
  		if( debug ) {
  			std::cerr << "PREDICTOR PROBS:\n"
  				<< predictor_probabilities->block(
  					0, 0,
  					std::min( 30, int( predictor_probabilities->rows() ) ), predictor_probabilities->cols()
  				).transpose() << "\n" ;
  			std::cerr << "PREDICTOR INCLUDED SAMPLES: " << (*nonmissing_predictor) << "\n" ;
  			std::cerr << "PREDICTOR CODINGS:\n" ;
  			for( std::size_t i = 1; i < designs->size(); ++i ) {
  				std::cerr << predictorCodings[i] << ".\n" ;
  			}
  		}
  
  		return std::make_pair( result, predictorCount ) ;
  	}
  
  	std::pair< bool, double > set_outcome(
  		genfile::VariantIdentifyingData const& variant,
  		genfile::VariantDataReader::SharedPtr outcomeReader,
  		Matrix* outcome,
  		Vector* outcome_ploidy,
  		std::vector< metro::SampleRange >* nonmissing_outcome,
  		boost::ptr_vector< metro::regression::Design >* designs
  	) {
  		bool result = true ;
  		bool const debug = options().check( "-debug" ) ;
  		AlleleCounter counter( outcome, outcome_ploidy, nonmissing_outcome ) ;
  		if( options().check( "-treat-outcome-as-haploid" )) {
  			counter.set_coerce_to_haploid() ;
  		}
  		if( outcomeReader->supports( "GT" )) {
  			outcomeReader->get( "GT", counter ) ;
  		} else if( outcomeReader->supports( "GP" )) {
  			outcomeReader->get(
  				"GP",
  				genfile::GPThresholdingGTSetter(
  					counter,
  					options().get< double >( "-outcome-genotype-call-threshold" )
  				)
  			) ;
  		} else {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::test()",
  				"outcome_source",
  				"Source must support hard genotype calls (GT) or genotype probabilities (GP) field."
  			) ;
  		}
  		if( debug ) {
  			std::cerr << "OUTCOME:\n" << outcome->block(
  				0, 0,
  				std::min( 30, int( outcome->rows() ) ), outcome->cols()
  			).transpose() << ".\n" ;
  			std::cerr << "OUTCOME INCLUDED SAMPLES: " << (*nonmissing_outcome) << "\n" ;
  		}
  
  		std::string const& outcomeName = options().get< std::string >( "-outcome-name" ) ;
  
  		for( std::size_t i = 0; i < designs->size(); ++i ) {
  			(*designs)[i].set_outcome(
  				*outcome,
  				*nonmissing_outcome,
  				std::vector< std::string >({outcomeName + "=0", outcomeName + "=1"} )
  			) ;
  		}
  
  		double const outcomeCount =  std::min( outcome->col(0).sum(), outcome->col(1).sum() ) ;
  		if( outcomeCount < options().get< double >( "-minimum-outcome-count" )) {
  			result = false ;
  		}
  
  		return std::make_pair( result, outcomeCount ) ;
  	}
  
  	Matrix cross_tabulate(
  		Matrix const& predictor,
  		Matrix const& outcome,
  		Vector const& outcome_ploidy,
  		std::vector< metro::SampleRange > const& ranges
  	) {
  		typedef Eigen::Block< Matrix > MatrixBlock ;
  		typedef Eigen::Block< Matrix const > ConstMatrixBlock ;
  
  		Matrix result = Matrix::Zero( outcome.cols(), predictor.cols() ) ;
  		for( std::size_t i = 0; i < ranges.size(); ++i ) {
  			metro::SampleRange const& range = ranges[i] ;
  			ConstMatrixBlock outcome_block = outcome.block( range.begin(), 0, range.size(), outcome.cols() ) ;
  			ConstMatrixBlock predictor_block = predictor.block( range.begin(), 0, range.size(), predictor.cols() ) ;
  			for( int j = 0; j < outcome.cols(); ++j ) {
  				for( int k = 0; k < predictor.cols(); ++k ) {
  					result(j, k) += ( outcome_block.col( j ).transpose() * predictor_block.col(k) ) ;
  				}
  			}
  		}
  		return result ;
  	}
  
  	void output_cross_counts(
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		qcdb::MultiVariantStorage& output,
  		Matrix const& cross_counts
  	) {
  
  		{
  			Matrix conversion( 3, 2 ) ;
  			conversion << 	2, 0,
  							1, 1,
  							0, 2 ;
  			Matrix predictorCount = cross_counts * conversion ;
  			Matrix::Index row, col ;
  			double const minimum_expected_predictor_count = predictorCount.minCoeff( &row, &col ) ;
  			output.store_data_for_key( variants, "minimum_expected_predictor_allele_count", minimum_expected_predictor_count ) ;
  			output.store_data_for_key(
  				variants,
  				"minimum_expected_predictor_allele_count_genotype",
  				( boost::format( "outcome=%d/predictor=%d" ) % row % col ).str()
  			) ;
  		}
  		
  		if( options().check( "-output-counts" )) {
  			for( int outcome_i = 0; outcome_i < cross_counts.rows(); ++outcome_i ) {
  				for( int predictor_i = 0; predictor_i < cross_counts.cols(); ++predictor_i ) {
  					output.store_data_for_key(
  						variants,
  						( boost::format( "outcome=%d/predictor=%d" ) % outcome_i % predictor_i ).str(),
  						cross_counts( outcome_i, predictor_i )
  					) ;
  				}
  			}
  		}
  	}
  
  	void build_models(
  		genfile::CohortIndividualSource const& samples,
  		std::vector< std::string >* names,
  		boost::ptr_vector< metro::regression::Design >* designs,
  		std::vector< Eigen::MatrixXd >* predictorCodings
  	) {
  		build_unadjusted_models( samples, names, designs, predictorCodings ) ;
  		if( options().check( "-covariates" )) {
  			add_covariates( samples, designs, options().get_values< std::string >( "-covariates" ) ) ;
  		}
  	}
  
  	void build_unadjusted_models(
  		genfile::CohortIndividualSource const& samples,
  		std::vector< std::string >* names,
  		boost::ptr_vector< metro::regression::Design >* designs,
  		std::vector< Eigen::MatrixXd >* predictorCodings
  	) {
  		typedef std::vector< metro::SampleRange > SampleRanges ;
  		typedef std::vector< std::string > Names ;
  
  		std::string const& outcomeName = options().get< std::string >( "-outcome-name" ) ;
  
  		// baseline model, no predictor
  		{
  			names->push_back( "null" ) ;
  			designs->push_back(
  				metro::regression::Design::create(
  					Matrix::Zero( samples.size(), 2 ), SampleRanges(), Names({ outcomeName + "=0", outcomeName + "=1" }),
  					Names() // only baseline
  				)
  			) ;
  			predictorCodings->push_back( Eigen::MatrixXd(3,0) ) ;
  		}
  
  		// alternative models
  		std::vector< std::string > const& models = options().get_values< std::string >( "-model" ) ;
  		// Null model encoding
  
  		using genfile::string_utils::slice ;
  		// Alternative models encoding
  		for( std::size_t i = 0; i < models.size(); ++i ) {
  			std::vector< slice > elts = parse_model_spec( models[i] ) ;
  			std::string const& predictorCoding = elts[1] ;
  		
  			Matrix predictor_levels ;
  			Names predictorNames ;
  			if( predictorCoding == "add+het" ) {
  				predictor_levels = Eigen::MatrixXd::Zero( 3, 2 ) ;
  				// add predictor
  				predictor_levels(0,0) = 0 ;
  				predictor_levels(1,0) = 1 ;
  				predictor_levels(2,0) = 2 ;
  				// het predictor
  				predictor_levels(1,1) = 1 ;
  				predictorNames = Names({"add", "overdominance"}) ;
  			} else if( predictorCoding == "add" ) {
  				predictor_levels = Eigen::MatrixXd::Zero( 3, 1 ) ;
  				// add predictor
  				predictor_levels(0,0) = 0 ;
  				predictor_levels(1,0) = 1 ;
  				predictor_levels(2,0) = 2 ;
  				predictorNames = Names({"add"}) ;
  			} else if( predictorCoding == "dom" ) {
  				predictor_levels = Eigen::MatrixXd::Zero( 3, 1 ) ;
  				// add predictor
  				predictor_levels(0,0) = 0 ;
  				predictor_levels(1,0) = 1 ;
  				predictor_levels(2,0) = 1 ;
  				predictorNames = Names({"dom"}) ;
  			} else if( predictorCoding == "rec" ) {
  				predictor_levels = Eigen::MatrixXd::Zero( 3, 1 ) ;
  				// add predictor
  				predictor_levels(0,0) = 0 ;
  				predictor_levels(1,0) = 0 ;
  				predictor_levels(2,0) = 1 ;
  				predictorNames = Names({"rec"}) ;
  			} else if( predictorCoding == "het" ) {
  				predictor_levels = Eigen::MatrixXd::Zero( 3, 1 ) ;
  				// add predictor
  				predictor_levels(0,0) = 0 ;
  				predictor_levels(1,0) = 1 ;
  				predictor_levels(2,0) = 0 ;
  				predictorNames = Names({"het"}) ;
  			} else {
  				assert(0) ;
  			}
  
  			if( m_stratification.size() > 0 ) {
  				stratify_model( &predictor_levels, &predictorNames, m_stratification ) ;
  			}
  
  			names->push_back( elts[0] ) ;
  			designs->push_back(
  				metro::regression::Design::create(
  					Matrix::Zero( samples.size(), 2 ), SampleRanges(), Names({ "outcome=0", "outcome=1" }),
  					predictorNames
  				)
  			) ;
  			predictorCodings->push_back( predictor_levels ) ;
  		}
  	}
  	
  	// This must be implemented to match stratify_predictors() below.
  	// Each column of the levels corresponds to a parameter and each row to a predictor level.
  	// We expand the predictors so that each parameter is split into one parameter per strata.
  	void stratify_model(
  		Eigen::MatrixXd* predictor_levels,
  		std::vector< std::string >* predictor_names,
  		metro::SampleStratification const& stratification
  	) {
  		int const L = predictor_levels->rows() ;
  		int const P = predictor_levels->cols() ;
  		Eigen::MatrixXd stratified(
  			L * stratification.number_of_strata(),
  			P * stratification.number_of_strata()
  		) ;
  		stratified.setZero() ;
  		std::vector< std::string > stratified_names ;
  		for( std::size_t stratum_i = 0; stratum_i < stratification.number_of_strata(); ++stratum_i ) {
  			for( std::size_t j = 0; j < stratification.stratum(stratum_i).size(); ++j ) {
  				stratified.block(
  					stratum_i * L, stratum_i,
  					L, P
  				) = (*predictor_levels) ;
  			}
  			for( std::size_t name_i = 0; name_i < predictor_names->size(); ++name_i ) {
  				stratified_names.push_back( predictor_names->at(name_i) + "/" + stratification.stratum_name( stratum_i )) ;
  			}
  		}
  		
  		(*predictor_levels) = stratified ;
  		(*predictor_names) = stratified_names ;
  	}
  	
  	// This must be implemented to match stratify_models() above.
  	void stratify_predictors(
  		Eigen::MatrixXd const& predictor_probs,
  		metro::SampleStratification const& stratification,
  		Eigen::MatrixXd* result
  	) {
  		if( stratification.size() == 0 ) {
  			*result = predictor_probs ;
  			return ;
  		}
  		int const N = predictor_probs.rows() ;
  		int const L = predictor_probs.cols() ;
  
  		Eigen::MatrixXd stratified_predictor_probs( N, L * stratification.size() ) ;
  		stratified_predictor_probs.setZero() ;
  		for( std::size_t i = 0; i < stratification.size(); ++i ) {
  			for( std::size_t j = 0; j < stratification.stratum(i).size(); ++j ) {
  				metro::SampleRange const& range = stratification.stratum(i)[j] ;
  				Eigen::Block< Eigen::MatrixXd const > const block = predictor_probs.block(
  					range.begin(), 0,
  					range.size(), L
  				) ;
  				stratified_predictor_probs.block(
  					range.begin(), L * i,
  					range.size(), L
  				) = block ;
  			}
  		}
  		
  		*result = stratified_predictor_probs ;
  	}
  	
  	std::vector< genfile::string_utils::slice > parse_model_spec( genfile::string_utils::slice const& spec ) {
  		std::vector< genfile::string_utils::slice > result = spec.split( ":" ) ;
  		if( result.size() == 1 ) {
  			// model name is same as spec
  			result.insert( result.begin(), result[0] ) ;
  		}
  		if( result.size() != 2 ) {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::parse_model_spec()",
  				"spec=\"" + spec + "\"",
  				"Expected a model spec in the form <name>:<model>"
  			) ;
  		}
  		return result ;
  	}
  
  	void add_covariates(
  		genfile::CohortIndividualSource const& samples,
  		boost::ptr_vector< metro::regression::Design >* designs,
  		std::vector< std::string > const& covariates
  	) {
  		get_ui_context().logger() << "Adding covariates...\n" ;
  
  		genfile::CohortIndividualSource::ColumnSpec const& spec = samples.get_column_spec() ;
  		for( std::size_t i = 0; i < covariates.size(); ++i ) {
  			std::string const& covariateName = covariates[i] ;
  			genfile::CohortIndividualSource::SingleColumnSpec columnSpec = spec[covariateName] ;
  			genfile::CrossCohortCovariateValueMapping::UniquePtr mapping
  				= genfile::CrossCohortCovariateValueMapping::create( columnSpec, true ) ;
  			mapping->add_source( samples ) ;
  		
  			if( mapping->get_number_of_distinct_mapped_values() == 0 ) {
  				throw genfile::BadArgumentError(
  					"add_covariates()",
  					"-covarariates"
  					"Covariate \"" + covariates[i] + "\" has no non-missing levels in this set of samples."
  				) ;
  			}
  			for( std::size_t model_i = 0; model_i < designs->size(); ++model_i ) {
  				metro::regression::Design& design = designs->at(model_i) ;
  				
  				if( columnSpec.is_discrete() ) {
  					design.add_discrete_covariate(
  						covariateName,
  						boost::bind(
  							&extract_mapped_categorical_value,
  							_1,
  							boost::cref( samples ),
  							spec.find_column( covariateName ),
  							boost::cref( *mapping )
  						),
  						boost::bind(
  							&get_unmapped_level,
  							boost::cref( *mapping ),
  							_1
  						),
  						mapping->get_number_of_distinct_mapped_values()
  					) ;
  				} else {
  					design.add_single_covariate(
  						covariateName,
  						boost::bind(
  							&extract_mapped_continuous_value,
  							_1,
  							boost::cref( samples ),
  							spec.find_column( covariateName ),
  							boost::cref( *mapping )
  						)
  					) ;
  				}
  			}
  			
  			get_ui_context().logger()
  				<< "++ Added covariate: \"" + covariateName + "\":\n"
  				<<  mapping->get_summary( "     " )
  				<< "\n\n" ;
  		}
  	}
  	
  	void summarise_models(
  		std::vector< std::string > const& model_names,
  		boost::ptr_vector< metro::regression::Design >& designs
  	) {
  		get_ui_context().logger() << "Models are:\n" ;
  		for( std::size_t i = 0; i < designs.size(); ++i ) {
  			metro::regression::LogLikelihood::UniquePtr model = create_loglikelihood( designs[i] ) ; 
  			get_ui_context().logger() << "- Model " << (i+1) << " (\"" + model_names[i] << "\"): "
  				<< model->get_summary() << "\n" ;
  		}
  		get_ui_context().logger() << "Model design for first test:\n" << designs.back().get_summary() << ".\n" ;
  		get_ui_context().logger() << "\n" ;
  	}
  	
  	void output_design(
  		metro::regression::Design& design,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		std::vector< metro::SampleRange > const& included_samples,
  		qcdb::MultiVariantStorage& output
  	) {
  		metro::regression::Design::Matrix const& outcome = design.outcome() ;
  #if DEBUG
  		//std::cerr << "OUTCOME:\n" << outcome << "\n" ;
  #endif
  		std::map< int, int > outcomeCounts ;
  		int const maxOutcome = std::max( outcome.array().maxCoeff(), 2.0 ) ;
  		for( int i = 0; i <= maxOutcome; ++i ) {
  			outcomeCounts[i] = 0 ;
  		}
  		// Count outcome levels, sensitive to the included samples
  		int total = 0 ;
  		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
  			total += (included_samples[i].end() - included_samples[i].begin()) ;
  			for( int j = included_samples[i].begin(); j < included_samples[i].end(); ++j ) {
  				++outcomeCounts[outcome(j,1)] ;
  			}
  		}
  
  		output.store_data_for_key( variants, "N", total ) ;
  		output.store_data_for_key( variants, "missing", int( design.outcome().rows() - total ) ) ;
  		
  		for(
  			std::map< int, int >::const_iterator i = outcomeCounts.begin();
  			i != outcomeCounts.end();
  			++i
  		) {
  			output.store_data_for_key(
  				variants,
  				( boost::format( "outcome=%d" ) % i->first ).str(),
  				i->second
  			) ;
  		}
  		
  		for( int g = 0; g < design.get_number_of_predictor_levels(); ++g ) {
  			double sum = 0.0 ;
  			for( std::size_t i = 0; i < included_samples.size(); ++i ) {
  				sum += design.get_predictor_level_probabilities().col(g).segment(
  					included_samples[i].begin(),
  					included_samples[i].end() - included_samples[i].begin() 
  				).sum() ;
  			}
  			
  			output.store_data_for_key(
  				variants,
  				( boost::format( "predictor=%d" ) % g ).str(),
  				sum
  			) ;
  		}
  	}
  
  	metro::regression::LogLikelihood::UniquePtr create_loglikelihood( metro::regression::Design& design ) {
  		uint32_t const threads = options().get< uint32_t >( "-threads" ) ;
  		metro::regression::LogLikelihood::UniquePtr ll ;
  		if( threads > 0 ) {
  			if( !m_pool.get() ) {
  				m_pool = metro::concurrency::threadpool::create( threads ) ;
  			}
  			ll.reset(
  				metro::regression::ThreadedLogLikelihood::create(
  					design,
  					[]( metro::regression::Design& design, std::vector< metro::SampleRange > range ) {
  						metro::regression::LogLikelihood::UniquePtr ll(
  							metro::regression::BinomialLogistic::create( design, range ).release()
  						) ;
  						return ll ;
  					},
  					*m_pool
  				).release()
  			) ;
  		} else {
  			ll.reset(
  				metro::regression::BinomialLogistic::create( design ).release()
  			) ;
  		}
  		if( !options().check_if_option_has_value( "-no-prior" )) {
  			ll = apply_priors( ll, m_prior_specs ) ;
  		}
  		return ll ;
  	}
  
  	std::vector< std::string > expand_prior_specs(
  		std::vector< std::string > const& specs
  	) const {
  		// Expand any specs by the stratification, if needed
  		std::string const& outcome_name = options().get< std::string >( "-outcome-name" ) ;
  		std::vector< std::string > result;
  		for( std::size_t i = 0; i < specs.size(); ++i ) {
  			std::string spec = specs[i] ;
  
  			// replace [outcome] with outcome name
  			{
  				std::size_t where = spec.find( "[outcome]" ) ;
  				if( where != std::string::npos ) {
  					spec = spec.substr( 0, where )
  						+ outcome_name
  						+ spec.substr( where + 9, spec.size() ) ;
  				}
  			}
  
  			std::size_t where = spec.find( "/[stratum]" ) ;
  			if( where == std::string::npos ) {
  				result.push_back( specs[i] ) ;
  			} else if( m_stratification.size() == 0 ) {
  				// no stratification, remove /[stratum] part of spec
  				result.push_back( 
  					spec.substr( 0, where )
  					+ spec.substr( where + 10, spec.size() )
  				) ;
  			} else {
  				for( std::size_t j = 0; j < m_stratification.size(); ++j ) {
  					result.push_back( 
  					spec.substr( 0, where )
  						+ "/" + m_stratification.stratum_name(j)
  						+ spec.substr( where + 10, spec.size() )
  					) ;
  				}
  			}
  		}
  		return result ;	
  	}
  
  	metro::regression::LogLikelihood::UniquePtr apply_priors(
  		metro::regression::LogLikelihood::UniquePtr ll,
  		std::vector< std::string > specs
  	) const {
  		std::vector< int > gaussian_parameter_indices ;
  		std::vector< double > gaussian_variances ;
  		std::vector< double > gaussian_means ;
  		std::vector< int > logF_parameter_indices ;
  		std::vector< double > logF_nu1 ;
  		std::vector< double > logF_nu2 ;
  		std::map< std::string, int > used_parameters ;
  
  		using genfile::string_utils::to_lower ;
  		using genfile::string_utils::to_repr ;
  
  		typedef std::map< std::string, int > ParameterNameMap ;
  		ParameterNameMap parameter_index_by_name ;
  		{
  			metro::regression::LogLikelihood::IntegerMatrix const identity = ll->identify_parameters() ;
  			for( int i = 0; i < identity.rows(); ++i ) {
  				parameter_index_by_name[ ll->get_parameter_name(i) ] = i ;
  			}
  		}
  	
  		// Gather and parse all the data
  		for( std::size_t i = 0; i < specs.size(); ++i ) {
  			//std::cerr << "SPEC: \"" << specs[i] << "\".\n" ;
  			std::string parameter_name ;
  			std::string distribution ;
  			std::vector< std::string > parameters ;
  			parse_prior_spec(
  				specs[i],
  				&parameter_name,
  				&distribution,
  				&parameters,
  				options().get< std::string >( "-outcome-name" )
  			) ;
  
  			++used_parameters[parameter_name] ;
  			if( used_parameters[parameter_name] > 1 ) {
  				throw genfile::BadArgumentError(
  					"PerVariantComputationManager::apply_priors()",
  					"spec=\"" + specs[i] + "\"",
  					(boost::format(
  						"More than one prior specified for parameter \"%s\"" ) % parameter_name ).str()
  				) ;
  			}
  			
  			ParameterNameMap::const_iterator where = parameter_index_by_name.find( parameter_name ) ;
  			if( where != parameter_index_by_name.end() ) {
  				if( to_lower( distribution ) == "gaussian" ) {
  					if( parameters.size() != 2 ) {
  						throw genfile::BadArgumentError(
  							"PerVariantComputationManager::apply_priors()",
  							"spec=\"" + specs[i] + "\"",
  							(boost::format(
  								"Wrong number of parameters found for gaussian distribution (\"%d\", expected 2) - please specify mean and variance."
  									) % parameters.size()).str()
  						) ;
  					}
  				
  					gaussian_means.push_back( to_repr< double >( parameters[0] )) ;
  					gaussian_variances.push_back( to_repr< double >( parameters[1] )) ;
  					gaussian_parameter_indices.push_back( where->second ) ;
  				} else if( to_lower( distribution ) == "logf" ) {
  					if( parameters.size() != 2 ) {
  						throw genfile::BadArgumentError(
  							"PerVariantComputationManager::apply_priors()",
  							"spec=\"" + specs[i] + "\"",
  							(boost::format(
  								"Wrong number of parameters found for log F distribution (\"%d\", expected 2) - please specify nu1 and nu2"
  									) % parameters.size()).str()
  						) ;
  					}
  					logF_nu1.push_back( to_repr< double >( parameters[0] )) ;
  					logF_nu2.push_back( to_repr< double >( parameters[1] )) ;
  					logF_parameter_indices.push_back( where->second ) ;
  				} else {
  					throw genfile::BadArgumentError(
  						"PerVariantComputationManager::apply_priors()",
  						"spec=\"" + specs[i] + "\"",
  						(boost::format( "Unrecognised distribution specified (\"%s\"), should be \"gaussian\" or \"logF\"" ) % distribution).str()
  					) ;
  				}
  			}
  		}
  		
  		if( !options().check( "-no-covariate-priors" )) {
  			metro::regression::LogLikelihood::IntegerMatrix const identity = ll->identify_parameters() ;
  			for( int i = 0; i < identity.rows(); ++i ) {
  				std::string const& parameter_name = ll->get_parameter_name(i) ;
  				if( identity(i,1) > ll->design().number_of_predictors() && used_parameters[parameter_name] == 0 ) {
  					logF_parameter_indices.push_back(i) ;
  					// Choose a logF distribution with 95% coverage interval
  					// equivalent to a gaussian(0,40^2) distribution.
  					logF_nu1.push_back( 0.07704070 ) ;
  					logF_nu2.push_back( 0.07704070 ) ;
  				}
  			}
  		}
  		
  		
  		// Now apply the priors:
  		metro::regression::LogLikelihood::UniquePtr result( ll.release() ) ;
  		if( gaussian_parameter_indices.size() > 0 ) {
  			result.reset(
  				metro::regression::IndependentNormalWeightedLogLikelihood::create(
  					result,
  					gaussian_parameter_indices,
  					gaussian_means, gaussian_variances
  				).release()
  			) ;
  		}
  		if( logF_parameter_indices.size() > 0 ) {
  			result.reset(
  				metro::regression::IndependentLogFWeightedLogLikelihood::create(
  					result,
  					logF_parameter_indices,
  					logF_nu1, logF_nu2
  				).release()
  			) ;
  		}
  		return( result ) ;
  	}
  
  	// parse a parameter prior distribution specification string
  	// The spec should be of the form:
  	// <parameter_name>~<family>(<parameters>)
  	// Currently acceptable families are:
  	// logf (for a logF distribution)
  	// gaussian (for a gaussian distribution)
  	void parse_prior_spec(
  		std::string const& spec,
  		std::string* parameter_name,
  		std::string* distribution,
  		std::vector< std::string >* parameters,
  		std::string const& outcome_name
  	) const {
  		using genfile::string_utils::split ;
  		using genfile::string_utils::to_repr ;
  		using genfile::string_utils::slice ;
  		std::vector< slice > const elts = slice( spec ).split( "~" ) ;
  		if( elts.size() != 2 ) {
  			throw genfile::BadArgumentError(
  				"PerVariantComputationManager::parse_prior_spec()",
  				"spec=\"" + spec + "\"",
  				(boost::format(
  					"Malformatted value (\"%s\") specified, should be of the form <parameter name>~<family>(<parameters>)\n"
  						" where family is one of \"gaussian\" or \"logF\"") % spec).str()
  			) ;
  		}
  		std::vector< slice > const distribution_spec = elts[1].split( "()" ) ;
  		if(
  			distribution_spec.size() != 3
  			|| distribution_spec[1].size() < 2
  			|| distribution_spec[2].size() > 0
  		) {
  			throw genfile::BadArgumentError(
  				"PerVariantComputationManager::parse_prior_spec()",
  				"spec=\"" + spec + "\"",
  				(boost::format(
  					"Malformatted distribution (\"%s\") specified, should be of the form <family>(<parameters>)\n"
  						" where family is one of \"gaussian\" or \"logF\"") % elts[1]).str()
  			) ;
  		}
  		std::vector< slice > const parameter_spec = distribution_spec[1].split( "," ) ;
  
  		{
  			std::string parameterSpec = elts[0] ;
  			std::size_t where = parameterSpec.find( "[outcome]" ) ;
  			if( where != std::string::npos ) {
  				parameterSpec = parameterSpec.substr( 0, where )
  					+ outcome_name
  					+ parameterSpec.substr( where + 9, parameterSpec.size() ) ;
  			}
  			*parameter_name = parameterSpec ;
  		}
  
  		{
  			std::string parameterSpec = elts[0] ;
  			std::size_t where = parameterSpec.find( "/[stratum]" ) ;
  			if( where != std::string::npos ) {
  				if( m_stratification.size() == 0 ) {
  					parameterSpec = parameterSpec.substr( 0, where )
  						+ ""
  						+ parameterSpec.substr( where + 9, parameterSpec.size() ) ;
  				}
  			}
  			*parameter_name = parameterSpec ;
  		}
  		
  		*distribution = distribution_spec[0] ;
  		parameters->assign( parameter_spec.begin(), parameter_spec.end() ) ;
  	}
  
  	void test(
  		std::vector< std::string > const& model_names,
  		boost::ptr_vector< metro::regression::Design >& designs,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		qcdb::MultiVariantStorage& output
  	) {
  		using namespace metro::regression ;
  		boost::ptr_vector< LogLikelihood > lls ;
  		for( std::size_t i = 0; i < designs.size(); ++i ) {
  			lls.push_back( create_loglikelihood( designs[i] )) ;
  		}
  
  		std::string const output_models = options().get< std::string > ( "-output-models" ) ;
  		if( output_models != "all" && output_models != "alternative" ) {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::test()",
  				"output_models=\"" + output_models + "\"",
  				"Please specify \"alternative\" or \"all\"."
  			) ;
  		}
  		std::vector< std::string > comments ;
  		for( std::size_t model = 0; model < lls.size(); ++model ) {
  			boost::timer::cpu_timer timer ;
  			LogLikelihood& ll = lls[model] ;
  			std::string const& model_name = model_names[model] ;
  
  			metro::CholeskyStepper::Tracer tracer ;
  			if( options().check( "-debug" )) {
  				tracer = [this] ( 
  					int iteration,
  					double ll,
  					double target_ll,
  					metro::CholeskyStepper::Vector const& point,
  					metro::CholeskyStepper::Vector const& first_derivative,
  					metro::CholeskyStepper::Vector const& step,
  					bool converged
  				) {
  					this->get_ui_context().logger()
  						<< " ITERATION:" << iteration << "\n"
  						<< "        LL:" << ll << "\n"
  						<< "    TARGET:" << target_ll << "\n"
  						<< "     POINT:" << point.transpose() << "\n"
  						<< "DERIVATIVE:" << first_derivative.transpose() << "\n"
  						<< "      STEP:" << step.transpose() << "\n"
  						<< " CONVERGED:" << ( converged ? "yes" : "no" ) << "\n" ;
  				} ;
  			}
  
  			metro::CholeskyStepper stopping_condition(
  					options().get_value< double >( "-tolerance" ),
  					options().get_value< double >( "-max-iterations" ),
  					tracer
  				) ;
  
  			std::pair< bool, int > result = metro::fit_model(
  				ll,
  				model_name,
  				Eigen::VectorXd::Zero( ll.parameters().size() ),
  				stopping_condition,
  				&comments
  			) ;
  
  #if DEBUG
  			std::cerr << "hptest::test(): fit model in " << result.second << " iterations with " << ( result.first ? "convergence" : "no convergence" ) << ".\n" ;
  			for( std::size_t i = 0; i < ll.identify_parameters().rows(); ++i ) {
  				std::cerr << "model " << "NAME OF PARAMETER " << i << ": \"" << ll.get_parameter_name(i) << "\".\n" ;
  			}
  #endif
  		
  			output.store_data_for_key( variants, model_name + ":converged", result.first ? "1" : "0" ) ;
  			output.store_data_for_key( variants, model_name + ":iterations", result.second ) ;
  			output.store_data_for_key( variants, model_name + ":fit_time", boost::timer::format( timer.elapsed(), 4, "%t" ) ) ;
  			output.store_data_for_key( variants, model_name + ":ll", ll.get_value_of_function() ) ;
  			output.store_data_for_key(
  				variants,
  				model_name + ":degrees_of_freedom",
  				int( ll.parameters().size() - lls[0].parameters().size() )
  			) ;
  		
  			if( ( model > 0 || output_models == "all" ) && result.first ) {
  				output_results( lls[0], ll, variants, output, model_name, &comments ) ;
  			}
  		}
  		
  		if( comments.size() > 0 ) {
  			output.store_data_for_key(
  				variants,
  				"comment",
  				genfile::string_utils::join( comments, "," )
  			) ;
  		} else {
  			output.store_data_for_key(
  				variants,
  				"comment",
  				genfile::MissingValue()
  			) ;
  		}
  	}
  
  	void output_results(
  		metro::regression::LogLikelihood const& null_ll,
  		metro::regression::LogLikelihood const& ll,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		qcdb::MultiVariantStorage& output,
  		std::string const& model_name,
  		std::vector< std::string >* comments
  	) {
  		Eigen::ColPivHouseholderQR< Matrix > null_solver ;
  		Eigen::ColPivHouseholderQR< Matrix > solver ;
  		Matrix const& information = -ll.get_value_of_second_derivative() ;
  		null_solver.compute( -null_ll.get_value_of_second_derivative() ) ;
  		solver.compute( information ) ;
  		Matrix const covariance = solver.solve( Matrix::Identity( information.rows(), information.cols() ) ) ;
  #if DEBUG
  		std::cerr << "Model " << model_name << ": information:\n"
  			<< information
  				<< ".\n" ;
  #endif
  		if( check_model_fit( information, covariance, model_name, comments ) ) {
  			std::vector< int > const parametersToOutput = get_parameters_to_output( ll ) ;
  			output_parameter_estimates( ll, variants, model_name, covariance, output, parametersToOutput ) ;
  		
  			metro::regression::LogPosteriorDensity const* posterior = dynamic_cast< metro::regression::LogPosteriorDensity const* >( &ll ) ;
  			if( posterior ) {	
  				compute_bayesian_summary(
  					null_ll, null_solver.logAbsDeterminant(),
  					*posterior, solver.logAbsDeterminant(), covariance,
  					model_name,
  					variants, output,
  					parametersToOutput,
  					comments
  				) ;
  			} else {
  				compute_frequentist_summary(
  					null_ll, ll,
  					model_name,
  					variants, output,
  					comments
  				) ;
  			}
  		}
  	}
  
  	bool check_model_fit(
  		Matrix const& information,
  		Matrix const& covariance,
  		std::string const& model_name,
  		std::vector< std::string >* comments
  	) const {
  		double const relative_error = (( information * covariance ) - Matrix::Identity( information.rows(), information.cols() ) ).norm() ;
  		if( relative_error != relative_error || relative_error > options().get< double >( "-tolerance" ) ) {
  			if( comments ) {
  				comments->push_back( model_name + ":singular" ) ;
  			}
  			return false ;
  		} else {
  			return true ;
  		}
  	}
  	
  	std::vector< int > get_parameters_to_output( metro::regression::LogLikelihood const& ll ) {
  		std::vector< int > result ;
  		int lower = 0 ;
  		int upper = 0 ;
  		std::string const output_parameters = options().get< std::string >( "-output-parameters" ) ;
  		if( options().get< std::string >( "-output-parameters" ) == "all" ) {
  			lower = 0 ;
  			upper = ll.parameters().size() ;
  		} else if( output_parameters == "genetic" ) {
  			lower = 1 ;
  			upper = lower + ll.design().number_of_predictors() ;
  		} else {
  			throw genfile::BadArgumentError(
  				"HPTestApplication::output_parameter_estimates()",
  				"-output-paramaters=\"" + output_parameters + "\"",
  				"Expected \"all\" or \"genetic\"."
  			) ;
  		}
  		for( int i = lower; i < upper; ++i ) {
  			result.push_back(i) ;
  		}
  		return result ;
  	}
  	
  	void output_parameter_estimates(
  		metro::regression::LogLikelihood const& ll,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		std::string const& model_name,
  		Matrix const& covariance,
  		qcdb::MultiVariantStorage& output,
  		std::vector< int > const parametersToOutput
  	) const {
  		// Output parameter estimates, standard errors, etc.
  		metro::regression::LogLikelihood::IntegerMatrix const parameter_identity = ll.identify_parameters() ;
  		using genfile::string_utils::to_string ;
  
  		Vector const& parameters = ll.parameters() ;
  
  		for( std::size_t parameter_i = 0; parameter_i < parametersToOutput.size(); ++parameter_i ) {
  			int const i = parametersToOutput[parameter_i] ;
  			std::string const parameter_name = ll.get_parameter_name(i) ;
  			output.store_data_for_key( variants, model_name + ":beta_" + to_string( i ) + ":" + parameter_name, parameters( i ) ) ;
  		}
  
  		for( std::size_t parameter_i = 0; parameter_i < parametersToOutput.size(); ++parameter_i ) {
  			int const i = parametersToOutput[parameter_i] ;
  			output.store_data_for_key( variants, model_name + ":sd_" + to_string( i ) /*+ ":" + parameter_naming_scheme( i ) */, std::sqrt( covariance( i, i ) ) ) ;
  		}
  		for( std::size_t parameter_i = 0; parameter_i < parametersToOutput.size(); ++parameter_i ) {
  			int const i = parametersToOutput[parameter_i] ;
  			for( std::size_t parameter_j = parameter_i+1; parameter_j < parametersToOutput.size(); ++parameter_j ) {
  				int const j = parametersToOutput[parameter_j] ;
  				output.store_data_for_key( variants, model_name + ":cov_" + to_string( i ) + "," + to_string( j ), covariance( i, j ) ) ;
  			}
  		}
  	}
  	
  	double compute_twosided_wald_pvalue( double const beta, double const se ) const {
  		double result = std::numeric_limits< double >::quiet_NaN() ;
  		if( beta == beta && se > 0 ) {
  			boost::math::normal N( 0, se ) ;
  			result = 2 * boost::math::cdf( N, -std::abs( beta ) ) ;
  		}
  		return result ;
  	}
  	
  	void compute_bayesian_summary(
  		metro::regression::LogLikelihood const& null_ll,
  		double const logAbsNullDeterminant,
  		metro::regression::LogPosteriorDensity const& ll,
  		double const logAbsFullDeterminant,
  		Eigen::MatrixXd const& posterior_covariance,
  		std::string const& model_name,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		qcdb::MultiVariantStorage& output,
  		std::vector< int > const parametersToOutput,
  		std::vector< std::string >* comments
  	) const {
  		//  P(alt|D)       P(D|alt)    P(alt)           P(alt)
  		// ---------   =   -------- x  ------    = BF x ------ 
  		// P(null|D)       P(D|null)   P(null)          P(null)
  		//
  		// where
  		//      ∫ P(theta|alt)P(D|theta)
  		// BF = -------------------------
  		//      ∫ P(theta|null)P(D|theta)
  		//
  		// We have only un-normalised estimates of the top and bottom integrand.
  		
  		// Use a laplace approximation to compute numerator and denominator of the BF
  		// The point is that we have an (approximate) estimate of the unnormalised posterior
  		// under the null and alternative models.  The normalising constant gives us the
  		// posterior.
  		Vector const& alternative_parameters = ll.parameters() ;
  		Vector const& null_parameters = null_ll.parameters() ;
  		double const log_2pi = std::log( 2.0 * 3.1415926535897932384626 ) ;
  		double const log_numerator
  			= ll.get_value_of_function()
  				+ 0.5 * ( (alternative_parameters.size() * log_2pi) - logAbsFullDeterminant ) ;
  		double const log_denominator
  			= null_ll.get_value_of_function()
  				+ 0.5 * ( (null_parameters.size() * log_2pi) - logAbsNullDeterminant ) ;
  		double const log10_bf = ( log_numerator - log_denominator ) / std::log(10) ;
  
  		output.store_data_for_key(
  			variants,
  			model_name + ":log10_bf",
  			log10_bf
  		) ;
  		
  		metro::regression::LogPosteriorDensity const* posterior = &ll ;
  		assert( posterior ) ;
  		if( posterior ) {
  			Eigen::VectorXd const prior_mode = posterior->get_prior_mode() ;
  			Matrix const ll_information = -posterior->get_loglikelihood_second_derivative() ;
  			Matrix const approxSamplingCovariance = posterior_covariance * ll_information * posterior_covariance ;
  			Eigen::ColPivHouseholderQR< Matrix > solver( approxSamplingCovariance ) ;
  			if( solver.info() == Eigen::Success ) {
  	#if 0
  				std::cerr << "      Posterior covariance:\n" << posterior_covariance << ",\n" ;
  				std::cerr << "Approx sampling covariance:\n" << approxSamplingCovariance << ".\n" ;
  				std::cerr << "                  Approx Z:\n" << z << ".\n" ;
  	#endif		
  				using genfile::string_utils::to_string ;
  				for( std::size_t parameter_i = 0; parameter_i < parametersToOutput.size(); ++parameter_i ) {
  					int const i = parametersToOutput[parameter_i] ;
  					output.store_data_for_key(
  						variants,
  						model_name + ":prior_mode_" + to_string(i),
  						prior_mode(i)
  					) ;
  					double const sampling_se = std::sqrt( approxSamplingCovariance(i,i)) ;
  					double const posterior_sd = std::sqrt( posterior_covariance(i,i)) ;
  					output.store_data_for_key(
  						variants,
  						model_name + ":se_" + to_string(i),
  						sampling_se
  					) ;
  			
  					std::string const parameter_name = ll.get_parameter_name(i) ;
  					if( sampling_se < posterior_sd / 10 ) {
  						output.store_data_for_key(
  							variants,
  							model_name + ":pvalue_" + to_string( i ),
  							metro::NA
  						) ;
  						comments->push_back(
  							parameter_name + ":sampling_se_<<_posterior_sd:_p-value_may_be_unreliable"
  						) ;
  					} else {
  						output.store_data_for_key(
  							variants,
  							model_name + ":pvalue_" + to_string( i ),
  							compute_twosided_wald_pvalue( alternative_parameters(i), std::sqrt( approxSamplingCovariance(i,i) ) )
  						) ;
  					}
  				}
  			}
  		}
  	}
  
  	void compute_frequentist_summary(
  		metro::regression::LogLikelihood const& null_ll,
  		metro::regression::LogLikelihood const& ll,
  		std::string const& model_name,
  		std::vector< genfile::VariantIdentifyingData > const& variants,
  		qcdb::MultiVariantStorage& output,
  		std::vector< std::string >* comments
  	) const {
  		Vector const& parameters = ll.parameters() ;
  		Vector const& null_parameters = null_ll.parameters() ;
  		// Now output the test statistic and p-value
  		double p_value = compute_lrt_pvalue(
  			null_ll.get_value_of_function(), ll.get_value_of_function(),
  			parameters.size() - null_parameters.size(),
  			options().get< double >( "-tolerance" )
  		) ;
  		if( p_value != p_value || p_value < 0 ) {
  			if( comments ) {
  				// This should not happen because the convergence criterion checks for this.
  				comments->push_back( model_name + ":numerical_error" ) ;
  			}
  			p_value = std::numeric_limits< double >::quiet_NaN() ;
  		}
  		output.store_data_for_key(
  			variants,
  			model_name + ":lrt_pvalue",
  			p_value
  		) ;
  	}	
  	
  	double compute_lrt_pvalue(
  		double const null_ll,
  		double const alternative_ll,
  		int degrees_of_freedom,
  		double const tolerance
  	) const {
  		double p_value = std::numeric_limits< double >::quiet_NaN() ;
  		double test_statistic = -2.0 * ( null_ll - alternative_ll ) ;
  		if( degrees_of_freedom > 0 && test_statistic >= -tolerance ) {
  			test_statistic = std::max( 0.0, test_statistic ) ;
  			boost::math::chi_squared_distribution< double > chi_squared_distribution( degrees_of_freedom ) ;
  			p_value = boost::math::cdf(
  				boost::math::complement(
  					chi_squared_distribution,
  					test_statistic
  				)
  			) ;
  		}
  		return p_value ;
  	}
  	
  } ;
  
  int main( int argc, char** argv ) {
      try {
  		HPTestApplication app( argc, argv ) ;
      }
  	catch( appcontext::HaltProgramWithReturnCode const& e ) {
  		return e.return_code() ;
  	}
  	return 0 ;
  }
  
  