text:
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
// Name        :
// Author      : Avi
// Revision    : $Revision: #10 $
//
// Copyright 2009- ECMWF.
// This software is licensed under the terms of the Apache Licence version 2.0
// which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
// In applying this licence, ECMWF does not waive the privileges and immunities
// granted to it by virtue of its status as an intergovernmental organisation
// nor does it submit to any jurisdiction.
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
#include <iostream>

#include "boost/filesystem/operations.hpp"
#include <boost/test/unit_test.hpp>

#include "Defs.hpp"
#include "Suite.hpp"
#include "Family.hpp"
#include "Task.hpp"
#include "EcfFile.hpp"
#include "Str.hpp"
#include "File.hpp"
#include "Pid.hpp"

using namespace std;
using namespace ecf;
namespace fs = boost::filesystem;

BOOST_AUTO_TEST_SUITE( NodeTestSuite )

void located_ecf_file(task_ptr task, EcfFile::Origin expected_origin, EcfFile::EcfFileSearchAlgorithm expected_search, int line) {
   try {
      task->update_generated_variables();
      EcfFile ecf_file = task->locatedEcfFile();
      //cout << "located_ecf_file task: " << task->absNodePath() << " " << ecf_file.ecf_file_origin_dump() << "\n";
      BOOST_REQUIRE_MESSAGE( ecf_file.valid(), "\nline:" << line << " Expected search to pass with " << EcfFile::search_algorithm_str(expected_search));
      BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == expected_origin, "\nline:" << line  << " Expected ecf_file to be obtained from(" << EcfFile::origin_str(expected_origin) << ") but found " << ecf_file.ecf_file_origin_dump());
      BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_search_algorithm() == expected_search, "\nline:" << line << " Expected search " << EcfFile::search_algorithm_str(expected_search) << " but found " << EcfFile::search_algorithm_str(ecf_file.ecf_file_search_algorithm()));
   }
   catch (std::exception& e) {
      BOOST_REQUIRE_MESSAGE(true,"Expected search to succeed " << e.what());
   }
}

void located_ecf_file_fail(task_ptr task, int line) {
   task->update_generated_variables();
   try {
      EcfFile ecf_file = task->locatedEcfFile();
      BOOST_REQUIRE_MESSAGE( !ecf_file.valid(), "\nline:" << line << " Expected search to fail with standard look up(i.e PRUNE_ROOT) "<< ecf_file.ecf_file_origin_dump() );
   }
   catch (std::exception& e) {
      BOOST_REQUIRE_MESSAGE(true,"\nline:" << line << "Expected search to fail " << e.what());
   }
}

void create_ecf_file(const std::string& ecf_file_location) {
   //cout << "create file  = " << ecf_file_location << "\n";
   // generate the ecf file;
   string header = "%include <simple_head.h>\n";
   string body = "#body\n";
   string tail = "%include <simple_tail.h>\n";
   string ecf_file = header;
   ecf_file += body;
   ecf_file += tail;

   BOOST_CHECK_MESSAGE(File::createMissingDirectories(ecf_file_location),"Could not create missing dir\n");
   string errormsg;
   BOOST_CHECK_MESSAGE(File::create(ecf_file_location, ecf_file , errormsg), errormsg);
   BOOST_CHECK_MESSAGE(fs::exists(ecf_file_location), "Expected File " << ecf_file_location << " to exist");
}

BOOST_AUTO_TEST_CASE(  test_ecf_file_search )
{
   cout << "ANode:: ...test_ecf_file_search\n";
   // Do test which fails first, i.e does not find the ecf file, then add variable ECF_FILES_LOOKUP(PRUNE_LEAF)
   // The test should then succeed
   // Note: Setting ECF_FILES_LOOKUP(PRUNE_LEAF) or ECF_FILES_LOOKUP(PRUNE_ROOT)/default affects the lookup for
   //       ecf files in BOTH ECF_HOME and ECF_LISTS
   // PRUNE_ROOT/defaults
   //    <root>/suite/f1/f2/f3/task.ecf
   //    <root>/f1/f2/f3/task.ecf
   //    <root>/f2/f3/task.ecf
   //    <root>/f3/task.ecf
   //    <root>/task.ecf
   // PRUNE_LEAF( new from ecflow 4.12.0)
   //    <root>/suite/f1/f2/f3/task.ecf
   //    <root>/suite/f1/f2/task.ecf
   //    <root>/suite/f1/task.ecf
   //    <root>/suite/task.ecf
   //    <root>/task.ecf
   // IMPORTANT: The first and last lookups are the *SAME* for both search lookups.

   // Create a defs file corresponding to:
   //suite test_ecf_file_search
   // edit ECF_INCLUDE $ECF_HOME/includes
   // family f1
   //    family f2
   //       family f3
   //          task task
   Defs theDefs;
   suite_ptr suite = theDefs.add_suite(Pid::unique_name("test_ecf_file_search"));
   suite->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
   family_ptr f1 = suite->add_family( "f1" ) ;
   family_ptr f2 = f1->add_family( "f2" ) ;
   family_ptr f3 = f2->add_family( "f3" ) ;
   task_ptr task = f3->add_task( "task" ) ;

   // Override ECF_HOME. ECF_HOME is as default location for .ecf files, when ECF_INCLUDE not specified
   // or when file does not exist in ECF_INCLUDE
   std::string ecf_home = File::test_data("ANode/test/data","ANode");
   std::string ecf_lists = File::test_data("ANode/test/data","ANode");
   suite->add_variable(Str::ECF_HOME(),ecf_home);
   //cerr << theDefs << "\n";

   /// begin , will cause creation of generated variables. The generated variables
   /// are use in client scripts and used to locate the ecf files
   theDefs.beginAll();

   {
      //cout << "Test ECF_HOME/ECF_LISTS at leaf location this should be using ECF_SCRIPT\n";
      string ecf_file_location  = ecf_home + task->absNodePath() +  File::ECF_EXTN();
      create_ecf_file(ecf_file_location);

      located_ecf_file(task,EcfFile::ECF_SCRIPT,EcfFile::PRUNE_ROOT,__LINE__);
      suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );                 // change look up method
      located_ecf_file(task,EcfFile::ECF_SCRIPT,EcfFile::PRUNE_ROOT,__LINE__); // ECF_FILES_LOOKUP make *NO* difference for lookup with EcfFile::ECF_SCRIPT

      suite->deleteVariable("ECF_FILES_LOOKUP"); // cleanup

      suite->add_variable(Str::ECF_FILES(),ecf_lists);
      located_ecf_file(task,EcfFile::ECF_SCRIPT ,EcfFile::PRUNE_ROOT,__LINE__);// ECF_FILES make no difference for lookup with EcfFile::ECF_SCRIPT
      suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );                 // change look up method
      located_ecf_file(task,EcfFile::ECF_SCRIPT,EcfFile::PRUNE_ROOT,__LINE__); // ECF_FILES_LOOKUP make *NO* difference for lookup with EcfFile::ECF_SCRIPT

      suite->deleteVariable("ECF_FILES_LOOKUP");   // cleanup
      suite->deleteVariable("ECF_FILES");          // cleanup
      boost::filesystem::remove(ecf_file_location);// cleanup
   }
   {
      //cout << "Test ECF_HOME/ECF_LISTS at root location\n";
      string ecf_file_location  = ecf_home + "/" + task->name() +  File::ECF_EXTN();
      create_ecf_file(ecf_file_location);

      located_ecf_file(task,EcfFile::ECF_HOME,EcfFile::PRUNE_ROOT,__LINE__);
      suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );               // change look up method
      located_ecf_file(task,EcfFile::ECF_HOME,EcfFile::PRUNE_LEAF,__LINE__);

      suite->deleteVariable("ECF_FILES_LOOKUP"); // cleanup

      suite->add_variable(Str::ECF_FILES(),ecf_lists);
      located_ecf_file(task,EcfFile::ECF_FILES,EcfFile::PRUNE_ROOT,__LINE__);
      suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );                 // change look up method
      located_ecf_file(task,EcfFile::ECF_FILES,EcfFile::PRUNE_LEAF,__LINE__);

      suite->deleteVariable("ECF_FILES_LOOKUP");   // cleanup
      suite->deleteVariable("ECF_FILES");          // cleanup
      boost::filesystem::remove(ecf_file_location);// cleanup
   }
   {
      //cout << "Test ECF_HOME at intermediate location prune_leaf\n";
      Node* node = task.get();
      node = node->parent();
      node = node->parent();
      while(node) {

         string ecf_file_location = ecf_home + node->absNodePath() + "/" + task->name() + File::ECF_EXTN();
         create_ecf_file(ecf_file_location);

         // cout << "Node: " << node->absNodePath() << "\n";
         located_ecf_file_fail(task,__LINE__);

         // change look up method
         suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );
         located_ecf_file(task,EcfFile::ECF_HOME,EcfFile::PRUNE_LEAF,__LINE__);

         // Remove this file, required for the following test
         boost::filesystem::remove(ecf_file_location);

         // Remove ECF_FILES_LOOKUP variable
         suite->deleteVariable("ECF_FILES_LOOKUP");

         node = node->parent();
      }
      /// Remove all the generated files
      boost::filesystem::remove_all( ecf_home + suite->absNodePath() );
   }
   {
      //cout << "Test ECF_FILES at intermediate location, by prune_leaf\n";
      suite->add_variable(Str::ECF_FILES(),ecf_lists);
      Node* node = task.get();
      node = node->parent();
      node = node->parent();
      while(node) {

         string ecf_file_location =  ecf_lists + node->absNodePath() + "/" + task->name() + File::ECF_EXTN();
         create_ecf_file(ecf_file_location);

         //cout << "Node: " << node->absNodePath() << "\n";
         located_ecf_file_fail(task,__LINE__);

         // change look up method
         suite->add_variable( "ECF_FILES_LOOKUP", "prune_leaf" );
         located_ecf_file(task,EcfFile::ECF_FILES,EcfFile::PRUNE_LEAF,__LINE__);

         // Remove this file, required for the following test
         boost::filesystem::remove(ecf_file_location);

         // Remove ECF_FILES_LOOKUP variable
         suite->deleteVariable("ECF_FILES_LOOKUP");

         node = node->parent();
      }

      /// Remove all the generated files
      boost::filesystem::remove_all( ecf_home + suite->absNodePath() );
   }
}

BOOST_AUTO_TEST_CASE( test_ecf_file_locator )
{
	cout << "ANode:: ...test_ecf_file_locator\n";

	// SET ECF_HOME
   std::string smshome = File::test_data("ANode/test/data/SMSHOME","ANode");

	// Create a defs file corresponding to:
 	//# Test the sms file can be found via ECF_SCRIPT
	//#
	//suite suite
	//	edit ECF_INCLUDE $ECF_HOME/includes
	//	edit SLEEPTIME 10
	//	family family
	//   		task t1
	//   		task t2
	//    		task t3
	//  	endfamily
	//endsuite
	//
	//#
	//# This test suite should force a backwards search since the ecf files
	//# are located in ECF_HOME
	//suite suite1
	//  	family family
	//   		task suite1_task1
	//   		task suite1_task2
	//    		task suite1_task3
	//    endfamily
	//endsuite
	//
	//#
	//# This suite is used to test command substition with ECF_FETCH command
	//#
	//suite suite2
	//	edit ECF_INCLUDE $ECF_HOME/includes
	//	edit ECF_FILES   $ECF_HOME
	//	family family
	//		edit SMSFETCH "smsfetch -F %ECF_FILES% -I %ECF_INCLUDE%"
	//		task t2
	//	endfamily
	//endsuite
 	Defs theDefs;
 	{
 	   suite_ptr suite = theDefs.add_suite("suite");
		suite->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
 		suite->addVariable( Variable( "SLEEPTIME", "10" ) );
 		family_ptr fam = suite->add_family( "family" );
		fam->add_task( "t1" );
		fam->add_task( "t2" );
		fam->add_task( "t3" );
 	}
 	{
      suite_ptr suite = theDefs.add_suite("suite1");
      family_ptr fam = suite->add_family( "family" );
      fam->add_task( "suite1_task1" );
      fam->add_task( "suite1_task2" );
      fam->add_task( "suite1_task3" );
 	}
 	{
      suite_ptr suite2 = theDefs.add_suite("suite2");
		suite2->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
 		suite2->addVariable( Variable( Str::ECF_FILES(), "$ECF_HOME" ) );
      family_ptr fam = suite2->add_family( "family" );
		fam->addVariable( Variable( Str::ECF_FETCH(), "smsfetch -F %ECF_FILES% -I %ECF_INCLUDE%" ) );
 		fam->add_task( "t2" );
 	}
   {
      suite_ptr suite = theDefs.add_suite("suite3");
      suite->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
      suite->addVariable( Variable( Str::ECF_FILES(), "$ECF_HOME" ) );
      family_ptr fam = suite->add_family( "family" );
      fam->addVariable( Variable( "ECF_SCRIPT_CMD", "script_cmd -F %ECF_FILES% -I %ECF_INCLUDE%" ) );
      fam->add_task( "t2" );
   }
// 	cerr << theDefs << "\n";

	// get all the task, assume non hierarchical families
	std::vector<Task*> theTasks;
 	theDefs.getAllTasks(theTasks);
	BOOST_REQUIRE_MESSAGE(theTasks.size() == 8, "Expected 8 tasks but found, " << theTasks.size());

	// Override ECF_HOME.   ECF_HOME is need to locate to the ecf files
	theDefs.set_server().add_or_update_user_variables(Str::ECF_HOME(),smshome);

	/// begin , will cause creation of generated variables. The generated variables
	/// are use in client scripts and used to locate the sms files
	theDefs.beginAll();

	// Test for ECF_ file location
  	for(Task* t: theTasks) {
  		try {
  			EcfFile ecf_file = t->locatedEcfFile();
  		   //cout << "Task: " << t->absNodePath() << " " << ecf_file.ecf_file_origin_dump() << "\n";
         BOOST_REQUIRE_MESSAGE( ecf_file.valid(), "Could not locate ecf file for task ");
         Suite* suite = t->suite();
         if (suite->name() == "suite")   BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_SCRIPT, "Expected ecf_file to be obtained from ECF_SCRIPT " << ecf_file.ecf_file_origin_dump());
         if (suite->name() == "suite1")  BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_HOME, "Expected ecf_file to be obtained from ECF_HOME "<< ecf_file.ecf_file_origin_dump());
         if (suite->name() == "suite2")  BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_FETCH_CMD, "Expected ecf_file to be obtained from ECF_FETCH_CMD "<< ecf_file.ecf_file_origin_dump());
         if (suite->name() == "suite3")  BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_SCRIPT_CMD, "Expected ecf_file to be obtained from ECF_SCRIPT_CMD "<< ecf_file.ecf_file_origin_dump());
  		}
  		catch (std::exception& e) {
  			BOOST_REQUIRE_MESSAGE(false,"Could not locate ecf file for task " << e.what());
  		}
	}
}

BOOST_AUTO_TEST_CASE( test_ecf_file_locator_using_ECF_FILES )
{
   cout << "ANode:: ...test_ecf_file_locator_using_ECF_FILES\n";

   // This test will check we can locate the ecf files in ECF_FILES directory

   // SET ECF_HOME. *** TO A DIRECTORY WHERE THERE ARE NO .ECF FILES **
   std::string smshome = File::test_data("ANode/test/data","ANode");

   // set ECF_FILES where we do have some .ecf files
   std::string ecf_files = File::test_data("ANode/test/data/SMSHOME/suite/family","ANode");

   // Create a defs file corresponding to:
   //# Test the sms file can be found via ECF_SCRIPT
   //#
   //edit ECF_HOME    ANode/test/data
   //edit ECF_FILES   ANode/test/data/SMSHOME/suite/family
   //suite suite
   // edit ECF_INCLUDE $ECF_HOME/includes
   // edit SLEEPTIME 10
   // family family
   //       task t1
   //       task t2
   //       task t3
   Defs theDefs;
   {
      suite_ptr suite = theDefs.add_suite("suite");
      suite->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
      suite->addVariable( Variable( "SLEEPTIME", "10" ) );
      family_ptr fam = suite->add_family( "family" );
      fam->add_task( "t1" );
      fam->add_task( "t2" );
      fam->add_task( "t3" );
   }

   // get all the task, assume non hierarchical families
   std::vector<Task*> theTasks;
   theDefs.getAllTasks(theTasks);
   BOOST_REQUIRE_MESSAGE(theTasks.size() == 3, "Expected 3 tasks but found, " << theTasks.size());

   // ECF_HOME, a directory with no .ecf files
   theDefs.set_server().add_or_update_user_variables(Str::ECF_HOME(),smshome);
   theDefs.set_server().add_or_update_user_variables(Str::ECF_FILES(),ecf_files);

//    cerr << theDefs << "\n";

   /// begin , will cause creation of generated variables. The generated variables
   /// are use in client scripts and used to locate the sms files
   theDefs.beginAll();

   // Test for ECF_ file location
   for(Task* t: theTasks) {
      try {
         EcfFile ecf_file = t->locatedEcfFile();
         //cout << "Task: " << t->absNodePath() << " " << ecf_file.ecf_file_origin_dump() << "\n";
         BOOST_REQUIRE_MESSAGE( ecf_file.valid(), "Could not locate ecf file for task " << t->absNodePath());
         BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_FILES, "Expected ecf_file to be obtained from ECF_FILES "<< ecf_file.ecf_file_origin_dump());
      }
      catch (std::exception& e) {
         BOOST_REQUIRE_MESSAGE(false,"Could not locate ecf file for task\n" << e.what());
      }
   }
}

BOOST_AUTO_TEST_CASE( test_ecf_file_locator_using_ECF_FILES_variable_substitution )
{
   cout << "ANode:: ...test_ecf_file_locator_using_ECF_FILES_variable_substitution\n";

   // ECFLOW-788
   // This test will check we can locate the ecf files in ECF_FILES directory
   // Will attempt variable substitution, on ECF_FILES is directory does not exist

   // SET ECF_HOME. ***TO A DIRECTORY WHERE THERE ARE NO .ECF FILE **
   std::string smshome = File::test_data("ANode/test/data","ANode");

   // set ECF_FILES where we do have some .ecf files
   std::string ecf_files = File::test_data("ANode/test/data/SMSHOME/suite/%FAMILY%","ANode");

   // Create a defs file corresponding to:
   //# Test the sms file can be found via ECF_SCRIPT
   //#
   //edit ECF_HOME    ANode/test/data
   //edit ECF_FILES   ANode/test/data/SMSHOME/suite/%FAMILY%   # test we variable substitute
   //suite suite
   // edit ECF_INCLUDE $ECF_HOME/includes
   // edit SLEEPTIME 10
   // family family
   //       task t1
   //       task t2
   //       task t3
   Defs theDefs;
   {
      suite_ptr suite = theDefs.add_suite("suite");
      suite->addVariable( Variable( Str::ECF_INCLUDE(), "$ECF_HOME/includes" ) );
      suite->addVariable( Variable( "SLEEPTIME", "10" ) );
      family_ptr fam = suite->add_family( "family" );
      fam->addVariable( Variable( "FAMILY", "family" ) );
      fam->add_task( "t1" );
      fam->add_task( "t2" );
      fam->add_task( "t3" );
   }

   // get all the task, assume non hierarchical families
   std::vector<Task*> theTasks;
   theDefs.getAllTasks(theTasks);
   BOOST_REQUIRE_MESSAGE(theTasks.size() == 3, "Expected 3 tasks but found, " << theTasks.size());

   // ECF_HOME, a directory with no .ecf files
   theDefs.set_server().add_or_update_user_variables(Str::ECF_HOME(),smshome);
   theDefs.set_server().add_or_update_user_variables(Str::ECF_FILES(),ecf_files);

//    cerr << theDefs << "\n";

   /// begin , will cause creation of generated variables. The generated variables
   /// are use in client scripts and used to locate the sms files
   theDefs.beginAll();

   // Test for ECF_ file location
   for(Task* t: theTasks) {
      try {
         EcfFile ecf_file = t->locatedEcfFile();
         //cout << "Task: " << t->absNodePath() << " " << ecf_file.ecf_file_origin_dump() << "\n";
         BOOST_REQUIRE_MESSAGE( ecf_file.valid(), "Could not locate ecf file for task " << t->absNodePath());
         BOOST_REQUIRE_MESSAGE( ecf_file.ecf_file_origin() == EcfFile::ECF_FILES, "Expected ecf_file to be obtained from ECF_FILES "<< ecf_file.ecf_file_origin_dump());
      }
      catch (std::exception& e) {
         BOOST_REQUIRE_MESSAGE(false,"Could not locate ecf file for task\n" << e.what());
      }
   }
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: mpartio/ecflow

max_issues_repo_name: mpartio/ecflow

max_forks_repo_name: mpartio/ecflow


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 <unordered_map>

#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/CompoundSampleFilter.hpp"
#include "genfile/SampleFilterNegation.hpp"
#include "genfile/VariableInSetSampleFilter.hpp"
#include "genfile/GPThresholdingGTSetter.hpp"
#include "genfile/db/Error.hpp"
#include "genfile/CartesianProductVisitor.hpp"

#include "metro/constants.hpp"
#include "metro/regression/Design.hpp"
#include "metro/SampleRange.hpp"
#include "metro/intersect_ranges.hpp"
#include "metro/count_range.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/FishersExactTest.hpp"

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

#include <Eigen/Core>
#include <Eigen/QR>
#include "Eigen/Eigenvalues"

// #define DEBUG 1

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

namespace {
	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 ;
	}

}

struct CorrelatorOptions: 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[ "-g1" ]
			.set_description( 	"Path to first genotype file."
								"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[ "-g1-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[ "-g1-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.option_implies_option( "-g1-incl-range", "-g1" ) ;
		options.option_implies_option( "-g1-incl-rsids", "-g1" ) ;

		options[ "-g2" ]
			.set_description( 	"Path to second genotype file.  If not given the first genotype file will be used."
								"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( 0 )
			.set_maximum_multiplicity( 1 ) ;

		options[ "-g2-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[ "-g2-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.option_implies_option( "-g2-incl-range", "-g2" ) ;
		options.option_implies_option( "-g2-incl-rsids", "-g2" ) ;

		options[ "-s" ]
			.set_description( "Path of sample file" )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 1 )
			.set_maximum_multiplicity( 1 ) ;

		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[ "-details" ]
			.set_description( "Output counts and a Fisher's exact test P-value. "
				"Currently only supported if genotypes are haploid." ) ;

		options.declare_group( "Model 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_single_value()
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-min-r2" ]
			.set_description( "Do not output results where r^2 is lower than this threshold. "
					" (These values are still computed and contribute to the histogram.)" )
			.set_takes_single_value()
			.set_default_value( 0.05 ) ;
		options[ "-min-maf" ]
			.set_description( "Do not compute correlation results where the maf of either variant is lower than this threshold." )
			.set_takes_single_value()
			.set_default_value( 0 ) ;
		options[ "-min-N" ]
			.set_description( "Do not compute results where the total number of pairwise non-missing genotypes is below this threshold." )
			.set_takes_single_value()
			.set_default_value( 0 ) ;
		options[ "-min-N-propn" ]
			.set_description( "Ignore (do not compute or output) results where the propn of pairwise non-missing genotypes is below this threshold." )
			.set_takes_single_value()
			.set_default_value( 0.0 ) ;
		options[ "-min-distance" ]
			.set_description( "Ignore (do not compute or output) results where the variants are on the same chromosome and closer than this physical distance." )
			.set_takes_single_value()
			.set_default_value( 0 ) ;
		options[ "-assume-haploid" ]
			.set_description( "Convert all data to haploid calls.  This converts homozygous calls"
				" to haploid calls, and treats any heterozygous calls as missing.")
		;
		options[ "-prior-weight" ]
			.set_description( "Specify a prior weight for samples.  "
				" This augments data with one of each of the four possible haploid genotype combinations "
				" for each pairwise comparison, with the total weight given."
				" For example, a value of 1 indicates each haplotype is given 1/4 weight."
			)
			.set_takes_single_value()
			.set_default_value( "0" )
		;
		options.declare_group( "Miscellaneous options" ) ;
		options[ "-analysis-name" ]
			.set_description( "Specify a name to label results from this analysis with." )
			.set_takes_single_value()
			.set_default_value( "ldbird 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[ "-table-prefix" ]
			.set_description( "Specify a prefix to add to tables.  They will be called <prefix>Frequency and <prefix>R." )
			.set_takes_single_value()
			.set_default_value( "" ) ;
		options[ "-threshold" ]
			.set_description( "The threshold to apply to genotype probabilities, if necessary, to make calls." )
			.set_takes_single_value()
			.set_default_value( 0.9 ) ;
		options.declare_group( "Miscellaneous options" ) ;
		options[ "-debug" ]
			.set_description( "Output debugging information." ) ;
	}
} ;

namespace {
	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_nonmissing_samples( nonmissing_samples ),
			m_sample_i(0)
		{
			assert( genotypes != 0 && nonmissing_samples != 0 && ploidy != 0 ) ;
		}

		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"
				) ;
			}
			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 n,
			genfile::OrderType const order_type,
			genfile::ValueType const value_type
		) {
			if( value_type != genfile::eProbability ) {
				throw genfile::BadArgumentError(
					"Callsetter::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(
					"Callsetter::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(
					"Callsetter::set_number_of_entries()",
					"order_type=" + genfile::string_utils::to_string( order_type ),
					"Expected a diploid sample."
				) ;
			}
			(*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.
			if( m_sample_i > m_last_nonmissing_sample_i ) {
				m_nonmissing_samples->push_back(
					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
				) ;
			}
			// Skip this sample for next range.
			m_last_nonmissing_sample_i = m_sample_i + 1 ;
		}

		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 ) {
				(*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 ;
		SampleRanges* m_nonmissing_samples ;
		std::size_t m_last_nonmissing_sample_i ;
		std::size_t m_sample_i ;
	} ;
}

struct LDbirdApplication: public appcontext::ApplicationContext
{
public:

	typedef Eigen::MatrixXd Matrix ;
	typedef Eigen::VectorXd Vector ;
	typedef Eigen::VectorXi IntegerVector ;
	typedef Eigen::RowVectorXd RowVector ;

public:
	LDbirdApplication( int argc, char** argv ):
		appcontext::ApplicationContext(
			globals::program_name,
			globals::program_version + ", revision " + globals::program_revision,
			std::auto_ptr< appcontext::OptionProcessor >( new CorrelatorOptions ),
			argc,
			argv,
			"-log"
		)
	{
		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 ) ;
		}
	}

	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" ) ) ;

		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()
			) ;
		}

		genfile::SNPDataSource::UniquePtr g1 = open_genotype_data_sources(
			options().get< std::string >( "-g1" ),
			get_variant_filter(
				"-g1-incl-range",
				"-g1-incl-rsids"
			),
			excluded_samples
		) ;

		genfile::SNPDataSource::UniquePtr g2 ;
		if( options().check( "-g2" )) {
			// Two files given, we'll do a full carThis mechanism chooses a full cartesian product or a lower triangle implementation.
			// TODO: make this more obvious / cleaner.
			g2 = open_genotype_data_sources(
				options().get< std::string >( "-g2" ),
				get_variant_filter(
					"-g2-incl-range",
					"-g2-incl-rsids"
				),
				excluded_samples
			) ;
		} else {
			g2 = open_genotype_data_sources(
				options().get< std::string >( "-g1" ),
				get_variant_filter(
					"-g1-incl-range",
					"-g1-incl-rsids"
				),
				excluded_samples
			) ;
		}

		std::size_t const N = samples->size() ;
		if( g1->number_of_samples() != N ) {
			throw genfile::BadArgumentError(
				"LDbirdApplication::unsafe_process()",
				"-g1 \"" + options().get< std::string >( "-g1" ) + "\"",
				"Wrong number of samples (" + to_string(g1->number_of_samples()) + ", expected " + to_string(N) + ")"
			) ;
		}
		if( g2->number_of_samples() != N ) {
			throw genfile::BadArgumentError(
				"LDbirdApplication::unsafe_process()",
				"-g2 \"" + options().get< std::string >( "-g2" ) + "\"",
				"Wrong number of samples (" + to_string(g2->number_of_samples()) + ", expected " + to_string(N) + ")"
			) ;
		}

		write_preamble( *g1, *g2, *samples ) ;

		std::vector< int64_t > histogram( 2049, 0 ) ;

		boost::optional< genfile::db::Connection::RowId > analysis_id ;
		if( options().check( "-analysis-id" ) ) {
			analysis_id = options().get< genfile::db::Connection::RowId >( "-analysis-id" ) ;
		}

		std::string const fileSpec = options().get< std::string > ( "-o" ) ;
		std::string const tablePrefix = options().get< std::string > ( "-table-prefix" ) ;

		{
			qcdb::Storage::UniquePtr frequencyStorage = qcdb::Storage::create(
				fileSpec + ":" + tablePrefix + "Frequency",
				options().get< std::string > ( "-analysis-name" ),
				options().get< std::string > ( "-analysis-chunk" ),
				get_application_metadata(),
				"position,alleles",
				analysis_id
			) ;
			analysis_id = frequencyStorage->analysis_id() ;

			frequencyStorage->add_variable( "number_of_haplotypes" ) ;
			frequencyStorage->add_variable( "frequency" ) ;
			//frequencyStorage->set_variant_names( std::vector< std::string >({ "g1", "g2" })) ;

			qcdb::MultiVariantStorage::UniquePtr correlationStorage = qcdb::MultiVariantStorage::create(
				fileSpec + ":" + tablePrefix + "R+without rowid",
				2,
				options().get< std::string > ( "-analysis-name" ),
				options().get< std::string > ( "-analysis-chunk" ),
				get_application_metadata(),
				"position,alleles",
				analysis_id
			) ;

			correlationStorage->set_variant_names( std::vector< std::string >({ "g1", "g2" })) ;
			correlationStorage->add_variable( "N" ) ;
			correlationStorage->add_variable( "encoded_r" ) ;

			if( options().check( "-details" )) {
				for( int i = 0; i < 2; ++i ) {
					for( int j = 0; j < 2; ++j ) {
						correlationStorage->add_variable( "n_" + to_string(i) + to_string(j) ) ;
					}
				}

				correlationStorage->add_variable( "fet_pvalue" ) ;
			}
			// This mechanism chooses a full cartesian product or a lower triangle implementation
			// TODO: make this more obvious / cleaner.
			genfile::CartesianProductVisitor visitor( !options().check( "-g2" ) ) ;
			visitor.add_source( "g1", g1.get() ) ;
			visitor.add_source( "g2", g2.get() ) ;

			if( options().check( "-min-distance" )) {
				visitor.set_min_distance( options().get< int64_t >( "-min-distance" )) ;
			}

			run( visitor, *samples, *frequencyStorage, *correlationStorage, &histogram ) ;

			frequencyStorage->finalise() ;
			correlationStorage->finalise() ;
		}

		{
			qcdb::MultiVariantStorage::UniquePtr histogramStorage = qcdb::MultiVariantStorage::create(
				fileSpec + ":" + tablePrefix + "Histogram",
				0,
				options().get< std::string > ( "-analysis-name" ),
				options().get< std::string > ( "-analysis-chunk" ),
				get_application_metadata(),
				"position,alleles",
				analysis_id
			) ;
			histogramStorage->add_variable( "encoded_r" ) ;
			histogramStorage->add_variable( "count" ) ;
			qcdb::MultiVariantStorage::Key const key ;
			for( std::size_t i = 0; i < histogram.size(); ++i ) {
				histogramStorage->create_new_key( key ) ;
				histogramStorage->store_data_for_key( key, "encoded_r", int64_t(i) ) ;
				histogramStorage->store_data_for_key( key, "count", histogram[i] ) ;
			}
			histogramStorage->finalise() ;
		}
	}

	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() ) ;
	}

	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() << "  First dataset:\n" << host.get_summary() << "\n" ;
		get_ui_context().logger() << " Second dataset:\n" << para.get_summary() << "\n" ;
	}

	void run(
		genfile::MultiSourceVisitor& visitor,
		genfile::CohortIndividualSource const& samples,
		qcdb::Storage& frequencyOutput,
		qcdb::MultiVariantStorage& correlationOutput,
		std::vector< int64_t >* histogram
	) {
		appcontext::UIContext::ProgressContext progress_context = get_ui_context().get_progress_context( "Testing" ) ;
		histogram->resize( 2049 ) ;
		std::fill( histogram->begin(), histogram->end(), 0 ) ;
		for(
			std::size_t count = 0;
			visitor.step(
				[this,&frequencyOutput,&correlationOutput,&histogram](
					std::vector< int > const& changed,
					std::vector< genfile::VariantIdentifyingData > const& variants,
					std::vector< genfile::VariantDataReader::SharedPtr > const& readers
				) {
					this->process_one(
						changed,
						variants,
						readers,
						frequencyOutput,
						correlationOutput,
						histogram
					) ;
				}
			) ;
			++count
		) {
			progress_context( count, visitor.count() ) ;
		} ;
		progress_context.finish() ;
	}

	void process_one(
		std::vector< int > const& changed,
		std::vector< genfile::VariantIdentifyingData > const& variants,
		std::vector< genfile::VariantDataReader::SharedPtr > const& readers,
		qcdb::Storage& frequencyOutput,
		qcdb::MultiVariantStorage& correlationOutput,
		std::vector< int64_t >* histogram
	) {
		using genfile::string_utils::to_string ;

		if( m_dosages.size() != changed.size() ) {
			m_dosages.resize( changed.size() ) ;
			m_ploidy.resize( changed.size() ) ;
			m_nonmissingness.resize( changed.size() ) ;
		}
		double const prior_weight = options().get< double >( "-prior-weight" ) ;
		double min_maf = 0.5 ;
		for( std::size_t i = 0; i < changed.size(); ++i ) {
			if( changed[i] == 1 ) {
				DosageSetter setter( &(m_dosages[i]), &(m_ploidy[i]), &(m_nonmissingness[i]), options().check( "-assume-haploid" )) ;
				readers[i]->get( ":genotypes:", genfile::to_GP_unphased( setter )) ;
			}

#if DEBUG
			std::cerr << globals::program_name + ":process_one(): variant " << i << ": " << variants[i] << ":\n" ;
			std::cerr << "Loaded data with nonmissingness: " << m_nonmissingness[i] << ".\n" ;
#endif

			FrequencyStore::const_iterator where = m_frequencies.find( variants[i] ) ;
			double frequency = 0.0 ;
			if( where != m_frequencies.end() ) {
				frequency = where->second ;
			} else {
				std::pair< double, int64_t > computedFrequency = compute_regularised_frequency( m_dosages[i], m_ploidy[i], m_nonmissingness[i], prior_weight ) ;
				frequency = computedFrequency.first ;

				// not yet computed, so store it
				m_frequencies[ variants[i] ] = frequency ;
				// Output the sample count and the frequency
				frequencyOutput.store_per_variant_data(
					variants[i],
					"number_of_haplotypes",
					computedFrequency.second
				) ;
				frequencyOutput.store_per_variant_data(
					variants[i],
					"frequency",
					frequency
				) ;
			}
			min_maf = std::min( min_maf, std::min( frequency, 1.0 - frequency )) ;
		}

		// Bail out if any of the variants are too rare
		if( min_maf < options().get< double >( "-min-maf" )) {
			return ;
		}

		int const total_number_of_samples = m_dosages[0].size() ;
		std::vector< metro::SampleRange > included_samples( 1, metro::SampleRange( 0, total_number_of_samples )) ;
		for( std::size_t i = 0; i < changed.size(); ++i ) {
			included_samples = metro::impl::intersect_ranges( included_samples, m_nonmissingness[i] ) ;
		}

		// Bail out if there aren't enough samples in the pairwise comparison
		{
			std::size_t count = metro::impl::count_range( included_samples ) ;
			if(
				(count < options().get< double >( "-min-N" ))
			 	   || (( double(count)/total_number_of_samples ) < options().get< double >( "-min-N-propn" ))
			) {
				return ;
			}
		}

		double covariance = 0.0, correlation = 0.0, N = 0.0 ;

		compute_regularised_correlation(
			m_dosages, m_ploidy, included_samples,
			&covariance, &correlation, &N,
			prior_weight
		) ;

		// ints in sqlite use only 2 bytes if up to +ve integer 2287.
		// We map -1...1 to 0...2048, such that the transformation
		//
		// correlation = (stored_value / 1024) - 1.0
		// or
		// correlation = (stored_value - 1024.0) / 1024.0
		// maps back to correlation space.
		//
		int64_t const encoded_r = ( std::round((correlation + 1.0) * 1024 )) ;
		++((*histogram)[ encoded_r ]) ;

		if( ( correlation * correlation ) >= options().get< double >( "-min-r2" ) ) {
			correlationOutput.store_data_for_key(
				variants,
				"N",
				int64_t(N)
			) ;
			// We store correlations as integers on the scale -16384 to 16384.
			// This is because sqlite has compression for integer storage, such that
			// reals take 8 bytes but integers take
			correlationOutput.store_data_for_key(
				variants,
				"encoded_r",
				encoded_r
			) ;

			if( options().check( "-details" )) {
				Eigen::MatrixXd const table = tabulate( m_dosages, m_ploidy, included_samples ) ;
				for( int i = 0; i < table.rows(); ++i ) {
					for( int j = 0; j < table.cols(); ++j ) {
						correlationOutput.store_data_for_key(
							variants,
							"n_" + to_string(i) + to_string(j),
							table(i,j)
						) ;
					}
				}

				if( m_ploidy[0].maxCoeff() == 1 ) {
					metro::FishersExactTest test( table ) ;
					double pvalue = test.get_pvalue( metro::FishersExactTest::eTwoSided ) ;
					correlationOutput.store_data_for_key(
						variants,
						"fet_pvalue",
						pvalue
					) ;
				}
			}
		}
#if DEBUG
		std::cerr << globals::program_name + ":process_one():\n" ;
		std::cerr << "dosage1: " << m_dosages[0].head( 15 ).transpose() << "...\n" ;
		std::cerr << "dosage2: " << m_dosages[1].head( 15 ).transpose() << "...\n" ;
		std::cerr << "N: " << N << ".\n" ;
		std::cerr << "covariance: " << covariance << ".\n" ;
		std::cerr << "correlation: " << correlation << ".\n" ;
		std::cerr << "included: " << included_samples << ".\n" ;
#endif

	}

private:
	std::vector< Eigen::VectorXd > m_dosages ;
	std::vector< Eigen::VectorXd > m_ploidy ;
	std::vector< std::vector< metro::SampleRange > > m_nonmissingness ;
	typedef std::map< genfile::VariantIdentifyingData, double > FrequencyStore ;
	FrequencyStore m_frequencies ;

private:

	std::pair< double, int64_t > compute_regularised_frequency(
		Eigen::VectorXd const& dosages,
		Eigen::VectorXd const& ploidy,
		std::vector< metro::SampleRange > const& included_samples,
		double prior_weight = 1.0
	) {
		double result = 0.0 ;
		double N = 0.0 ;
		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
			metro::SampleRange const& range = included_samples[i] ;
			result += dosages.segment( range.begin(), range.size() ).sum() ;
			N += ploidy.segment( range.begin(), range.size() ).sum() ;
		}
		// Assume we have additional data adding up to prior_weight haploid samples
		// evenly split between both alleles
		result += prior_weight / 2 ;
		N += prior_weight ;
		return std::make_pair( result / N, int64_t( N ) ) ;
	}

	void compute_regularised_correlation(
		std::vector< Eigen::VectorXd > const& dosages,
		std::vector< Eigen::VectorXd > const& ploidy,
		std::vector< metro::SampleRange > const& included_samples,
		double* covariance,
		double* correlation,
		double* number_of_samples,
		double const prior_weight = 1.0
	) {
		assert( dosages.size() == 2 ) ;
		assert( covariance ) ;
		assert( correlation ) ;
		assert( number_of_samples ) ;
		assert( ploidy[0] == ploidy[1] ) ;
		std::vector< double > frequencies( dosages.size(), 0.0 ) ;

		typedef Eigen::VectorBlock< Eigen::VectorXd const > ConstBlock ;
		typedef Eigen::VectorBlock< Eigen::VectorXi const > ConstIntegerBlock ;

		// Although frequencies are computed above, we recompute here because the
		// number of included samples may differ for the pairwise test.
		for( std::size_t v = 0; v < frequencies.size(); ++v ) {
			frequencies[v] = compute_regularised_frequency(
				dosages[v], ploidy[v], included_samples, prior_weight
			).first ;
		}

		double result = 0.0 ;
		double N = 0.0 ;
		std::vector< double > variances( dosages.size(), 0.0 ) ;
		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
			metro::SampleRange const& range = included_samples[i] ;
			ConstBlock d0 = dosages[0].segment( range.begin(), range.size() ) ;
			ConstBlock d1 = dosages[1].segment( range.begin(), range.size() ) ;
			ConstBlock p = ploidy[0].segment( range.begin(), range.size() ) ;
			result += (
				( d0 - p * frequencies[0] ).array()
				* ( d1 - p * frequencies[1] ).array()
			).sum() ;
			N += p.array().sum() ;
			for( std::size_t v = 0; v < frequencies.size(); ++v ) {
				ConstBlock d = dosages[v].segment( range.begin(), range.size() ) ;
				ConstBlock p = ploidy[v].segment( range.begin(), range.size() ) ;
				variances[v] += ( d - p * frequencies[v] ).array().square().sum() ;
			}
		}

		// Add regularising information consisting of a total of
		// prior_weight haploid samples, split evenly between
		// the four possible genotype combinations
		result += (0 - frequencies[0]) * (0 - frequencies[1] ) * (prior_weight / 4) ;
		result += (0 - frequencies[0]) * (1 - frequencies[1] ) * (prior_weight / 4) ;
		result += (1 - frequencies[0]) * (0 - frequencies[1] ) * (prior_weight / 4) ;
		result += (1 - frequencies[0]) * (1 - frequencies[1] ) * (prior_weight / 4) ;
		for( std::size_t v = 0; v < frequencies.size(); ++v ) {
			variances[v] += ( 0 - frequencies[v] ) * ( 0 - frequencies[v] ) * (prior_weight / 2.0) ;
			variances[v] += ( 1 - frequencies[v] ) * ( 1 - frequencies[v] ) * (prior_weight / 2.0) ;
		}
		N += prior_weight ;

		*covariance = result / N ;
		*correlation = result / (std::sqrt( variances[0] ) * std::sqrt( variances[1] )) ;
		*number_of_samples = N ;
	}

	Eigen::MatrixXd tabulate(
		std::vector< Eigen::VectorXd > const& dosages,
		std::vector< Eigen::VectorXd > const& ploidy,
		std::vector< metro::SampleRange > const& included_samples
	) const {
		typedef Eigen::VectorBlock< Eigen::VectorXd const > ConstBlock ;
		typedef Eigen::VectorBlock< Eigen::VectorXi const > ConstIntegerBlock ;

		int maxPloidy0 = 0 ;
		int maxPloidy1 = 0 ;
		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
			metro::SampleRange const& range = included_samples[i] ;
			maxPloidy0 = std::max( maxPloidy0, int( ploidy[0].segment( range.begin(), range.size() ).array().maxCoeff()) ) ;
			maxPloidy1 = std::max( maxPloidy1, int( ploidy[1].segment( range.begin(), range.size() ).array().maxCoeff()) ) ;
		}
		Eigen::MatrixXd result = Eigen::MatrixXd::Zero( maxPloidy0 + 1, maxPloidy1 + 1 ) ;
		for( std::size_t i = 0; i < included_samples.size(); ++i ) {
			metro::SampleRange const& range = included_samples[i] ;
			for( int j = range.begin(); j < range.end(); ++j ) {
				++result( dosages[0][j], dosages[1][j] ) ;
			}
		}
		return result ;
	}

	struct DosageSetter: public genfile::VariantDataReader::PerSampleSetter {
		typedef std::vector< metro::SampleRange > SampleRanges ;

		DosageSetter(
			Eigen::VectorXd* dosages,
			Eigen::VectorXd* ploidy,
			SampleRanges* nonmissing_samples,
			bool treat_as_haploid = false
		):
			m_dosages( dosages ),
			m_ploidy( ploidy ),
			m_total_prob(0.0),
			m_treat_as_haploid( treat_as_haploid ),
			m_number_of_entries( 0 ),
			m_nonmissing_samples( nonmissing_samples ),
			m_sample_i(0)
		{
			assert( dosages != 0 && nonmissing_samples != 0 && ploidy != 0 ) ;
		}

		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"
				) ;
			}
			m_dosages->resize( nSamples ) ;
			m_dosages->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::eProbability ) {
				throw genfile::BadArgumentError(
					"Callsetter::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(
					"Callsetter::set_number_of_entries()",
					"order_type=" + genfile::string_utils::to_string( order_type ),
					"Expected genotype call probabilities (e.g. GP field)."
				) ;
			}
			if( m_treat_as_haploid ) {
				(*m_ploidy)(m_sample_i) = 1 ;
			} else {
				(*m_ploidy)(m_sample_i) = ploidy ;
			}
			m_number_of_entries = n ;
			m_total_prob = 0.0 ;
		}

		void set_value( std::size_t entry_i, genfile::MissingValue const value ) {
			// This sample has missing data, end our current sample range here.
			if( m_sample_i > m_last_nonmissing_sample_i ) {
				m_nonmissing_samples->push_back(
					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
				) ;
#if DEBUG > 2
				std::cerr << "DosageSetter::set_value(): added sample range: " << m_nonmissing_samples->back() << ".\n" ;
#endif
			}
			// Skip this sample for next range.
			m_last_nonmissing_sample_i = m_sample_i + 1 ;
		}

		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( !sample_recorded_missing() ) {
				bool const last_entry = ((entry_i+1) == m_number_of_entries ) ;
				// we assume two alleles, so entry_i is also the count of the B allele
				// in the genotype.
				if( m_treat_as_haploid ) {
					// only homozygous genotypes are counted, others are assumed missing
					// entry_i == 0 contributes 0 dosage anyway so just count the last entry
					//
					// Sample considered missing if has nonzero probability on non-homozygous genotype
					// (caught here), or if it has total 0 probability on homozygous genotypes (handled below).
					if( entry_i == 0 || last_entry ) {
						(*m_dosages)(m_sample_i) += ( last_entry ? value : 0.0 );
						m_total_prob += value ;
					} else if( value != 0 ) {
						record_missing_sample() ;
					}
				} else {
					(*m_dosages)(m_sample_i) += value * entry_i ;
					m_total_prob += value ;
				}

				// handle missing data encoded as zeroes
				if( last_entry && m_total_prob == 0.0 ) {
					record_missing_sample() ;
				}

#if DEBUG > 2
				std::cerr << "DosageSetter::set_value(): added value " << (value * entry_i) << " to entry " << entry_i << " for sample " << m_sample_i << ".\n" ;
#endif
			}
		}

		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 )
				) ;
#if DEBUG > 2
				std::cerr << "DosageSetter::finalise(): added sample range: " << m_nonmissing_samples->back() << ".\n" ;
#endif
			}
#if DEBUG > 2
				std::cerr << "DosageSetter::finalise(): complete.\n" ;
#endif
		}

	private:
		Eigen::VectorXd* m_dosages ;
		Eigen::VectorXd* m_ploidy ;
		double m_total_prob ;
		bool m_treat_as_haploid ;

		std::size_t m_number_of_entries ;
		SampleRanges* m_nonmissing_samples ;
		std::size_t m_last_nonmissing_sample_i ;
		std::size_t m_sample_i ;
	private:

		bool sample_recorded_missing() const {
			return m_sample_i < m_last_nonmissing_sample_i ;
		}

		void record_missing_sample() {
			// missing
			m_nonmissing_samples->push_back(
				metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
			) ;
			// Skip this sample for next range.
			m_last_nonmissing_sample_i = m_sample_i + 1 ;
		}
	} ;
} ;

int main( int argc, char** argv ) {
    try {
		LDbirdApplication app( argc, argv ) ;
    }
	catch( appcontext::HaltProgramWithReturnCode const& e ) {
		return e.return_code() ;
	}
	return 0 ;
}

ext: cpp

max_stars_repo_name: gavinband/qctool

max_issues_repo_name: gavinband/qctool

max_forks_repo_name: gavinband/qctool


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)

#ifndef GENFILE_SNP_IDENTIFYING_DATA2_HPP
#define GENFILE_SNP_IDENTIFYING_DATA2_HPP

#include <string>
#include <vector>
#include <boost/function.hpp>
#include "genfile/GenomePosition.hpp"
#include "genfile/string_utils/slice.hpp"
#include "genfile/SNPIdentifyingData.hpp"

namespace genfile {
	struct VariantIdentifyingData
	{
	public:
		VariantIdentifyingData( std::string const& rsid = "(unknown variant)" ) ;
		VariantIdentifyingData(
			std::string const& RSID,
			GenomePosition const& position,
			std::string const& first_allele,
			std::string const& second_allele
		) ;
	 	VariantIdentifyingData(
			std::string const& SNPID,
			std::string const& rsid,
			GenomePosition const& position,
			std::string const& first_allele,
			std::string const& second_allele
		) ;
		VariantIdentifyingData(
			SNPIdentifyingData const& snp
		) ;
		VariantIdentifyingData( VariantIdentifyingData const& other ) ;
		VariantIdentifyingData& operator=( VariantIdentifyingData const& other ) ;
		VariantIdentifyingData& operator=( SNPIdentifyingData const& other ) ;

		operator SNPIdentifyingData() const ;

		typedef string_utils::slice slice ;

		void clear_identifiers() ;
		std::size_t number_of_identifiers() const ;
		void set_primary_id( slice const& id ) ;
		void add_identifier( slice const& id ) ;
		std::vector< slice > get_identifiers(
			std::size_t start = 0,
			std::size_t end = std::string::npos
		) const ;
		void get_identifiers(
			boost::function< void( slice ) >,
			std::size_t start = 0,
			std::size_t end = std::string::npos
		) const ;
		std::string get_identifiers_as_string(
			std::string const& separator,
			std::size_t start = 0,
			std::size_t end = std::string::npos
		) const ;
		slice get_primary_id() const { return slice( m_data, m_rsid_start, m_allele_starts[0] ) ; }

		void set_position( GenomePosition const& position ) { m_position = position ;}
		GenomePosition const& get_position() const { return m_position ; }

		void clear_alleles() ;
		void set_allele( std::size_t i, std::string const& allele ) ;
		void set_allele( std::size_t i, slice const& allele ) ;
		void set_allele( std::size_t i, char const* allele ) ;
		void add_allele( slice const& allele ) ;
		std::size_t number_of_alleles() const { return m_allele_starts.size() - 1 ; }
		slice get_allele( std::size_t i ) const ;
		std::vector< slice > get_alleles( std::size_t start = 0, std::size_t end = std::string::npos ) const ;
		void get_alleles( boost::function< void( slice ) >, std::size_t start = 0, std::size_t end = std::string::npos ) const ;
		void swap_alleles() ;

		std::size_t estimate_bytes_used() const ;
	public:
		struct CompareFields {
			CompareFields(
				std::string const& fields_to_compare = "position,rsid,IDs,alleles",
				bool flip_alleleles_if_necessary = false
			) ;
			CompareFields( CompareFields const& other ) ;
			CompareFields& operator=( CompareFields const& other ) ;

			bool get_flip_alleles_if_necessary() const { return m_flip_alleles_if_necessary ; }

			enum { eIDs = 0x1, eRSID = 0x2, ePosition = 0x4, eAlleles = 0x8, eMask = 0xF } ;
			bool operator()( VariantIdentifyingData const& left, VariantIdentifyingData const& right ) const ;
			bool are_equal( VariantIdentifyingData const& left, VariantIdentifyingData const& right ) const ;
			bool check_if_comparable_fields_are_known( VariantIdentifyingData const& value ) const ;
			std::vector< int > const& get_compared_fields() const { return m_fields_to_compare ;}

			std::string get_summary() const ;

		private:
			static std::vector< int > parse_fields_to_compare( std::string const& field_spec ) ;
			std::vector< int > m_fields_to_compare ;
			bool m_flip_alleles_if_necessary ;
		} ;

		friend bool operator==( VariantIdentifyingData const& left, VariantIdentifyingData const& right ) ;
		friend bool operator!=( VariantIdentifyingData const& left, VariantIdentifyingData const& right ) ;
		friend bool operator<( VariantIdentifyingData const& left, VariantIdentifyingData const& right ) ;

	private:
		// we store all IDs and alleles in a string.
		// and hold slices to the string.
		std::string m_data ;
		typedef uint32_t Size ;
		Size m_rsid_start ; // always equals 0 in current implementation.
		std::vector< Size > m_allele_starts ;
		GenomePosition m_position ;

		CompareFields& operator=( CompareFields const& other ) ;
	} ;

	std::ostream& operator<<( std::ostream&, VariantIdentifyingData const& ) ;
	std::ostream& operator<<( std::ostream&, std::vector< VariantIdentifyingData > const& ) ;
}

#endif

ext: hpp

max_stars_repo_name: gavinband/bingwa

max_issues_repo_name: gavinband/bingwa

max_forks_repo_name: gavinband/bingwa


text:
//  Copyright (c) 2013-2014 Hartmut Kaiser
//  Copyright (c) 2015 Andreas Schaefer
//
//  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)

#if !defined(HPX_SERIALIZATION_SERIALIZE_BUFFER_APR_05_2013_0312PM)
#define HPX_SERIALIZATION_SERIALIZE_BUFFER_APR_05_2013_0312PM

#include <hpx/config.hpp>
#include <hpx/runtime/serialization/array.hpp>
#include <hpx/runtime/serialization/serialize.hpp>
#include <hpx/throw_exception.hpp>
#include <hpx/traits/supports_streaming_with_any.hpp>
#include <hpx/util/bind.hpp>

#include <boost/shared_array.hpp>

#include <algorithm>
#include <cstddef>
#include <type_traits>

namespace hpx { namespace serialization
{
    ///////////////////////////////////////////////////////////////////////////
    template <typename T, typename Allocator = std::allocator<T> >
    class serialize_buffer
    {
    private:
        typedef Allocator allocator_type;

        static void no_deleter(T*) {}

        template <typename Deallocator>
        static void deleter(T* p, Deallocator dealloc, std::size_t size)
        {
            dealloc.deallocate(p, size);
        }

    public:
        enum init_mode
        {
            copy = 0,       // constructor copies data
            reference = 1,  // constructor does not copy data and does not
                            // manage the lifetime of it
            take = 2        // constructor does not copy data but does take
                            // ownership and manages the lifetime of it
        };

        typedef T value_type;

        explicit serialize_buffer(allocator_type const& alloc = allocator_type())
          : size_(0)
          , alloc_(alloc)
        {}

        explicit serialize_buffer(std::size_t size,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            using util::placeholders::_1;
            data_.reset(alloc_.allocate(size),
                        util::bind(&serialize_buffer::deleter<allocator_type>,
                                   _1, alloc_, size_));
        }

        // The default mode is 'copy' which is consistent with the constructor
        // taking a T const * below.
        serialize_buffer (T* data, std::size_t size, init_mode mode = copy,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            if (mode == copy) {
                using util::placeholders::_1;
                data_.reset(alloc_.allocate(size),
                    util::bind(&serialize_buffer::deleter<allocator_type>,
                        _1, alloc_, size_));
                if (size != 0)
                    std::copy(data, data + size, data_.get());
            }
            else if (mode == reference) {
                data_ = boost::shared_array<T>(data,
                    &serialize_buffer::no_deleter);
            }
            else {
                // take ownership
                using util::placeholders::_1;
                data_ = boost::shared_array<T>(data,
                    util::bind(&serialize_buffer::deleter<allocator_type>,
                        _1, alloc_, size_));
            }
        }

        template <typename Deallocator>
        serialize_buffer (T* data, std::size_t size,
                allocator_type const& alloc, Deallocator const& dealloc)
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            // if 2 allocators are specified we assume mode 'take'
            using util::placeholders::_1;
            data_ = boost::shared_array<T>(data,
                util::bind(&serialize_buffer::deleter<Deallocator>,
                    _1, dealloc, size_));
        }

        template <typename Deleter>
        serialize_buffer (T* data, std::size_t size, init_mode mode,
                Deleter const& deleter,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            if (mode == copy) {
                data_.reset(alloc_.allocate(size), deleter);
                if (size != 0)
                    std::copy(data, data + size, data_.get());
            }
            else {
                // reference or take ownership, behavior is defined by deleter
                data_ = boost::shared_array<T>(data, deleter);
            }
        }

        template <typename Deleter>
        serialize_buffer (T const* data, std::size_t size, init_mode mode, //-V659
                Deleter const& deleter,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            if (mode == copy) {
                data_.reset(alloc_.allocate(size), deleter);
                if (size != 0)
                    std::copy(data, data + size, data_.get());
            }
            else if (mode == reference) {
                data_ = boost::shared_array<T>(const_cast<T*>(data), deleter);
            }
            else {
                // can't take ownership of const buffer
                HPX_THROW_EXCEPTION(bad_parameter,
                    "serialize_buffer::serialize_buffer",
                    "can't take ownership of const data");
            }
        }

        // Deleter needs to use deallocator
        template <typename Deallocator, typename Deleter>
        serialize_buffer (T* data, std::size_t size,
                allocator_type const& alloc, Deallocator const& dealloc,
                Deleter const& deleter)
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            // if 2 allocators are specified we assume mode 'take'
            data_ = boost::shared_array<T>(data, deleter);
        }

        // same set of constructors, but taking const data
        serialize_buffer (T const* data, std::size_t size,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            // create from const data implies 'copy' mode
            using util::placeholders::_1;
            data_.reset(alloc_.allocate(size),
                util::bind(&serialize_buffer::deleter<allocator_type>,
                    _1, alloc_, size_));
            if (size != 0)
                std::copy(data, data + size, data_.get());
        }

        template <typename Deleter>
        serialize_buffer (T const* data, std::size_t size,
                Deleter const& deleter,
                allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            // create from const data implies 'copy' mode
            data_.reset(alloc_.allocate(size), deleter);
            if (size != 0)
                std::copy(data, data + size, data_.get());
        }

        serialize_buffer (T const* data, std::size_t size,
                init_mode mode, allocator_type const& alloc = allocator_type())
          : data_()
          , size_(size)
          , alloc_(alloc)
        {
            if (mode == copy) {
                using util::placeholders::_1;
                data_.reset(alloc_.allocate(size),
                    util::bind(&serialize_buffer::deleter<allocator_type>,
                        _1, alloc_, size_));
                if (size != 0)
                    std::copy(data, data + size, data_.get());
            }
            else if (mode == reference) {
                data_ = boost::shared_array<T>(
                    const_cast<T*>(data),
                    &serialize_buffer::no_deleter);
            }
            else {
                // can't take ownership of const buffer
                HPX_THROW_EXCEPTION(bad_parameter,
                    "serialize_buffer::serialize_buffer",
                    "can't take ownership of const data");
            }
        }

        // accessors enabling data access
        T* data() { return data_.get(); }
        T const* data() const { return data_.get(); }

        T* begin() { return data(); }
        T* end() { return data() + size_; }

        T& operator[](std::size_t idx) { return data_[idx]; }
        T operator[](std::size_t idx) const { return data_[idx]; }

        boost::shared_array<T> data_array() const { return data_; }

        std::size_t size() const { return size_; }

    private:
        // serialization support
        friend class hpx::serialization::access;

        ///////////////////////////////////////////////////////////////////////
        template <typename Archive>
        void save(Archive& ar, const unsigned int version) const
        {
            ar << size_ << alloc_; //-V128

            if (size_ != 0)
            {
                ar << hpx::serialization::make_array(data_.get(), size_);
            }
        }

        ///////////////////////////////////////////////////////////////////////
        template <typename Archive>
        void load(Archive& ar, const unsigned int version)
        {
            using util::placeholders::_1;
            ar >> size_ >> alloc_; //-V128

            data_.reset(alloc_.allocate(size_),
                util::bind(&serialize_buffer::deleter<allocator_type>, _1,
                    alloc_, size_));

            if (size_ != 0)
            {
                ar >> hpx::serialization::make_array(data_.get(), size_);
            }
        }

        HPX_SERIALIZATION_SPLIT_MEMBER()

        // this is needed for util::any
        friend bool
        operator==(serialize_buffer const& rhs, serialize_buffer const& lhs)
        {
            return rhs.data_.get() == lhs.data_.get() && rhs.size_ == lhs.size_;
        }

    private:
        boost::shared_array<T> data_;
        std::size_t size_;
        Allocator alloc_;
    };
}}

namespace hpx { namespace traits
{
    ///////////////////////////////////////////////////////////////////////////
    // Customization point for streaming with util::any, we don't want
    // serialization::serialize_buffer to be streamable
    template <typename T, typename Allocator>
    struct supports_streaming_with_any<serialization::serialize_buffer<T, Allocator> >
      : std::false_type
    {};
}}

#endif

max_stars_repo_name: brycelelbach/hpx

max_issues_repo_name: atrantan/hpx

max_forks_repo_name: atrantan/hpx


text:
/*
 * Boost.Extension / common:
 *         common include files
 *
 * (C) Copyright Jeremy Pack 2008
 * 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)
 *
 * See http://www.boost.org/ for latest version.
 */

#ifndef BOOST_EXTENSION_COMMON_HPP
#define BOOST_EXTENSION_COMMON_HPP

#include <boost/preprocessor/arithmetic/inc.hpp>
#include <boost/preprocessor/if.hpp>
#include <boost/preprocessor/punctuation/comma_if.hpp>
#include <boost/preprocessor/repetition.hpp>
#include <boost/preprocessor/iteration/iterate.hpp>
#include <boost/concept_check.hpp>
/** This determines the maximum number of parameters that a constructor
  * or exported shared library function can have. 10 is the same default
  * as Boost.Function.
  */
#ifndef BOOST_EXTENSION_MAX_FUNCTOR_PARAMS
#define BOOST_EXTENSION_MAX_FUNCTOR_PARAMS 10
#endif

#endif  // BOOST_EXTENSION_COMMON_HPP

max_stars_repo_name: loafofpiecrust/lfant


text:
#include "stdafx.h"
#include "Library.h"

#include <boost/test/unit_test.hpp>
#ifdef UNIT_TEST
BOOST_AUTO_TEST_SUITE( test_library )

	BOOST_AUTO_TEST_CASE(null_arg_ctor)
	{
		CLibrary library(NULL) ;
		BOOST_CHECK(! library.is_loaded()) ;
	}
	BOOST_AUTO_TEST_CASE(null_arg_load)
	{
		CLibrary library(NULL) ;
		BOOST_CHECK(! library.load(NULL)) ;
	}
	BOOST_AUTO_TEST_CASE( rich_edit)
	{
		CLibrary library(_T("RICHED20.DLL")) ;
		BOOST_CHECK(library.is_loaded()) ;
	}
	BOOST_AUTO_TEST_CASE( get_proc_address)
	{
		CLibrary library(_T("ArmAccess.DLL")) ;
		BOOST_CHECK(library.is_loaded()) ;
		FARPROC proc_address = library.get_proc_address("InstallKeyLater") ;
		BOOST_CHECK(proc_address) ;
	}

BOOST_AUTO_TEST_SUITE_END()

#endif

max_stars_repo_name: ultimatezen/felix

max_issues_repo_name: ultimatezen/felix

max_forks_repo_name: ultimatezen/felix


text:
//
// Copyright (c) 2019 INRIA
//

#ifndef __pinocchio_joint_revolute_unbounded_unaligned_hpp__
#define __pinocchio_joint_revolute_unbounded_unaligned_hpp__

#include "pinocchio/fwd.hpp"
#include "pinocchio/spatial/inertia.hpp"
#include "pinocchio/math/rotation.hpp"
#include "pinocchio/math/matrix.hpp"

#include "pinocchio/multibody/joint/joint-revolute-unaligned.hpp"

namespace pinocchio
{

  template<typename Scalar, int Options = 0> struct JointRevoluteUnboundedUnalignedTpl;

  template<typename _Scalar, int _Options>
  struct traits< JointRevoluteUnboundedUnalignedTpl<_Scalar,_Options> >
  {
    enum {
      NQ = 2,
      NV = 1
    };
    typedef _Scalar Scalar;
    enum { Options = _Options };

    typedef Eigen::Matrix<Scalar,NQ,1,Options> ConfigVector_t;
    typedef Eigen::Matrix<Scalar,NV,1,Options> TangentVector_t;

    typedef JointDataRevoluteUnboundedUnalignedTpl<Scalar,Options> JointDataDerived;
    typedef JointModelRevoluteUnboundedUnalignedTpl<Scalar,Options> JointModelDerived;
    typedef ConstraintRevoluteUnalignedTpl<Scalar,Options> Constraint_t;
    typedef SE3Tpl<Scalar,Options> Transformation_t;
    typedef MotionRevoluteUnalignedTpl<Scalar,Options> Motion_t;
    typedef BiasZeroTpl<Scalar,Options> Bias_t;
    typedef Eigen::Matrix<Scalar,6,NV,Options> F_t;

    // [ABA]
    typedef Eigen::Matrix<Scalar,6,NV,Options> U_t;
    typedef Eigen::Matrix<Scalar,NV,NV,Options> D_t;
    typedef Eigen::Matrix<Scalar,6,NV,Options> UD_t;

    PINOCCHIO_JOINT_DATA_BASE_ACCESSOR_DEFAULT_RETURN_TYPE
  };

  template<typename Scalar, int Options>
  struct traits< JointDataRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  { typedef JointRevoluteUnboundedUnalignedTpl<Scalar,Options> JointDerived; };

  template<typename Scalar, int Options>
  struct traits <JointModelRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  { typedef JointRevoluteUnboundedUnalignedTpl<Scalar,Options> JointDerived; };

  template<typename _Scalar, int _Options>
  struct JointDataRevoluteUnboundedUnalignedTpl
  : public JointDataBase< JointDataRevoluteUnboundedUnalignedTpl<_Scalar,_Options> >
  {
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    typedef JointRevoluteUnboundedUnalignedTpl<_Scalar,_Options> JointDerived;
    PINOCCHIO_JOINT_DATA_TYPEDEF_TEMPLATE(JointDerived);
    PINOCCHIO_JOINT_DATA_BASE_DEFAULT_ACCESSOR

    Transformation_t M;
    Constraint_t S;
    Motion_t v;
    Bias_t c;

    // [ABA] specific data
    U_t U;
    D_t Dinv;
    UD_t UDinv;

    JointDataRevoluteUnboundedUnalignedTpl() {}

    template<typename Vector3Like>
    JointDataRevoluteUnboundedUnalignedTpl(const Eigen::MatrixBase<Vector3Like> & axis)
    : M(1)
    , S(axis)
    , v(axis,(Scalar)NAN)
    , U(), Dinv(), UDinv()
    {}

    static std::string classname() { return std::string("JointDataRevoluteUnboundedUnalignedTpl"); }
    std::string shortname() const { return classname(); }

  }; // struct JointDataRevoluteUnboundedUnalignedTpl

  PINOCCHIO_JOINT_CAST_TYPE_SPECIALIZATION(JointModelRevoluteUnboundedUnalignedTpl);

  template<typename _Scalar, int _Options>
  struct JointModelRevoluteUnboundedUnalignedTpl
  : public JointModelBase< JointModelRevoluteUnboundedUnalignedTpl<_Scalar,_Options> >
  {
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    typedef JointRevoluteUnboundedUnalignedTpl<_Scalar,_Options> JointDerived;
    PINOCCHIO_JOINT_TYPEDEF_TEMPLATE(JointDerived);
    typedef Eigen::Matrix<Scalar,3,1,Options> Vector3;

    typedef JointModelBase<JointModelRevoluteUnboundedUnalignedTpl> Base;
    using Base::id;
    using Base::idx_q;
    using Base::idx_v;
    using Base::setIndexes;

    JointModelRevoluteUnboundedUnalignedTpl() {}

    JointModelRevoluteUnboundedUnalignedTpl(const Scalar & x,
                                            const Scalar & y,
                                            const Scalar & z)
    : axis(x,y,z)
    {
      axis.normalize();
      assert(isUnitary(axis) && "Rotation axis is not unitary");
    }

    template<typename Vector3Like>
    JointModelRevoluteUnboundedUnalignedTpl(const Eigen::MatrixBase<Vector3Like> & axis)
    : axis(axis)
    {
      EIGEN_STATIC_ASSERT_VECTOR_ONLY(Vector3Like);
      assert(isUnitary(axis) && "Rotation axis is not unitary");
    }

    JointDataDerived createData() const { return JointDataDerived(axis); }

    using Base::isEqual;
    bool isEqual(const JointModelRevoluteUnboundedUnalignedTpl & other) const
    {
      return Base::isEqual(other) && axis == other.axis;
    }

    template<typename ConfigVector>
    void calc(JointDataDerived & data,
              const typename Eigen::MatrixBase<ConfigVector> & qs) const
    {
      typedef typename ConfigVector::Scalar OtherScalar;
      typename ConfigVector::template ConstFixedSegmentReturnType<NQ>::Type
      & q = qs.template segment<NQ>(idx_q());

      const OtherScalar & ca = q(0);
      const OtherScalar & sa = q(1);

      toRotationMatrix(axis,ca,sa,data.M.rotation());
    }

    template<typename ConfigVector, typename TangentVector>
    void calc(JointDataDerived & data,
              const typename Eigen::MatrixBase<ConfigVector> & qs,
              const typename Eigen::MatrixBase<TangentVector> & vs) const
    {
      calc(data,qs.derived());
      data.v.w = (Scalar)vs[idx_v()];
    }

    template<typename Matrix6Like>
    void calc_aba(JointDataDerived & data,
                  const Eigen::MatrixBase<Matrix6Like> & I,
                  const bool update_I) const
    {
      data.U.noalias() = I.template middleCols<3>(Motion::ANGULAR) * axis;
      data.Dinv[0] = (Scalar)(1)/axis.dot(data.U.template segment<3>(Motion::ANGULAR));
      data.UDinv.noalias() = data.U * data.Dinv;

      if (update_I)
        PINOCCHIO_EIGEN_CONST_CAST(Matrix6Like,I) -= data.UDinv * data.U.transpose();
    }

    static std::string classname() { return std::string("JointModelRevoluteUnboundedUnaligned"); }
    std::string shortname() const { return classname(); }

    /// \returns An expression of *this with the Scalar type casted to NewScalar.
    template<typename NewScalar>
    JointModelRevoluteUnboundedUnalignedTpl<NewScalar,Options> cast() const
    {
      typedef JointModelRevoluteUnboundedUnalignedTpl<NewScalar,Options> ReturnType;
      ReturnType res(axis.template cast<NewScalar>());
      res.setIndexes(id(),idx_q(),idx_v());
      return res;
    }

    // data
    ///
    /// \brief 3d main axis of the joint.
    ///
    Vector3 axis;
  }; // struct JointModelRevoluteUnboundedUnalignedTpl

} //namespace pinocchio

#include <boost/type_traits.hpp>

namespace boost
{
  template<typename Scalar, int Options>
  struct has_nothrow_constructor< ::pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  : public integral_constant<bool,true> {};

  template<typename Scalar, int Options>
  struct has_nothrow_copy< ::pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  : public integral_constant<bool,true> {};

  template<typename Scalar, int Options>
  struct has_nothrow_constructor< ::pinocchio::JointDataRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  : public integral_constant<bool,true> {};

  template<typename Scalar, int Options>
  struct has_nothrow_copy< ::pinocchio::JointDataRevoluteUnboundedUnalignedTpl<Scalar,Options> >
  : public integral_constant<bool,true> {};
}

#endif // ifndef __pinocchio_joint_revolute_unbounded_unaligned_hpp__

max_stars_repo_name: andreadelprete/pinocchio

max_issues_repo_name: andreadelprete/pinocchio

max_forks_repo_name: andreadelprete/pinocchio


text:
#ifndef OSRM_EXTRACTION_TURN_HPP
#define OSRM_EXTRACTION_TURN_HPP

#include "extractor/road_classification.hpp"
#include "extractor/travel_mode.hpp"

#include <boost/numeric/conversion/cast.hpp>

#include <cstdint>

namespace osrm
{
namespace extractor
{

struct ExtractionTurnLeg
{
    ExtractionTurnLeg(bool is_restricted,
                      bool is_motorway,
                      bool is_link,
                      int number_of_lanes,
                      int highway_turn_classification,
                      int access_turn_classification,
                      int speed,
                      RoadPriorityClass::Enum priority_class,
                      bool is_incoming,
                      bool is_outgoing)
        : is_restricted(is_restricted), is_motorway(is_motorway), is_link(is_link),
          number_of_lanes(number_of_lanes),
          highway_turn_classification(highway_turn_classification),
          access_turn_classification(access_turn_classification), speed(speed),
          priority_class(priority_class), is_incoming(is_incoming), is_outgoing(is_outgoing)
    {
    }

    const bool is_restricted;
    const bool is_motorway;
    const bool is_link;
    const int number_of_lanes;
    const int highway_turn_classification;
    const int access_turn_classification;
    const int speed;
    const RoadPriorityClass::Enum priority_class;
    const bool is_incoming;
    const bool is_outgoing;
};

struct ExtractionTurn
{
    ExtractionTurn(double angle,
                   int number_of_roads,
                   bool is_u_turn,
                   bool has_traffic_light,
                   bool is_left_hand_driving,
                   bool source_restricted,
                   TravelMode source_mode,
                   bool source_is_motorway,
                   bool source_is_link,
                   int source_number_of_lanes,
                   int source_highway_turn_classification,
                   int source_access_turn_classification,
                   int source_speed,
                   RoadPriorityClass::Enum source_priority_class,

                   bool target_restricted,
                   TravelMode target_mode,
                   bool target_is_motorway,
                   bool target_is_link,
                   int target_number_of_lanes,
                   int target_highway_turn_classification,
                   int target_access_turn_classification,
                   int target_speed,
                   RoadPriorityClass::Enum target_priority_class,

                   const std::vector<ExtractionTurnLeg> &roads_on_the_right,
                   const std::vector<ExtractionTurnLeg> &roads_on_the_left)
        : angle(180. - angle), number_of_roads(number_of_roads), is_u_turn(is_u_turn),
          has_traffic_light(has_traffic_light), is_left_hand_driving(is_left_hand_driving),

          source_restricted(source_restricted), source_mode(source_mode),
          source_is_motorway(source_is_motorway), source_is_link(source_is_link),
          source_number_of_lanes(source_number_of_lanes),
          source_highway_turn_classification(source_highway_turn_classification),
          source_access_turn_classification(source_access_turn_classification),
          source_speed(source_speed), source_priority_class(source_priority_class),

          target_restricted(target_restricted), target_mode(target_mode),
          target_is_motorway(target_is_motorway), target_is_link(target_is_link),
          target_number_of_lanes(target_number_of_lanes),
          target_highway_turn_classification(target_highway_turn_classification),
          target_access_turn_classification(target_access_turn_classification),
          target_speed(target_speed), target_priority_class(target_priority_class),

          roads_on_the_right(roads_on_the_right), roads_on_the_left(roads_on_the_left), weight(0.),
          duration(0.)

    {
    }
    const double angle;
    const int number_of_roads;
    const bool is_u_turn;
    const bool has_traffic_light;
    const bool is_left_hand_driving;

    // source info
    const bool source_restricted;
    const TravelMode source_mode;
    const bool source_is_motorway;
    const bool source_is_link;
    const int source_number_of_lanes;
    const int source_highway_turn_classification;
    const int source_access_turn_classification;
    const int source_speed;
    const RoadPriorityClass::Enum source_priority_class;

    // target info
    const bool target_restricted;
    const TravelMode target_mode;
    const bool target_is_motorway;
    const bool target_is_link;
    const int target_number_of_lanes;
    const int target_highway_turn_classification;
    const int target_access_turn_classification;
    const int target_speed;
    const RoadPriorityClass::Enum target_priority_class;

    const std::vector<ExtractionTurnLeg> roads_on_the_right;
    const std::vector<ExtractionTurnLeg> roads_on_the_left;

    double weight;
    double duration;
};
}
}

#endif

max_stars_repo_name: asaveljevs/osrm-backend

max_issues_repo_name: serarca/osrm-backend

max_forks_repo_name: serarca/osrm-backend