text:
/**
@file

Copyright John Reid 2006, 2007, 2013
*/

#include "bio_test_defs.h"
#include "bio/defs.h"
#include "bio/lcs.h"
#include "bio/max_chain.h"

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/parameterized_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/assign/list_of.hpp>
#include <boost/range/algorithm/equal.hpp>

using namespace boost;
using namespace boost::assign;
using boost::unit_test::test_suite;
using namespace std;
USING_BIO_NS;

//#define VERBOSE_CHECKING

struct random_gen
{
	typedef boost::mt19937 rng;                 // produces randomness out of thin air
	typedef boost::uniform_int< int > int_dist;
	typedef boost::uniform_int< unsigned > unsigned_dist;
	typedef boost::uniform_real< double > weight_dist;
	typedef boost::variate_generator< rng &, int_dist > int_gen;
	typedef boost::variate_generator< rng &, unsigned_dist > unsigned_gen;
	typedef boost::variate_generator< rng &, weight_dist > weight_gen;

	unsigned_dist _num_points_dist;
	unsigned_gen _num_points_gen;

	unsigned_dist _char_dist;
	unsigned_gen _char_gen;

	int_dist _coord_dist;
	int_gen _coord_gen;

	weight_dist _weight_dist;
	weight_gen _weight_gen;

	random_gen(
		rng & _rng,
		unsigned max_num_points,
		int max_coord,
		unsigned alphabet_size = 26 )
		: _num_points_dist( 0, max_num_points )
		, _num_points_gen( _rng, _num_points_dist )
		, _char_dist( 0, alphabet_size - 1 )
		, _char_gen( _rng, _char_dist )
		, _coord_dist( 0, max_coord )
		, _coord_gen( _rng, _coord_dist )
		, _weight_dist( 0.0, 1.0 )
		, _weight_gen( _rng, _weight_dist )
	{
	}

	char character( )
	{
		return 'a' + _char_gen();
	}

	int coord( )
	{
		return _coord_gen();
	}

	std::pair< int, int > coord_range( )
	{
		std::pair< int, int > result( coord(), coord() );
		if( result.first > result.second )
		{
			std::swap( result.first, result.second );
		}
		return result;
	}

	int num_points( )
	{
		return _num_points_gen();
	}

	double weight( )
	{
		return _weight_gen();
	}

	template< unsigned d >
	boost::array< int, d > point()
	{
		typedef boost::array< int, d > point;
		point p;
		std::generate( p.begin(), p.end(), boost::bind< int >( &random_gen::coord, this ) );
		return p;
	}

	template< unsigned d >
	std::pair< boost::array< int, d >, boost::array< int, d > > range()
	{
		typedef std::pair< boost::array< int, d >, boost::array< int, d > > point_pair;
		point_pair points( point< d >(), point< d >() );
		//make sure first <= second in every dimension
		for( unsigned i = 0; d != i; ++i )
		{
			if( points.first[ i ] > points.second[ i ] )
			{
				std::swap( points.first[ i ], points.second[ i ] );
			}
		}
		return points;
	}

};

namespace range_tree {

template< unsigned d >
struct test
{
	struct traits : point_traits< int, d >
	{
		typedef double weight;
	};

	typedef std::pair< typename traits::coord_array, typename traits::weight > test_point;
	typedef std::vector< test_point > test_point_vec;

	typedef RangeTree< d, traits > range_tree;
	typedef typename range_tree_ns::node< d, traits, d > node;
	typedef typename node::ptr tree;

	random_gen::rng _rng; //random number generator

	test() : _rng( 42u ) { } //seed the rng

	void generate_and_check( unsigned num_test_cases, unsigned max_points, unsigned max_coord )
	{
		cout
			<< "******* lcs::range_tree::test::check(): "
			<< num_test_cases << " test cases of up to "
			<< max_points << " points of dimension " << d << "\n";

		for( unsigned i = 0; num_test_cases != i; ++i )
		{
			check( create_case( max_points, max_coord ) );
		}
	}

	test_point_vec create_case( unsigned max_points, unsigned max_coord )
	{
		random_gen gen( _rng, max_points, max_coord );

		const unsigned num_points = gen.num_points();				// generate a random number of points

		test_point_vec test_case;
		for( unsigned i = 0; num_points != i; ++i )
		{
			test_point p( gen.point< d >(), gen.weight() );
			test_case.push_back( p );
		}

		return test_case;
	}

	void print_point( typename traits::point p )
	{
		for( unsigned i = 0; d != i; ++i )
		{
			std::cout << ( *p )[ i ];
			if( d - 1 != i )
			{
				std::cout << ",";
			}
		}
	}

	void
	check( const test_point_vec & test_points )
	{
		using namespace range_tree;
		using namespace boost::lambda;
		using boost::lambda::_1;

		//get the set of all points...
		typedef std::set< typename traits::point, typename traits::point_less > point_set;
		point_set points;
		for( unsigned i = 0; test_points.size() != i; ++i )
		{
			typename traits::point p = boost::addressof( test_points[ i ].first );
			points.find( p );
			points.insert( p );
#ifdef VERBOSE_CHECKING
			print_point( p );
			std::cout << std::endl;
#endif //VERBOSE_CHECKING
		}
#ifdef VERBOSE_CHECKING
		std::cout << std::endl;
#endif //VERBOSE_CHECKING

		//build tree
		tree t;
		typename range_tree_ns::create_tree< d, traits, d >()( t, points );
		typename range_tree_ns::check_tree< d, traits, d >()( t );
#ifdef VERBOSE_CHECKING
		typename range_tree_ns::visit_tree< d, traits, d >()( t, typename range_tree::print_point< d, traits, d >() );
#endif //VERBOSE_CHECKING

		//check tree has the correct points
		point_set tree_points;
		typename range_tree_ns::visit_tree< d, traits, d >()(
			t,
			range_tree_ns::make_get_points< d, traits, d >(
				std::inserter( tree_points, tree_points.begin() ) ) );
		BOOST_ASSERT( points == tree_points );

		//insert points and weights
		BOOST_FOREACH( const typename test_point_vec::value_type & v, test_points )
		{
			typename range_tree_ns::insert_point< d, traits, d >()( t, boost::addressof( v.first ), v.second );
			//std::cout << v.first[ 0 ] << "," << v.first[ 1 ] << " - " << max_chain::max_weight< 2 >()( t, addressof( v.first ) ) << "\n";
		}
		//std::cout << "\n";

#ifdef VERBOSE_CHECKING
		//query points
		BOOST_FOREACH( traits::point p, points )
		{
			std::copy( p->begin(), p->end(), std::ostream_iterator< traits::coord >( std::cout, "," ) );
			std::cout << " - " << range_tree::max_weight< d >()( t, p ) << "\n";
		}
#endif //VERBOSE_CHECKING

		const typename range_tree::heaviest_point null_point =
			typename range_tree_ns::max_weight< d, traits, d >()( t, traits::always_dominated_point() );
		BOOST_CHECK_EQUAL( null_point.first, traits::always_dominated_point() );
		BOOST_CHECK_EQUAL( null_point.second, 0.0 );

		BOOST_FOREACH( typename traits::point p1, points )
		{
			const typename range_tree::heaviest_point h1 =
				typename range_tree_ns::max_weight< d, traits, d >()( t, p1 );

			BOOST_FOREACH( typename traits::point p2, points )
			{
				if( typename traits::point_dominate()( p1, p2 ) )
				{
					const typename range_tree::heaviest_point h2 =
						range_tree_ns::max_weight< d, traits, d >()( t, p2 );

					BOOST_CHECK( h1.second >= h2.second );
				}
			}
		}
	}
};

void
generate_and_check( )
{
	const unsigned num_tests = 10;
	const std::vector< unsigned > max_points = list_of( 4 )( 15 )( 300 );

	BOOST_FOREACH( unsigned mp, max_points )
	{
		const unsigned max_coord = mp;

		test<  1 >().generate_and_check( num_tests, mp, max_coord );
		test<  2 >().generate_and_check( num_tests, mp, max_coord );
#ifndef _MSC_VER
		test<  3 >().generate_and_check( num_tests, mp, max_coord );
		test<  4 >().generate_and_check( num_tests, mp, max_coord );

#if 0
		test<  5 >().generate_and_check( num_tests, mp, max_coord );
		test<  6 >().generate_and_check( num_tests, mp, max_coord );
		test<  7 >().generate_and_check( num_tests, mp, max_coord );
		test<  7 >().generate_and_check( num_tests, mp, max_coord );
		test<  8 >().generate_and_check( num_tests, mp, max_coord );
		test<  9 >().generate_and_check( num_tests, mp, max_coord );
#endif
#endif //_MSC_VER
	}
}

}

namespace max_chain_test {

struct V
{
	typedef char character;

	V( char c = ' ', int start = -1, int end = 0, double score = 1.0 )
		: _start( start )
		, _end( end )
		, _c( c )
		, _score( score )
	{
		if( _start >= end )
		{
			throw std::invalid_argument( "start must be < end" );
		}
	}

	int _start;
	int _end;
	char _c;
	double _score;
};

struct VTraits
{
    typedef V hit;
    typedef const V * data;
	typedef char character;
	typedef double weight;
	typedef int coord;

	static data get_data( const V & v ) { return boost::addressof( v ); }
	static character get_char( const V & v ) { return v._c; }
	static weight get_weight( const V & v ) { return v._score; }
	static coord get_start( const V & v ) { return v._start; }
	static coord get_end( const V & v ) { return v._end; }
};

struct CharExtractor
{
	char operator()( const V & v ) const
	{
		return v._c;
	}
};

struct StartExtractor
{
	int operator()( const V & v ) const
	{
		return v._start;
	}
};

struct EndExtractor
{
	int operator()( const V & v ) const
	{
		return v._end;
	}
};

struct ScoreExtractor
{
	double operator()( const V & v ) const
	{
		return v._score;
	}
};

template< unsigned d >
struct VBox
{
	struct traits : point_traits< int, d >
	{
		typedef VBox< d > box;
		typedef const box * box_ptr;
		typedef double weight;
		typedef const V * data;
		typedef std::vector< box > box_vec;					/**< Vector of boxes. */
		typedef std::list< box_ptr > box_ptr_list;			/**< List of box pointers. */
		typedef typename point_traits< int, d >::coord_array coord_array;
		typedef typename point_traits< int, d >::point point;

		static unsigned get_dimensions() { return d; }

		static box make_box( data _d, weight w, coord_array s, coord_array e ) { return box( _d, w, s, e ); }
		static point get_start( const box & b ) { return boost::addressof( b._start ); }
		static point get_end( const box & b ) { return boost::addressof( b._end ); }
	};

	typedef typename traits::box_ptr box_ptr;
	typedef typename traits::weight weight;
	typedef typename traits::data data;
	typedef typename traits::coord_array point;

	VBox(
		typename traits::data _d,
		typename traits::weight w,
		typename traits::coord_array s,
		typename traits::coord_array e )
		: _data( _d )
		, _weight( w )
		, _start( s )
		, _end( e )
	{
	}

    typename traits::data _data;
	typename traits::weight _weight;
	typename traits::coord_array _start;
	typename traits::coord_array _end;

};

typedef LCS< V, char, CharExtractor, StartExtractor, EndExtractor, ScoreExtractor > lcs;

typedef std::vector< V > v_vec;
typedef std::vector< v_vec > v_array;
typedef boost::tuple< std::string, v_array, lcs::string, double > lcs_test_case;
typedef std::vector< lcs_test_case > lcs_test_case_vec;

} // namespace max_chain_test

bool
operator==( const max_chain_test::lcs::string & lhs, const max_chain_test::lcs::string & rhs )
{
    return boost::equal( lhs, rhs );
}

std::ostream &
operator<<( std::ostream & os, const max_chain_test::lcs::string & s )
{
    os << "\"";
    BOOST_FOREACH( max_chain_test::lcs::character c, s )
    {
        os << c;
    }
    os << "\"";
    return os;
}

namespace max_chain_test {

struct generate_lcs_test_case
{
	random_gen::rng _rng;
	generate_lcs_test_case() : _rng( 42u ) { }

	template< typename Int >
	Int gen_int( Int max )
	{
		typedef uniform_int< Int > dist;
		return boost::variate_generator< random_gen::rng &, dist >( _rng, dist( 0, max ) )();
	};

	template< typename Real >
	Real gen_real( )
	{
		typedef uniform_real< Real > dist;
		return boost::variate_generator< random_gen::rng &, dist >( _rng, dist( 0.0, 1.0 ) )();
	};

	v_array operator()( unsigned num_sequences, unsigned max_sequence_length, unsigned max_coord, unsigned alphabet_size )
	{
		v_array result( num_sequences );
		for( unsigned i = 0; num_sequences != i; ++i )
		{
			const unsigned length = gen_int< unsigned >( max_sequence_length );
			result[ i ].resize( length );
			for( unsigned j = 0; length != j; ++j )
			{
				const char c = 'a' + gen_int< unsigned >( alphabet_size - 1 );
				int start = gen_int< int >( max_coord );
				int end = gen_int< int >( max_coord );
				if( start > end )
				{
					std::swap( start, end );
				}
				else if( start == end )
				{
					++end;
				}
				const double score = gen_real< double >();
				result[ i ][ j ] = V( c, start, end, score );
			}
		}
		return result;
	}
};

//using boost::assign::tuple_list_of;
const lcs_test_case_vec test_cases = tuple_list_of
    (
        "empty"
        ,
        v_array()
        ,
        lcs::string()
        ,
        0.0
    )
    (
        "singular - 2D"
        ,
        list_of< v_vec >
            ( list_of ( V( 'a',  0, 10, 1. ) ) )
            ( list_of ( V( 'a', 10, 20, 2. ) ) )
        ,
        list_of( 'a' )
        ,
        1.5
    )
    (
        "singular - 3D"
        ,
        list_of< v_vec >
            ( list_of ( V( 'a',  0, 10, 1. ) ) )
            ( list_of ( V( 'a', 10, 20, 2. ) ) )
            ( list_of ( V( 'a', 20, 30, 3. ) ) )
        ,
        list_of( 'a' )
        ,
        2.
    )
    (
        "singular - 4D"
        ,
        list_of< v_vec >
            ( list_of ( V( 'a',  0, 10, 1. ) ) )
            ( list_of ( V( 'a', 10, 20, 2. ) ) )
            ( list_of ( V( 'a', 20, 30, 3. ) ) )
            ( list_of ( V( 'a', 30, 40, 4. ) ) )
        ,
        list_of( 'a' )
        ,
        2.5
    )
	(
		"problem from random tests - max_chain swaps first two matches "
		,
		list_of< v_vec >
			( list_of
                ( V( 'd', 17, 20, 0.996334 ) )( V( 'f', 6, 14,  0.448611 ) )( V( 'h', 9, 14, 0.886196 ) )
                ( V( 'e', 7, 16,  0.391526 ) )( V( 'h', 15, 20, 0.695618 ) )( V( 'f', 5, 7, 0.619589 ) )
                ( V( 'b', 14, 15, 0.794197 ) )( V( 'd', 5, 19,  0.588202 ) )( V( 'a', 13, 20, 0.642326 ) )
                ( V( 'b', 6, 13,  0.579984 ) )( V( 'b', 0, 4,    0.56066 ) )( V( 'i', 8, 17, 0.676468 ) )
                ( V( 'h', 0, 18,  0.269821 ) )( V( 'g', 13, 17, 0.498256 ) )( V( 'd', 1, 2, 0.0585509 ) )
                ( V( 'e', 8, 20,  0.784897 ) ) )
			( list_of
                ( V( 'j', 13, 18, 0.825131 ) )( V( 'j', 15, 20, 0.373463 ) )( V( 'b', 0, 15, 0.735506 ) )
                ( V( 'i', 0, 11,  0.371422 ) )( V( 'b', 7, 8,   0.801367 ) )( V('g', 0, 9, 0.733616 ) )
                ( V( 'b', 13, 17,  0.17348 ) )( V( 'c', 3, 14,  0.496166 ) )( V( 'j', 3, 11, 0.685911 ) )
                ( V( 'j', 12, 13,  0.74722 ) )( V( 'g', 1, 13,  0.950313 ) )( V( 'd', 3, 20, 0.56555 ) )
                ( V( 'b', 3, 10,  0.165294 ) )( V( 'd', 3, 9,    0.58581 ) )( V( 'd', 10, 16, 0.697671 ) )
                ( V( 'b', 11, 13, 0.782582 ) )( V( 'g', 3, 18,  0.852728 ) )( V( 'j', 9, 20, 0.376914 ) )
                ( V( 'e', 2, 12,  0.480906 ) )( V( 'd', 15, 16, 0.937353 ) )( V( 'j', 1, 8, 0.631171 ) )
                ( V( 'f', 9, 20,  0.903913 ) )( V( 'g', 8, 15,  0.780213 ) )( V( 'e', 13, 15, 0.0572803 ) )
                ( V( 'f', 0, 18,  0.131349 ) )( V( 'h', 2, 12,  0.922623 ) )( V( 'f', 3, 4, 0.0134026 ) )
                ( V( 'b', 5, 6,  0.0445949 ) )( V( 'h', 1, 3,    0.27572 ) )( V( 'f', 10, 11, 0.787127 ) ) )
		,
		list_of( 'b' )( 'f' )( 'b' )( 'd' )
		,
		3.1396045
	)
	(
		"problem with start == end in first V"
		,
		list_of< v_vec >
			( list_of( V( 'b', 5, 6, 0.386 ) )( V( 'f', 3, 11, 0.191 ) )( V( 'h', 4, 11, 0.619 ) )( V( 'b', 12, 14, 0.484 ) )( V( 'c', 14, 18, 0.602 ) ) )
		,
		list_of( 'h' )( 'b' )( 'c' )
		,
		1.705
	)
	(
		"one good hit score"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1, 0.1 ) )( V( 'b', 1, 2, 0.1 ) )( V( 'c', 2, 3, 0.1 ) )( V( 'd', 3, 4, 0.1 ) ) )
			( list_of( V( 'd', 0, 1, 0.9 ) )( V( 'a', 0, 1, 0.1 ) )( V( 'b', 1, 2, 0.1 ) )( V( 'c', 2, 3, 0.1 ) )( V( 'd', 3, 4, 0.1 ) ) )
		,
		list_of( 'd' )
		,
		0.5
	)
	(
		"long chain score"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1, 0.1 ) )( V( 'b', 1, 2, 0.1 ) )( V( 'c', 2, 3, 0.1 ) )( V( 'd', 3, 4, 0.1 ) ) )
			( list_of( V( 'd', 0, 1, 0.3 ) )( V( 'a', 0, 1, 0.1 ) )( V( 'b', 1, 2, 0.1 ) )( V( 'c', 2, 3, 0.1 ) )( V( 'd', 3, 4, 0.1 ) ) )
		,
		list_of( 'a' )( 'b' )( 'c' )( 'd' )
		,
		0.4
	)
	(
		"next to"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( list_of( V( 'a', 0, 5 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) ) )
		,
		list_of( 'c' )( 'g' )
		,
		2.0
	)
	(
		"overlapping"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 3 ) )( V( 'c', 1, 3 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( list_of( V( 'a', 0, 9 ) )( V( 'c', 1, 9 ) )( V( 'g', 2, 9 ) ) )
		,
		list_of( 'a' )
		,
		1.0
	)
	(
		"empty sequence"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( v_vec() )
		,
		lcs::string()
		,
		0.0
	)
	(
		"empty sequence 2"
		,
		list_of< v_vec >
			( v_vec() )
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
		,
		lcs::string()
		,
		0.0
	)
	(
		"identical"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
			( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
		,
		list_of( 'a' )( 'c' )( 'g' )( 't' )
		,
		4.0
	)
	(
		"score"
		,
		list_of< v_vec >
			( list_of( V( 'a', 0, 1, 0.1 ) )( V( 'c', 1, 2, 1.0 ) ) )
			( list_of( V( 'c', 0, 1, 1.0 ) )( V( 'c', 1, 2, 0.1 ) ) )
		,
		list_of( 'c' )
		,
		1.0
	)
	;

struct int_cmp
{
	bool operator()( int lhs, int rhs ) { return lhs < rhs; }
};

typedef std::pair< lcs::string, double > lcs_result;

template< unsigned d >
lcs_result
get_lcs_using_max_chain_d( const v_array & sequences, bool print = false )
{
    typedef typename VBox< d >::traits box_traits;
    typename box_traits::box_vec boxes;
    auto values = box_generator< VTraits >::template boxes_from_sequences< box_traits >(
        sequences,
        std::back_inserter( boxes ),
        50000
    );

	typedef typename VBox< d >::traits traits;
	typedef MaxChain< d, VBox< d >, traits > max_chain;
	max_chain _max_chain( boxes );

	lcs::string _longest;
	double weight = 0.0;
	BOOST_FOREACH( typename max_chain::box_ptr b, _max_chain.maximal_chain )
	{
		_longest.push_back( b->_data->_c );
#if 0
		std::cout
			<< b->_data->_start << ","
			<< b->_data->_end << ","
			<< b->_data->_c << ","
			<< b->_data->_score << endl;
#endif
		weight += b->_weight;
	}
	BOOST_CHECK_CLOSE( _max_chain.maximal_weight, weight, 0.001 );

	return std::make_pair( _longest, _max_chain.maximal_weight );
}

lcs_result
get_lcs_using_max_chain( const v_array & sequences, bool print = false )
{
	switch( boost::size( sequences ) )
	{
	case 0: return get_lcs_using_max_chain_d< 0 >( sequences, print ); break;
	case 1: return get_lcs_using_max_chain_d< 1 >( sequences, print ); break;
	case 2: return get_lcs_using_max_chain_d< 2 >( sequences, print ); break;
	case 3: return get_lcs_using_max_chain_d< 3 >( sequences, print ); break;
	case 4: return get_lcs_using_max_chain_d< 4 >( sequences, print ); break;
	case 5: return get_lcs_using_max_chain_d< 5 >( sequences, print ); break;
	case 6: return get_lcs_using_max_chain_d< 6 >( sequences, print ); break;
	case 7: return get_lcs_using_max_chain_d< 7 >( sequences, print ); break;
	//case 8: return get_lcs_using_max_chain_d< 8 >( sequences, print ); break;
	default: throw std::logic_error( "Too many test input sequences" );
	}
}

void
check_max_chain( const lcs_test_case & test_case )
{
    using namespace boost;

	cout << "******* check_max_chain(): \"" << test_case.get< 0 >() << "\"\n";

	if( 1 != test_case.get< 1 >().size() ) //does not work with one sequence at the moment
	{
		lcs_result result = get_lcs_using_max_chain( test_case.get< 1 >() );
//        std::cout
//            << size( test_case.get< 2 >() ) << " : " << size( result.first ) << ", "
//            << test_case.get< 2 >() << " : " << result.first
//            << ", weights " << test_case.get< 3 >() << " : " << result.second << "\n";

//        if( ! boost::equal( test_case.get< 2 >(), result.first ) ) {
//            std::cout << "Not equal under boost range equality\n";
//        }
//        if( test_case.get< 2 >() != result.first) {
//            std::cout << "Not equal under != operator\n";
//        }
		BOOST_CHECK_EQUAL( test_case.get< 2 >(), result.first );
		BOOST_CHECK_CLOSE( test_case.get< 3 >(), result.second, 0.001 );
	}
}

lcs_result
get_lcs_using_lcs_algorithm( const v_array & sequences )
{
	lcs _lcs( sequences );
	_lcs.calculate_best();

	return lcs_result( _lcs.get_best().get_string(), _lcs.get_best().get_score() );
}

void
check_random_lcs_test_cases( )
{
	typedef boost::array< unsigned, 5 > test_params;
	typedef std::vector< test_params > test_params_vec;

	// # tests, # seqs, max seq length, alphabet_size, max coord
	const test_params_vec _params = list_of
		( list_of( 300 )( 0 )( 100 )( 20 )( 100 ) )
		( list_of(   1 )( 1 )( 100 )( 20 )( 100 ) )
		( list_of( 300 )( 2 )( 100 )( 20 )( 100 ) )
		( list_of( 300 )( 3 )( 100 )( 20 )( 100 ) )
		( list_of( 300 )( 4 )(  30 )( 20 )( 100 ) )
		( list_of( 300 )( 5 )(  30 )( 20 )( 100 ) )
		( list_of( 300 )( 5 )(  30 )(  4 )( 100 ) )
		( list_of( 300 )( 5 )(  30 )(  4 )(  10 ) )
		( list_of( 300 )( 6 )(  30 )(  4 )(  10 ) )
		( list_of( 300 )( 7 )(  30 )(  4 )(  10 ) )
		;

	generate_lcs_test_case gen;

	BOOST_FOREACH( const test_params & p, _params )
	{
		cout
			<< "******* check_random_lcs_test_cases(): checking " << p[ 0 ] << " random cases\n"
			<< "# seqs: " << p[ 1 ] << "\n"
			<< "max length: " << p[ 2 ] << "\n"
			<< "alphabet size: " << p[ 3 ] << "\n"
			<< "max coord: " << p[ 4 ] << "\n"
			;

		double lcs_cumulative_time = 0.0;
		double mc_cumulative_time = 0.0;
		for( unsigned i = 0; p[ 0 ] != i; ++i )
		{
			//generate a test case
			v_array sequences( gen( p[ 1 ], p[ 2 ], p[ 3 ], p[ 4 ] ) );
#ifdef VERBOSE_CHECKING
			cout << "# seqs: " << sequences.size();
#endif //VERBOSE_CHECKING

			boost::timer _timer;
			const lcs_result _lcs_result = get_lcs_using_lcs_algorithm( sequences );
			lcs_cumulative_time += _timer.elapsed();
#ifdef VERBOSE_CHECKING
			cout  << ", LCS: " << _timer.elapsed();
#endif //VERBOSE_CHECKING

			if( 1 == sequences.size() )
			{
				//max_chain does not work for single sequences...
#ifdef VERBOSE_CHECKING
				cout << "\n";
#endif //VERBOSE_CHECKING
				continue;
			}

			_timer.restart();
			const lcs_result max_chain_result = get_lcs_using_max_chain( sequences );
			mc_cumulative_time += _timer.elapsed();
#ifdef VERBOSE_CHECKING
			cout  << ", MC: " << _timer.elapsed() << "\n";
#endif //VERBOSE_CHECKING

			//do we have a mismatch?
			if( max_chain_result.first != _lcs_result.first )
			{
				get_lcs_using_max_chain( sequences, true ); //run again with print = true
				BOOST_FOREACH( const v_vec & seq, sequences )
				{
					//( list_of( V( 'a', 0, 1 ) )( V( 'c', 1, 2 ) )( V( 'g', 2, 3 ) )( V( 't', 3, 4 ) ) )
					cout << "( list_of";
					BOOST_FOREACH( const V & v, seq )
					{
						cout << "( V( '" << v._c << "', " << v._start << ", " << v._end << ", " << v._score << " ) )";
					}
					cout << " )\n";
				}
				//throw std::logic_error( "Results don't agree" );
			}

			//check
			BOOST_CHECK_EQUAL( max_chain_result.first, _lcs_result.first );
			BOOST_CHECK_CLOSE( max_chain_result.second, _lcs_result.second, 0.001 );
		}
		std::cout << "LCS:" << lcs_cumulative_time << " secs\n";
		std::cout << "MC :" << mc_cumulative_time << " secs\n";
	}
}

void
check_lcs( const lcs_test_case & test_case )
{
	cout << "******* check_lcs(): \"" << test_case.get< 0 >() << "\"\n";

	lcs _lcs( test_case.get< 1 >() );
	_lcs.calculate_best();

#ifdef VERBOSE_CHECKING
	std::cout << "Universe:\n";
	BOOST_FOREACH( lcs::character c, _lcs._universe )
	{
		std::cout << c << "\n";
	}

	std::cout << "Sequences:\n";
	BOOST_FOREACH( lcs::seq s, _lcs._sequences )
	{
		BOOST_FOREACH( lcs::value c, s )
		{
			std::cout << _lcs.get_char( c  )<< " ";
		}
		std::cout << "\n";
	}

	std::cout << "End positions:\n";
	BOOST_FOREACH( lcs::seq s, _lcs._sequences )
	{
		BOOST_FOREACH( lcs::value v, s )
		{
			std::cout << _lcs.get_end( v  )<< " ";
		}
		std::cout << "\n";
	}

	std::cout << "End position\n";
	BOOST_FOREACH( unsigned i, _lcs._end_position )
	{
		std::cout << i << " ";
	}
	std::cout << "\n";

	std::cout << "End value\n";
	BOOST_FOREACH( unsigned i, _lcs._end_value )
	{
		std::cout << i << " ";
	}
	std::cout << "\n";

	std::cout << "Indices\n";
	lcs::multi_iterator end_it = _lcs.get_end_iterator();
	for( lcs::multi_iterator end_it = _lcs.get_end_iterator(); ! end_it.at_end(); end_it.next() )
	{
		BOOST_FOREACH( unsigned i, end_it._ind )
		{
			std::cout << i << " ";
		}
		std::cout << "\n";

		for( lcs::multi_iterator value_it = _lcs.get_value_range_for_end( end_it._ind ); ! value_it.at_end(); value_it.next() )
		{
			for( unsigned d = 0; _lcs.get_num_seqs() != d; ++d )
			{
				std::cout << " " << _lcs.get_end( _lcs.get_value( d, value_it._ind[ d ] ) );
			}
			std::cout << "\n";
		}
	}
#endif

	BOOST_CHECK_EQUAL( test_case.get< 2 >(), _lcs.get_best().get_string() );
	BOOST_CHECK_CLOSE( test_case.get< 3 >(), _lcs.get_best().get_score(), 0.001 );

	lcs::value_holder::vec values = _lcs.get_best().get_values();

#ifdef VERBOSE_CHECKING
	BOOST_FOREACH( const lcs::value_holder & v, values )
	{
		std::cout
			<< v._start
			<< "," << v._end
			<< "," << v._score
			<< "," << v._c
			<< "\n";
	}
#endif //VERBOSE_CHECKING
}

}

BOOST_TEST_DONT_PRINT_LOG_VALUE( max_chain_test::lcs::string )

#ifdef TEST_AGAINST_CGAL
#include <CGAL/Cartesian.h>
#include <CGAL/Range_segment_tree_traits.h>
#include <CGAL/Range_tree_k.h>

void
check_cgal_range_tree()
{
	cout << "******* check_cgal_range_tree()\n";

	typedef CGAL::Cartesian< int > K;
	typedef CGAL::Range_segment_tree_set_traits_2< K > Traits;
	typedef CGAL::Range_tree_2< Traits > Range_tree_2_type;

	typedef Traits::Key Key;
	typedef std::vector<Key> KeyVec;
	typedef Traits::Interval Interval;
	std::vector<Key>::iterator first, last, current;

	KeyVec InputList = list_of
		( Key( 1, 1) )
		( Key( 1, 2) )
		( Key( 1, 3) )
		( Key( 2, 1) )
		( Key( 2, 2) )
		( Key( 2, 3) )
		( Key( 3, 1) )
		( Key( 3, 2) )
		( Key( 3, 3) )
		;

	//create a range tree
	Range_tree_2_type Range_tree_2( InputList.begin(), InputList.end() );

	const Key l( 1, 1 );
	const Key u( 3, 3 );
	const Interval win = Interval( l, u );
	std::cout
		<< "Window Query: lower left point: "
		<< l
		<< "\nupper right point: "
		<< u
		<< "\n";

	KeyVec OutputList;
	Range_tree_2.window_query( win, std::back_inserter(OutputList) );
	current = OutputList.begin();
	while( current != OutputList.end())
	{
		std::cout << (*current).x()<< "-" << (*current).y() << std::endl;
		current++;
	}
}

#endif //TEST_AGAINST_CGAL

void
register_lcs_tests( boost::unit_test::test_suite * test )
{
	using namespace max_chain_test;
	using namespace range_tree;

    test->add( BOOST_TEST_CASE( &check_random_lcs_test_cases ) );
	test->add( BOOST_PARAM_TEST_CASE( &check_lcs, test_cases.begin(), test_cases.end() ) );
	test->add( BOOST_PARAM_TEST_CASE( &check_max_chain, test_cases.begin(), test_cases.end() ) );
    test->add( BOOST_TEST_CASE( &range_tree::generate_and_check ) );
#ifdef TEST_AGAINST_CGAL
	test->add( BOOST_TEST_CASE( &check_cgal_range_tree ) );
#endif //TEST_AGAINST_CGAL
}

//
// Are we going to compile this test into its own executable?
//
#ifdef BIO_STANDALONE_TEST
BIO_DEFINE_STANDALONE_TEST( "lcs", register_lcs_tests )
#endif //BIO_STANDALONE_TEST

ext: cpp

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing,
// software distributed under the License is distributed on an
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.  See the License for the
// specific language governing permissions and limitations
// under the License.

#include <memory>
#include <unordered_map>
#include <boost/thread/mutex.hpp>
#include <boost/shared_ptr.hpp>
#include "arrow/buffer.h"
#include "dataset/dataset_cache_manager.h"
#include "runtime/worker_exec_env.h"
#include "parquet/parquet_reader.h"
#include "cache/cache_memory_pool.h"
#include "common/logging.h"
#include "cache/lru_cache.h"
#include "rpc/file_batch_reader.h"
#include "ipc/malloc.h"
#include "rpc/types.h"

DEFINE_int64(chunk_size, 2048, "The maximum chunk size of record batches");
DECLARE_bool(cache_format_arrow);

namespace pegasus {

DatasetCacheManager::DatasetCacheManager()
  : cache_engine_manager_(nullptr),
    cache_block_manager_(nullptr) {
}

DatasetCacheManager::~DatasetCacheManager() {
  if (cache_engine_manager_ != nullptr) {
    delete cache_engine_manager_;
    cache_engine_manager_ = nullptr;
  }

  if (cache_block_manager_ != nullptr) {
    delete cache_block_manager_;
    cache_block_manager_ = nullptr;
  }
}

Status DatasetCacheManager::Init() {
  ExecEnv* env =  ExecEnv::GetInstance();
  storage_factory_ = env->get_storage_factory();

  cache_block_manager_ = new DatasetCacheBlockManager();
  cache_engine_manager_ = new DatasetCacheEngineManager();

  cache_metrics_.ResetCacheMetrics();

  RETURN_IF_ERROR(cache_block_manager_->Init());
  RETURN_IF_ERROR(cache_engine_manager_->Init());

  return Status::OK();
}

CacheEngine::CachePolicy DatasetCacheManager::GetCachePolicy(
  RequestIdentity* request_identity) {
  // TODO Choose the CachePolicy based on the data type in Identity
  return CacheEngine::CachePolicy::LRU;
}

Status DatasetCacheManager::GetPartition(
  RequestIdentity* request_identity,
  std::shared_ptr<CachedPartition>* new_partition) {
    std::shared_ptr<CachedDataset> dataset;
    RETURN_IF_ERROR(cache_block_manager_->GetCachedDataSet(
      request_identity->dataset_path(), &dataset));

    RETURN_IF_ERROR(dataset->GetCachedPartition(
     request_identity->partition_path(), new_partition));
    return Status::OK();
}

Status DatasetCacheManager::WrapDatasetStream(
  RequestIdentity* request_identity,
  unordered_map<int, std::shared_ptr<CachedColumn>> request_columns,
  std::unique_ptr<rpc::FlightDataStream>* data_stream) {
  LOG(WARNING) << "Wrap the dataset into flight data stream";

  std::vector<std::shared_ptr<ChunkedArray>> chunked_arrays;
  std::vector<std::shared_ptr<CachedColumn>> columns;

  std::vector<int> col_ids = request_identity->column_indices();
  for(int index : col_ids) {
    auto iter = request_columns.find(index);
    if (iter != request_columns.end()) {
      std::shared_ptr<CachedColumn> cache_column = iter->second;
      CacheRegion* cache_region = cache_column->GetCacheRegion();
      std::shared_ptr<ChunkedArray> chunked_array = cache_region->chunked_array();
      chunked_arrays.push_back(chunked_array);
      columns.push_back(cache_column);
    } else {
      return Status::ObjectNotFound("can't find the cached column.");
    }
  }

  std::shared_ptr<arrow::Schema> schema;
  request_identity->get_schema(&schema);

  if (FLAGS_cache_format_arrow) {
    // arrow data format
    std::shared_ptr<arrow::Table> table = arrow::Table::Make(schema, chunked_arrays);

    std::shared_ptr<arrow::TableBatchReader> reader =
      std::make_shared<arrow::TableBatchReader>(*table);
    reader->set_chunksize(FLAGS_chunk_size);
    *data_stream = std::unique_ptr<rpc::FlightDataStream>(
      new rpc::TableRecordBatchStream(reader, columns, table));
 } else {
   // file data format
   //TO DO: passed the cached column data to the reader
   std::shared_ptr<rpc::CachedFileBatchReader> reader =
    std::make_shared<rpc::CachedFileBatchReader>(columns, schema);
   *data_stream = std::unique_ptr<rpc::FlightDataStream>(
      new rpc::FileBatchStream(reader, columns));
 }
 LOG(WARNING) << "Successfully wrap the dataset into flight data stream";
 return Status::OK();
}

Status DatasetCacheManager::GetDatasetStreamWithMissedColumns(
  RequestIdentity* request_identity,
  std::vector<int> col_ids,
  unordered_map<int, std::shared_ptr<CachedColumn>> cached_columns,
  std::shared_ptr<CacheEngine> cache_engine,
  std::unique_ptr<rpc::FlightDataStream>* data_stream) {

    unordered_map<int, std::shared_ptr<CachedColumn>> retrieved_columns;
    RETURN_IF_ERROR(RetrieveColumns(request_identity, col_ids, cache_engine, retrieved_columns));

    for (auto iter = cached_columns.begin(); iter != cached_columns.end(); iter ++) {
      int col_id = iter->first;
      std::shared_ptr<CachedColumn> cached_column = iter->second;
      retrieved_columns.insert(std::make_pair(col_id, cached_column));
    }

    if (data_stream == nullptr) {
      return Status::OK();
    } else {
      return WrapDatasetStream(request_identity, retrieved_columns, data_stream);
    }
}

Status DatasetCacheManager::RetrieveColumns(
  RequestIdentity* request_identity,
  const std::vector<int>& col_ids,
  std::shared_ptr<CacheEngine> cache_engine,
  unordered_map<int, std::shared_ptr<CachedColumn>>& retrieved_columns) {
    LOG(WARNING) << "Retrieve the columns from storage and insert the "
     << "retrieved columns into dataset block manager and cache engine";
    std::string dataset_path = request_identity->dataset_path();
    std::string partition_path = request_identity->partition_path();

    // Read the columns into ChunkArray.
    // Asumming the cache memory pool is only in same store.
    std::shared_ptr<CacheMemoryPool> memory_pool(new CacheMemoryPool(cache_engine));
    RETURN_IF_ERROR(memory_pool->Create());

    std::unique_ptr<ParquetReader> parquet_reader;
    std::unique_ptr<ParquetRawDataReader> parquet_raw_reader;

    // Get the ReadableFile
    std::shared_ptr<Storage> storage;
    RETURN_IF_ERROR(storage_factory_->GetStorage(partition_path, &storage));
    if (storage->GetStorageType() == Storage::HDFS) {
      std::shared_ptr<HdfsReadableFile> file;
      RETURN_IF_ERROR(std::dynamic_pointer_cast<HDFSStorage>(storage)
          ->GetReadableFile(partition_path, &file));
      if (FLAGS_cache_format_arrow) {
        parquet::ArrowReaderProperties properties(parquet::default_arrow_reader_properties());
        parquet_reader = std::unique_ptr<ParquetReader>(
          new ParquetReader(file, memory_pool.get(), properties));
      } else {
        parquet::ReaderProperties properties(memory_pool.get());

        // only when the flag of buffered_stream_enabled_ is enabled, the buffer can be allocated from the memory pool.
        properties.enable_buffered_stream();
        parquet_raw_reader = std::unique_ptr<ParquetRawDataReader>(
          new ParquetRawDataReader(file, properties));
      }
    }

    std::shared_ptr<CachedPartition> partition;
    GetPartition(request_identity, &partition);

    int64_t occupied_size = 0;
    for(auto iter = col_ids.begin(); iter != col_ids.end(); iter ++) {
      int colId = *iter;

      unordered_map<int, std::shared_ptr<BufferEntry>> object_buffers;
      std::shared_ptr<arrow::ChunkedArray> chunked_array;
      unordered_map<int, std::shared_ptr<ObjectEntry>> object_entries;
      int64_t row_counts_per_rowgroup = 0;
      if (FLAGS_cache_format_arrow) {

        LOG(INFO) << "Begin read the column chunk with col ID " << colId << " partition path " << partition_path;
        RETURN_IF_ERROR(parquet_reader->ReadColumnChunk(*iter, &chunked_array));
      } else {

        std::shared_ptr<parquet::RowGroupReader> row_group_reader;
        int row_group_counts = parquet_raw_reader->RowGroupsNum();

        for(int i = 0; i < row_group_counts; i ++) {
          std::shared_ptr<Buffer> buffer;
          LOG(INFO) << "Begin read the raw column chunk with row group ID " << i << " col ID " << colId << " partition path " << partition_path;
          RETURN_IF_ERROR(parquet_raw_reader->GetColumnBuffer(i, colId, &buffer));

          RETURN_IF_ERROR(parquet_raw_reader->GetRowGroupReader(i, &row_group_reader));
          row_counts_per_rowgroup = row_group_reader->metadata()->num_rows();

          std::shared_ptr<BufferEntry> buffer_entry =
            std::shared_ptr<BufferEntry>(new BufferEntry(buffer, row_counts_per_rowgroup));
          object_buffers[i] = std::move(buffer_entry);

          int fd = -1;
          int64_t map_size = 0;
          ptrdiff_t offset = 0;
          uint8_t* pointer = const_cast< uint8_t*>(buffer->data());

          GetMallocMapinfo(pointer, &fd, &map_size, &offset);

          std::shared_ptr<ObjectEntry> entry = std::shared_ptr<ObjectEntry>(new ObjectEntry(fd,
           offset, map_size, row_counts_per_rowgroup, buffer->size()));

          object_entries[i] = std::move(entry);
        }
      }

      int64_t column_size = memory_pool->bytes_allocated() - occupied_size;
      occupied_size += column_size;

      CacheRegion* cache_region = new CacheRegion(memory_pool,
        chunked_array, column_size, object_buffers, object_entries, row_counts_per_rowgroup);
      std::shared_ptr<CachedColumn> column = std::shared_ptr<CachedColumn>(
        new CachedColumn(partition_path, colId, cache_region));

      cache_metrics_.cached_size += column_size;
      // Insert the column if already inserted,
      // we do not put the column into the LRU cache.
      // And because this column is shared ptr, so out this for clause ,
      // it will be delete automatically.
      std::shared_ptr<CachedColumn> cached_column;
      bool is_inserted = partition->InsertColumn(colId, column, &cached_column);
      if (is_inserted && cached_column == nullptr) {
        LRUCache::CacheKey* key = new LRUCache::CacheKey(dataset_path, partition_path, colId);
        LOG(WARNING) << "Put the cached column into cache engine";

        RETURN_IF_ERROR(cache_engine->PutValue(key, column_size));

        retrieved_columns.insert(std::make_pair(colId, column));
      } else {
        retrieved_columns.insert(std::make_pair(colId, cached_column));
      }
    }
    LOG(INFO) << "the cached size is " << cache_metrics_.cached_size;
    return Status::OK();
}

std::vector<int> DatasetCacheManager::GetMissedColumnsIds(std::vector<int> col_ids,
  unordered_map<int, std::shared_ptr<CachedColumn>> cached_columns) {
   LOG(WARNING) << "Get the missed column IDs ";
   std::vector<int> missed_col_ids;
    for(auto iter = col_ids.begin(); iter != col_ids.end(); iter ++) {
        auto entry = cached_columns.find(*iter);
        if (entry == cached_columns.end()) {
            missed_col_ids.push_back(*iter);
        }
    }
    return missed_col_ids;
}

Status DatasetCacheManager::DropCachedDataset(
  std::vector<rpc::PartitionDropList> drop_lists) {
  for (auto iter = drop_lists.begin(); iter != drop_lists.end(); iter ++) {
     std::string dataset_path = iter->get_dataset_path();
     std::vector<std::string> partitions = iter->get_partitions();
     std::shared_ptr<CachedDataset> cached_dataset;
     cache_block_manager_->GetCachedDataSet(dataset_path, &cached_dataset);
     for(auto partition = partitions.begin(); partition != partitions.end(); partition ++) {
       std::string partition_path = *partition;

       // Get the cache store and delete the columns cached in lru cache when delete the partition
       std::shared_ptr<CacheEngine> cache_engine;
       CacheEngine::CachePolicy cache_policy = CacheEngine::CachePolicy::LRU; // TODO get the cache policy based on the data set type.
       RETURN_IF_ERROR(cache_engine_manager_->GetCacheEngine(cache_policy, &cache_engine));
       cached_dataset->DeletePartition(partition_path, cache_engine);
       LOG(INFO) << "Drop the cached dataset: " << dataset_path << " the partition path: " << partition_path;
     }
  }
  return Status::OK();
}

// Wrap the data to flight data stream.
// According DatasetCacheBlockManager to chech whether sotred.
// If yes, read it based on the address stored in DatasetCacheBlockManager and then return.
// If not, get dataset from hdfs and then put the dataset into CacheEngine.
//         1. Choose the CachePolicy based on the Identity.
//         2. Call DatasetCacheEngineManager#GetCacheEngine method to get CacheEngine;
Status DatasetCacheManager::GetDatasetStream(
  RequestIdentity* request_identity,
  std::unique_ptr<rpc::FlightDataStream>* data_stream) {
  // Get cache engine.
  std::shared_ptr<CacheEngine> cache_engine;
  CacheEngine::CachePolicy cache_policy = GetCachePolicy(request_identity);
  RETURN_IF_ERROR(cache_engine_manager_->GetCacheEngine(cache_policy, &cache_engine));

  cache_metrics_.total_cacherd_cnt++;
  // Check whether the dataset is cached.
  std::vector<int> col_ids = request_identity->column_indices();
  unordered_map<int, std::shared_ptr<CachedColumn>> cached_columns;

  std::shared_ptr<CachedDataset> dataset;
  cache_block_manager_->GetCachedDataSet(request_identity->dataset_path(), &dataset);
  if (dataset->GetCachedPartitions().size() == 0) {
    LOG(WARNING) << "The dataset "<< request_identity->dataset_path()
    <<" is new added. We will get all the columns from storage and"
     << " then insert the column into dataset cache block manager";

    return GetDatasetStreamWithMissedColumns(request_identity,
     col_ids, cached_columns, cache_engine, data_stream);
  } else {
    // dataset is cached
    cache_metrics_.ds_cacherd_cnt++;
    std::shared_ptr<CachedPartition> partition;
    dataset->GetCachedPartition(
     request_identity->partition_path(), &partition);
    if (partition->GetCachedColumns().size() == 0) {
      LOG(WARNING) << "The partition "<< request_identity->partition_path()
      <<" is new added. We will get all the columns from storage and"
       << "then insert the column into dataset cache block manager";
      return GetDatasetStreamWithMissedColumns(request_identity, col_ids,
       cached_columns, cache_engine, data_stream);
    } else {
      // partition is cached.
      cache_metrics_.pt_cacherd_cnt++;
      // Check which column is cached.
      partition->GetCachedColumns(request_identity->column_indices(),
       cache_engine, &cached_columns);
      if (col_ids.size() == cached_columns.size()) {

        cache_metrics_.col_cacherd_cnt++;

        if (data_stream == nullptr) {
          return Status::OK();
        } else {
          LOG(WARNING) << "All the columns are cached. And we will wrap the columns into Flight data stream";

          return WrapDatasetStream(request_identity, cached_columns, data_stream);
        }

      } else {
        // Not all columns cached.
        // Get the not cached col_ids.
        std::vector<int> missed_col_ids = GetMissedColumnsIds(col_ids, cached_columns);
        LOG(WARNING) << "Partial columns is cached and we will get the missed columns from storage";
        return GetDatasetStreamWithMissedColumns(request_identity, missed_col_ids,
         cached_columns, cache_engine, data_stream);
      }
   }
  }
}

Status DatasetCacheManager::GetLocalData(
  RequestIdentity* request_identity, std::unique_ptr<rpc::LocalPartitionInfo>* result) {
  // get the missed columns from hdfs and put the columns into cache and block manager.
  //  Then all the request columns are cached.
  GetDatasetStream(request_identity, nullptr);

  // get the columns and put the columns into in_used_columns_.
  std::shared_ptr<CachedDataset> dataset;
  cache_block_manager_->GetCachedDataSet(request_identity->dataset_path(), &dataset);

  std::shared_ptr<CachedPartition> partition;
  dataset->GetCachedPartition(
     request_identity->partition_path(), &partition);

  std::vector<int> col_ids = request_identity->column_indices();
  unordered_map<int, std::shared_ptr<CachedColumn>> cached_columns;
  // here no need to touch the value, because it already done when
  // call GetDatasetStream method. So we set the cache_engine to nullptr to skip.
  partition->GetCachedColumns(request_identity->column_indices(),
       nullptr, &cached_columns);

  std::string dataset_path = request_identity->dataset_path();
  std::string partition_path = request_identity->partition_path();

  // store the requested columns into in_used_columns_.
  for (auto iter = col_ids.begin(); iter != col_ids.end(); iter++) {
    int column_id = *iter;
    std::string key = dataset_path.append(partition_path).append(to_string(column_id));

    auto entry = cached_columns.find(column_id);
    // assert(entry != cached_columns.end());

    std::shared_ptr<CachedColumn> column = entry->second;
    {
      boost::lock_guard<boost::mutex> l(in_used_columns_lock_);
      auto in_use_entry = in_used_columns_.find(key);

      if (in_use_entry == in_used_columns_.end()) {
        // store the key firstly
        std::vector<std::shared_ptr<CachedColumn>> columns;
        columns.push_back(column);
        in_used_columns_[key] = columns;
      } else {
        // directly put the value into the in_used_columns_.
        in_use_entry->second.push_back(column);
      }
   }
  }

  // wrap LocalPartitionInfo
  std::vector<rpc::LocalColumnInfo> columns;
  for (auto iter = col_ids.begin(); iter != col_ids.end(); iter++) {
    int column_id = *iter;
    auto entry = cached_columns.find(column_id);
    // assert(entry != cached_columns.end());

    std::shared_ptr<CachedColumn> column = entry->second;
    CacheRegion* region = column->GetCacheRegion();
    unordered_map<int, std::shared_ptr<ObjectEntry>> object_entries = region->object_entries();

    int row_group_counts = object_entries.size();
    std::vector<rpc::LocalColumnChunkInfo> chunks;

    for (int i =0; i < row_group_counts; i++) {
      auto entry = object_entries.find(i);
     // assert(entry != cached_columns.end());
      std::shared_ptr<ObjectEntry> object_entry = entry->second;

      rpc::LocalColumnChunkInfo chunk;
      chunk.chunk_index = i;
      chunk.data_offset = object_entry.get()->offset_;
      chunk.data_size = object_entry.get()->data_size_;
      chunk.mmap_fd = object_entry.get()->fd_;
      chunk.mmap_size = object_entry.get()->map_size_;
      chunk.row_counts = object_entry.get()->row_counts_;

      chunks.push_back(chunk);
    }

    rpc::LocalColumnInfo column_entry;
    column_entry.column_index = column_id;
    column_entry.chunks = chunks;

    columns.push_back(column_entry);
  }

  *result = std::unique_ptr<rpc::LocalPartitionInfo>(new rpc::LocalPartitionInfo(columns));
  return Status::OK();
}

Status DatasetCacheManager::ReleaseLocalData(
  RequestIdentity* request_identity, std::unique_ptr<rpc::LocalReleaseResult>* result) {
   // get the columns and release the columns in in_used_columns_.
  std::shared_ptr<CachedDataset> dataset;
  cache_block_manager_->GetCachedDataSet(request_identity->dataset_path(), &dataset);

  std::shared_ptr<CachedPartition> partition;
  dataset->GetCachedPartition(
     request_identity->partition_path(), &partition);

  std::vector<int> col_ids = request_identity->column_indices();
  unordered_map<int, std::shared_ptr<CachedColumn>> cached_columns;
  // here no need to touch the value, because it already done when
  // call GetDatasetStream method. So we set the cache_engine to nullptr to skip.
  partition->GetCachedColumns(request_identity->column_indices(),
       nullptr, &cached_columns);

  std::string dataset_path = request_identity->dataset_path();
  std::string partition_path = request_identity->partition_path();

  for (auto iter = col_ids.begin(); iter != col_ids.end(); iter++) {
    int column_id = *iter;
    std::string key = dataset_path.append(partition_path).append(to_string(column_id));

    auto entry = cached_columns.find(column_id);
    assert(entry != cached_columns.end());

    std::shared_ptr<CachedColumn> column = entry->second;
    {
      boost::lock_guard<boost::mutex> l(in_used_columns_lock_);
      auto in_use_entry = in_used_columns_.find(key);

      if (in_use_entry == in_used_columns_.end()) {
        LOG(WARNING) << "The released key is not in the global map of in_used_columns_. Please check the key";

        return Status::KeyError("the key of ", key, "is not valid");
      } else {
        // directly delete the last value in in_used_columns_.
        in_use_entry->second.pop_back();
      }
    }
  }

  rpc::LocalReleaseResult* result_tmp = new rpc::LocalReleaseResult();
  result_tmp->result_code = 0;
  *result = std::unique_ptr<rpc::LocalReleaseResult>(result_tmp);

  return Status::OK();
}

} // namespace pegasus

ext: cc

max_stars_repo_name: jerrychenhf/pegasus

max_issues_repo_name: jerrychenhf/pegasus

max_forks_repo_name: jerrychenhf/pegasus


text:
#define BOOST_TEST_DYN_LINK

#include <iostream>
#include <boost/test/unit_test.hpp>
#include <string>
#include "cutvm_tokenizer.hpp"
#include "generator_methods.hpp"
#include "cutvm_context.hpp"
#include "cutvm_generator.hpp"
#include "value_methods.hpp"
#include "context_methods.hpp"

using namespace cuttle;

struct cutvm_generator_init_fixture {
    context_t context;
    generator_state_t state;
    tokenizer_config_t tokenizer_config;
    generator_config_t generator_config;

    void setup() {
        initialize(context);
        initialize(state);
        lang::get_cutvm_context(context);
        lang::get_cutvm_tokenizer_config(tokenizer_config);
        lang::get_cutvm_generator_config(generator_config);
    }
};

BOOST_FIXTURE_TEST_SUITE(generates_basic_cutvm_code_suite, cutvm_generator_init_fixture)

    BOOST_AUTO_TEST_CASE(case1) {
        values_t values = {
                {
                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"+", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"1", value_type::number},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"2", value_type::number},

                        {"c", value_type::func_name},
                        {"3", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name}
                }};
        call_tree_t tree = {
                {
                        {1, 4, 7, 10},
                        {2}, {3}, {},
                        {5}, {6}, {},
                        {8}, {9}, {},
                        {11, 12, 13}, {}, {}, {}, {0}
                },
                {}};

        generate(tokenizer_config, generator_config, context, values, tree, state);
        BOOST_CHECK_EQUAL(state.output, "b s +\n"
                                        "b i 1\n"
                                        "b i 2\n"
                                        "c 3 0 array");
    }

    BOOST_AUTO_TEST_CASE(case2) {
        values_t values = {
                {
                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"+", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"1", value_type::number},

                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"foo", value_type::string},
                        {"c", value_type::func_name},
                        {"1", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name},

                        {"c", value_type::func_name},
                        {"3", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name}
                }};
        call_tree_t tree = {
                {
                        {1, 4, 7, 15},
                        {2}, {3}, {},
                        {5}, {6}, {},

                        {8, 11},
                        {9}, {10}, {},
                        {12, 13, 14}, {}, {}, {},

                        {16, 17, 18}, {}, {}, {},
                        {0}
                },
                {}};

        generate(tokenizer_config, generator_config, context, values, tree, state);
        BOOST_CHECK_EQUAL(state.output, "b s +\n"
                                        "b i 1\n"
                                        "b s foo\n"
                                        "c 1 0 array\n"
                                        "c 3 0 array");
    }

    BOOST_AUTO_TEST_CASE(case3) {
        values_t values = {
                {
                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"+", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"1", value_type::number},

                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"foo", value_type::string},

                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"-", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"1", value_type::number},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"2", value_type::number},
                        {"c", value_type::func_name},
                        {"3", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name},

                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"3", value_type::number},

                        {"c", value_type::func_name},
                        {"3", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name},

                        {"c", value_type::func_name},
                        {"3", value_type::number},
                        {"0", value_type::number},
                        {"array", value_type::func_name}
                }};
        call_tree_t tree = {
                {
                        {1, 4, 7, 32},
                        {2}, {3}, {},
                        {5}, {6}, {},

                        {8, 11, 25, 28},
                        {9}, {10}, {},

                        {12, 15, 18, 21},
                        {13}, {14}, {},
                        {16}, {17}, {},
                        {19}, {20}, {},
                        {22, 23, 24}, {}, {}, {},

                        {26}, {27}, {},
                        {29, 30, 31}, {}, {}, {},

                        {33, 34, 35}, {}, {}, {},
                        {0}
                },
                {}};

        generate(tokenizer_config, generator_config, context, values, tree, state);
        BOOST_CHECK_EQUAL(state.output, "b s +\n"
                                        "b i 1\n"
                                        "b s foo\n"
                                        "b s -\n"
                                        "b i 1\n"
                                        "b i 2\n"
                                        "c 3 0 array\n"
                                        "b i 3\n"
                                        "c 3 0 array\n"
                                        "c 3 0 array");
    }

    BOOST_AUTO_TEST_CASE(case4) {
        values_t values = {
                {
                        {"call", value_type::func_name},
                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"+", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"1", value_type::number},
                        {"call", value_type::func_name},
                        {"c", value_type::func_name}, {"3", value_type::number}, {"0", value_type::number}, {"array", value_type::func_name},

                        {"b", value_type::func_name}, {"s", value_type::func_name}, {"-", value_type::string},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"2", value_type::number},
                        {"b", value_type::func_name}, {"i", value_type::func_name}, {"2", value_type::number},
                        {"c", value_type::func_name}, {"3", value_type::number}, {"0", value_type::number}, {"array", value_type::func_name}
                }};
        call_tree_t tree = {
                {
                        {1, 4, 7, 8},
                        {2}, {3}, {},
                        {5}, {6}, {},
                        {12, 15, 18, 21},
                        {9, 10, 11}, {}, {}, {},

                        {13}, {14}, {},
                        {16}, {17}, {},
                        {19}, {20}, {},
                        {22, 23, 24}, {}, {}, {},
                        {0}
                },
                {}};

        generate(tokenizer_config, generator_config, context, values, tree, state);
        BOOST_CHECK_EQUAL(state.output, "b s +\n"
                                        "b i 1\n"
                                        "b s -\n"
                                        "b i 2\n"
                                        "b i 2\n"
                                        "c 3 0 array\n"
                                        "c 3 0 array");
    }

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: cuttle-system/cuttle-lang

max_issues_repo_name: cuttle-system/cuttle-lang

max_forks_repo_name: cuttle-system/cuttle-lang


text:
/* --------------------------------------------------------------------------
*
* (C) Copyright …
*
* --------------------------------------------------------------------------
*/

/*!
 * @file MatToTransform3DConverter.hpp
 * @date 20/02/2018
 * @author Alessandro Bianco
 */

/*!
 * @addtogroup Converters
 *
 *  This is the class for type conversion from cv Matrix to ASN Transform.
 *
 *
 * @{
 */

#ifndef MAT_TO_TRANSFORM_3D_CONVERTER_HPP
#define MAT_TO_TRANSFORM_3D_CONVERTER_HPP

/* --------------------------------------------------------------------------
 *
 * Includes
 *
 * --------------------------------------------------------------------------
 */
#include <Types/CPP/Pose.hpp>
#include <Eigen/Core>
#include <opencv2/core.hpp>

namespace Converters {

/* --------------------------------------------------------------------------
 *
 * Class definition
 *
 * --------------------------------------------------------------------------
 */
class MatToTransform3DConverter
	{
	/* --------------------------------------------------------------------
	 * Public
	 * --------------------------------------------------------------------
	 */
	public:
		virtual const PoseWrapper::Transform3DConstPtr Convert(const cv::Mat transform);
		const PoseWrapper::Transform3DSharedConstPtr ConvertShared(const cv::Mat transform);

	/* --------------------------------------------------------------------
	 * Protected
	 * --------------------------------------------------------------------
	 */
        protected:

	/* --------------------------------------------------------------------
	 * Private
	 * --------------------------------------------------------------------
	 */
	private:
	};

typedef MatToTransform3DConverter MatToPose3DConverter;

}

#endif

/* MatToTransform3DConverter.hpp */
/** @} */

ext: hpp


text:
/*
 * Copyright (c) 2020 Haresh Bhachandani
 *
 * Permission is hereby granted, free of charge, to any person
 * obtaining a copy of this software and associated documentation
 * files (the "Software"), to deal in the Software without
 * restriction, including without limitation the rights to use,
 * copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following
 * conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 * OTHER DEALINGS IN THE SOFTWARE.
 */

/* Entry point and facade controller of the translator
 */
#include "Translator.h"
#include "Parser.h"

#define BOOST_FILESYSTEM_NO_DEPRECATED
#include <boost/filesystem.hpp>

#include <iostream>
#include <fstream>
#include <exception>

using namespace std;
using namespace boost::filesystem;

const string OUTPUT_PREFIX("asm");

Translator::Translator(vector<string> arguments)
{
    if(arguments.size() != 2)
        throw runtime_error("Error: Program only takes 1 argument");

    inputFilename = arguments[1];
}

void Translator::translate(string inputFilename, CodeWriter &writer)
{
    // prepare file for parsing
    std::ifstream input(inputFilename);
    Parser parser(input);

    writer.setFilename(path(inputFilename).filename().string());

    try
    {
        // while there is more to parse
        while(parser.hasMoreCommands())
        {
            // parse it
            parser.advance();

            CommandType type = parser.getCommandType();

            writer.writeAnnotation(type, parser.getArg1(), parser.getArg2());

            // write the correct command
            switch(type)
            {
                case CommandType::Label:
                {
                    writer.writeLabel(parser.getArg1());
                    break;
                }

                case CommandType::Goto:
                {
                    writer.writeGoto(parser.getArg1());
                    break;
                }

                case CommandType::If:
                {
                    writer.writeIf(parser.getArg1());
                    break;
                }

                case CommandType::Call:
                {
                    int args = stoi(parser.getArg2());
                    writer.writeCall(parser.getArg1(), args);
                    break;
                }

                case CommandType::Function:
                {
                    int args = stoi(parser.getArg2());
                    writer.writeFunction(parser.getArg1(), args);
                    break;
                }

                case CommandType::Return:
                {
                    writer.writeReturn();
                    break;
                }

                case CommandType::Arithmetic:
                {
                    writer.writeArithmetic(parser.getArg1());
                    break;
                }

                case CommandType::Push:
                case CommandType::Pop:
                {
                    try
                    {
                        int index = stoi(parser.getArg2());
                        writer.writePushPop(type, parser.getArg1(), index);
                    }
                    catch(invalid_argument &e)
                    {
                        throw runtime_error("Could not convert '" + parser.getArg2() + "' to integer");
                    }
                    catch(out_of_range &e)
                    {
                        throw runtime_error("'" + parser.getArg2() + "' does not fit in an integer");
                    }

                    break;
                }

                default:
                {
                    throw runtime_error("Unrecognized command");
                    break;
                }
            }
        }
    }
    catch(runtime_error &e)
    {
        throw runtime_error(inputFilename + " (" + to_string(parser.getLineNumber()) + "): " + e.what());
    }
}

void Translator::run()
{
    // if the file is a directory, parse and translate all .vm files
    if(is_directory(inputFilename))
    {
        string outputFilename = path(inputFilename).filename().string() + "." + OUTPUT_PREFIX;
        std::ofstream out(outputFilename);
        CodeWriter cw(out);

        bool empty = true;

        cw.writeInit();

        for(auto&& f : directory_iterator(path(inputFilename)))
        {
            if(is_vm(f.path().filename().string()))
            {
                empty = false;
                string input_name = path(f).string();
                translate(input_name, cw);
            }
        }

        if(empty)
        {
            out.close();
            remove(outputFilename);
            throw runtime_error("Error: Directory '" + inputFilename + "' has no .vm files");
        }
    }
    // otherwise, if the file is a regular file, make sure it's a vm file and translate it
    else if(is_regular_file(inputFilename))
    {
        if(is_vm(path(inputFilename).filename().string()))
        {
            std::ofstream out(path(inputFilename).filename().stem().string() + + "." + OUTPUT_PREFIX);
            CodeWriter cw(out);

            cw.writeInit();

            translate(inputFilename, cw);
        }
        else
            throw runtime_error("Error: '" + inputFilename + "' does not end in a .vm extension");
    }
    // otherwise it's not a supported file
    else
        throw runtime_error("Error: '" + inputFilename + "' is not a .vm file or directory");

}

int main(int argc, char *argv[])
{
    // collect arguments passed to the program
    vector<string> arguments;

	for(int i = 0; i < argc; i++)
    {
        arguments.push_back(argv[i]);
    }

    try
    {
        // pass them to the assembler
        Translator program(arguments);
        program.run();
    }
    catch(exception &e)
    {
        cerr << e.what() << endl;
        return 1;
    }
    catch(...)
    {
        cerr << "Unknown error!" << endl;
        return 2;
    }

    return 0;
}


text:
/*
Copyright (c) 2019, NVIDIA CORPORATION. All rights reserved.

NVIDIA CORPORATION and its licensors retain all intellectual property
and proprietary rights in and to this software, related documentation
and any modifications thereto. Any use, reproduction, disclosure or
distribution of this software and related documentation without an express
license agreement from NVIDIA CORPORATION is strictly prohibited.
*/

#include "packages/pnp/gems/epnp/epnp.hpp"

#include <Eigen/SVD>

#include <cmath>
#include <utility>

namespace isaac {
namespace pnp {
namespace epnp {

// High-level function that implements the full EPnP pipeline
pnp::Status ComputeCameraPose(double focal_u, double focal_v, double principal_u,
                              double principal_v, const Matrix3Xd& points3,
                              const Matrix2Xd& points2, epnp::Result* result) {
  if (result == nullptr) {
    return pnp::Status::kErrorBadInputParams;
  }

  // Make sure focal lengths in pixels are non-zero.
  // There is no restriction on the principal point location.
  if (focal_u == 0 || focal_v == 0) {
    return pnp::Status::kErrorBadInputParams;
  }

  // At least 6 2D-3D point correspondences are required currently.
  // TODO: extend to 5 and 4 (requires complicated relinearization)
  if (points2.cols() != points3.cols() || points3.cols() < 6) {
    return pnp::Status::kErrorBadNumInputPoints;
  }

  // Compute control points in world coordinates. They represent a chosen 3D basis.
  constexpr double tol = 1e-3;
  result->ctl_points_world = ChooseBasis(points3, tol);

  // The dimensionality of the point cloud is the number of control points minus one because
  // the first control point is always the centroid.
  // Except when zero input points passed to ChooseBasis().
  result->input_dims = result->ctl_points_world.cols();
  if (result->input_dims) {
    result->input_dims--;
  }

  // Only planar and 3D input point cloud shapes are useful for camera pose estimation.
  if (result->input_dims < 2 || result->input_dims > 3) {
    return pnp::Status::kErrorDegenerateGeom;
  }

  // Compute barycentric coordinates of all input 3D points with respect to chosen basis.
  result->bary_coeffs = ComputeBaryCoords(points3, result->ctl_points_world);
  if (result->bary_coeffs.cols() == 0) {
    return pnp::Status::kErrorDegenerateGeom;
  }

  // Calculate space of the control point camera coordinates such that all
  // projection constraints are approximately satisified.
  if (!SolveProjConstraints(focal_u, focal_v, principal_u, principal_v, result->bary_coeffs,
                            points2, &result->proj_coeffs, &result->solution_basis,
                            &result->singular_values)) {
    return pnp::Status::kErrorDegenerateGeom;
  }

  // Precompute distances between control points in the world frame.
  VectorXd ctl_distances = ComputeDistances(result->ctl_points_world);

  // Construct calibration matrix needed for reprojections at residual calculation.
  Matrix3d calib_matrix;
  calib_matrix << focal_u, 0, principal_u, 0, focal_v, principal_v, 0, 0, 1;

  // Given the solution basis, solve for control point camera coordinates.
  // Solve for various possible dimensions of the solution space and retain
  // the resulting pose that has the least reprojection error.
  const int max_solution_dims = 3;  // Only 1,2,3 are supported
  double min_rss_error = Inf;
  result->rss_repr_errors = VectorXd::Constant(max_solution_dims, Inf);
  for (int k = 0; k < max_solution_dims; k++) {
    const int sol_dims = k + 1;

    // Compute weights that describe the concrete solution in the solution space.
    Vector4d weights = SolveControlPoints(result->solution_basis, sol_dims, ctl_distances);
    if (weights == Vector4d::Zero()) {
      continue;
    }

    // Calculate control points in camera frame given the weights and the base vectors.
    Matrix3Xd ctl_points_cam = ReshapeToMatrix3xN(result->solution_basis * weights);

    // Compute rigid-body transform between control points in world and in camera frame.
    // The resulting transformation is the sought camera pose.
    Matrix3d rotation;
    Vector3d translation;
    if (!pnp::RegisterPointsEucl(result->ctl_points_world, ctl_points_cam, &rotation,
                                 &translation)) {
      continue;
    }

    // Compute squared reprojection errors of the input points, given the computed camera pose.
    VectorXd repr_errors = pnp::ColwiseSquaredNorms(
        points2 - pnp::EuclideanFromHomogeneous(calib_matrix *
                                                ((rotation * points3).colwise() + translation)));

    // Calculate reprojection Residual Sum of Squares (RSS).
    double rss_error = repr_errors.sum();
    result->rss_repr_errors(k) = rss_error;

    // Keep the best solution in terms of reprojection RSS.
    if (rss_error < min_rss_error) {
      min_rss_error = rss_error;
      result->ctl_points_cam = ctl_points_cam;
      result->rotation = rotation;
      result->translation = translation;
      result->solution_dims = sol_dims;
    }
  }

  // Solution could not be found with any dimensionality of the solution space.
  if (min_rss_error == Inf) {
    return pnp::Status::kErrorDegenerateGeom;
  }

  return pnp::Status::kSuccess;
}

// Choose a suitable 3D basis to represent the input 3D points in (via barycentric coordinates).
// Principal Component Analysis of the 3D point cloud.
Matrix3Xd ChooseBasis(const Matrix3Xd& points3, double tol) {
  // The output basis will be represented by a small set of 3D control points (up to 4).
  Matrix3Xd ctl_points;
  ctl_points.resize(3, 0);

  // At least 1 input point is required for further calculations to make sense.
  const int num_points = points3.cols();
  if (num_points < 1) {
    return ctl_points;
  }

  // Calculate centroid of the point cloud - requires at least 1 point.
  Vector3d centroid = points3.rowwise().sum() / num_points;

  // Center the point cloud - requires at least 1 point
  Matrix3Xd centered_points3 = points3;
  for (int i = 0; i < num_points; i++) {
    centered_points3.col(i) -= centroid;
  }

  // Compute principal components.
  // Singular values are sorted in decreasing order.
  Eigen::JacobiSVD<Eigen::MatrixXd> svd;
  svd.compute(centered_points3, Eigen::ComputeFullU | Eigen::ComputeThinV);
  const VectorXd& singular_values = svd.singularValues();
  const MatrixXd& principal_directions = svd.matrixU();

  // Determine the number of non-vanishing principal components given a tolerance.
  int dims = 0;
  for (int i = 0; i < singular_values.rows(); i++) {
    if (singular_values(i) > tol) {
      dims++;
    }
  }

  // The output matrix is always of size 3-by-(dims+1) with the first column being the centroid
  // of the input point cloud.
  ctl_points = Matrix3Xd::Zero(3, dims + 1);
  ctl_points.col(0) = centroid;

  // Scale the other control points with respect to the centroid along the principal directions
  // such that their distance from the centroid equals the square root of the point variance.
  // This is only done for non-vanishing principal directions.
  const double factor = 1.0 / std::sqrt(static_cast<double>(num_points));
  for (int i = 0; i < dims; i++) {
    ctl_points.col(i + 1) = principal_directions.col(i) * factor * singular_values(i) + centroid;
  }

  return ctl_points;
}

// Compute barycentric coords of points with respect to a basis defined by 3 or 4 control points.
// Solves the linear problem A*X = B for X where:
// B is the 4xN matrix of homogeneous 3D points (3xN matrix points3 extended by a row of 1's)
// A is the 4x3 or 4x4 homogenenous matrix of control point coordinates,
// X is the 3xN or 4xN matrix of barycentric coordinates.
// Since the last columns of A and B are all 1's, the column sums of X are 1.
MatrixXd ComputeBaryCoords(const Matrix3Xd& points3, const Matrix3Xd& ctl_points) {
  MatrixXd bary_coords;
  bary_coords.resize(ctl_points.cols(), 0);
  if (points3.cols() && (ctl_points.cols() == 3 || ctl_points.cols() == 4)) {
    Matrix4Xd basis = pnp::HomogeneousFromEuclidean(ctl_points);

    // Solve the problem A*X = B
    bary_coords = basis.householderQr().solve(pnp::HomogeneousFromEuclidean(points3));
  }
  return bary_coords;
}

// Solve the projection equations via Singular-Value Decomposition and
// find the solution space of the control point camera coordinates.
bool SolveProjConstraints(double focal_u, double focal_v, double principal_u, double principal_v,
                          const MatrixXd& bary_coords, const Matrix2Xd points2,
                          MatrixXd* proj_coeffs, MatrixXd* sol_basis, VectorXd* singular_values) {
  if (proj_coeffs == nullptr || sol_basis == nullptr || singular_values == nullptr) {
    return false;
  }

  proj_coeffs->resize(0, 0);

  // At least 6 input points are required, each generating 2 equations.
  int num_points = points2.cols();
  if (num_points < 6) {
    return false;
  }
  if (bary_coords.cols() != num_points) {
    return false;
  }

  // Focal lengths should be non-zero.
  if (focal_u == 0 || focal_v == 0) {
    return false;
  }

  // Shorthands to make it easier to visually verify projection coefficients below.
  const double fu = focal_u;
  const double fv = focal_v;
  const double uc = principal_u;
  const double vc = principal_v;

  // Non-planar case with 4 control points (12 unknowns)
  if (bary_coords.rows() == 4) {
    // Construct the 2Nx12 coefficient matrix
    proj_coeffs->resize(2 * num_points, 12);
    for (int i = 0; i < num_points; i++) {
      // Shorthands to make it easier to check coefficient matrix
      const double u = points2(0, i);
      const double v = points2(1, i);
      const double b1 = bary_coords(0, i);
      const double b2 = bary_coords(1, i);
      const double b3 = bary_coords(2, i);
      const double b4 = bary_coords(3, i);
      proj_coeffs->row(2 * i + 0) << b1 * fu, 0, b1 * (uc - u), b2 * fu, 0, b2 * (uc - u), b3 * fu,
          0, b3 * (uc - u), b4 * fu, 0, b4 * (uc - u);
      proj_coeffs->row(2 * i + 1) << 0, b1 * fv, b1 * (vc - v), 0, b2 * fv, b2 * (vc - v), 0,
          b3 * fv, b3 * (vc - v), 0, b4 * fv, b4 * (vc - v);
    }

    // Planar case with 3 control points (9 unknowns)
  } else if (bary_coords.rows() == 3) {
    // Construct the 2Nx9 coefficient matrix
    proj_coeffs->resize(2 * num_points, 9);
    for (int i = 0; i < num_points; i++) {
      // Shorthands to make it easier to check coefficient matrix
      const double u = points2(0, i);
      const double v = points2(1, i);
      const double b1 = bary_coords(0, i);
      const double b2 = bary_coords(1, i);
      const double b3 = bary_coords(2, i);
      proj_coeffs->row(2 * i + 0) << b1 * fu, 0, b1 * (uc - u), b2 * fu, 0, b2 * (uc - u), b3 * fu,
          0, b3 * (uc - u);
      proj_coeffs->row(2 * i + 1) << 0, b1 * fv, b1 * (vc - v), 0, b2 * fv, b2 * (vc - v), 0,
          b3 * fv, b3 * (vc - v);
    }

    // Cases that are not solved with the current implementation.
  } else {
    return false;
  }

  // Solve the homogeneous linear problem by SVD
  // Matrix V is 9x9 (planar case) or 12x12 (non-planar case)
  // There are 9 (planar case) or 12 singular values  (non-planar case)
  Eigen::JacobiSVD<Eigen::MatrixXd> svd;
  svd.compute(*proj_coeffs, Eigen::ComputeThinU | Eigen::ComputeFullV);
  *singular_values = svd.singularValues();
  const MatrixXd& singular_vectors = svd.matrixV();

  // Output the last 4 columns of V in reverse order
  sol_basis->resize(singular_vectors.rows(), 4);
  for (int i = 0; i < 4; i++) {
    sol_basis->col(i) = singular_vectors.col(singular_vectors.cols() - i - 1);
  }

  // Put singular values in increasing order
  for (int i = 0; i < singular_values->size() / 2; i++) {
    std::swap((*singular_values)(i), (*singular_values)(singular_values->size() - 1 - i));
  }

  return true;
}

// Convert a solution vector of control point coordinates of length 3*C to a 3xC matrix for C=3,4
// such that the input vector contains the elements of the output matrix column-wise.
Matrix3Xd ReshapeToMatrix3xN(const VectorXd& vec) {
  Matrix3Xd matrix;
  if (vec.size() % 3) {
    return matrix;  // return a 3x0 matrix
  }
  matrix.resize(3, vec.size() / 3);
  for (int i = 0; i < matrix.cols(); i++) {
    matrix.col(i) = vec.block(3 * i, 0, 3, 1);
  }

  return matrix;
}

// Project points directly given in camera coordinates to the image
Matrix2Xd ProjectPoints(const Matrix3d& calib_matrix, const Matrix3Xd& points3) {
  Matrix2Xd points2;
  points2.resize(2, points3.cols());
  for (int i = 0; i < points3.cols(); i++) {
    Vector3d proj = calib_matrix * points3.col(i);
    if (proj(2) == 0) {
      // Avoids NaNs due to 0/0.
      points2(0, i) = Inf;
      points2(1, i) = Inf;
    } else {
      // Division by homogeneous coordinate.
      points2(0, i) = proj(0) / proj(2);
      points2(1, i) = proj(1) / proj(2);
    }
  }
  return points2;
}

// Compute distances between all possible pairs (i,j) of input points (in matrix columns).
VectorXd ComputeDistances(const MatrixXd& points) {
  int num_points = points.cols();
  if (num_points < 2) {
    return VectorXd(0);
  }

  VectorXd distances(num_points * (num_points - 1) / 2);
  int k = 0;
  for (int i = 0; i < points.cols(); i++) {
    for (int j = i + 1; j < points.cols(); j++) {
      distances[k++] = (points.col(i) - points.col(j)).norm();
    }
  }

  return distances;
}

// Compute the weights that correspond to a concrete solution for the control point camera
// coordinates in the basis given as input.
Vector4d SolveControlPoints(const MatrixXd& sol_basis, int sol_dims, const VectorXd& distances) {
  Vector4d weights = Vector4d::Zero();
  if (sol_basis.cols() != 4) {
    return weights;
  }
  if (sol_basis.rows() != 12 && sol_basis.rows() != 9) {
    return weights;
  }

  // non-planar case
  if (sol_basis.rows() == 12 && distances.size() != 6) {
    return weights;
  }

  // planar case
  if (sol_basis.rows() == 9 && distances.size() != 3) {
    return weights;
  }

  // Planar and non-planar case for 1-dimensional solution space
  if (sol_dims == 1) {
    // Only use singular vector corresponding to the least singular value.
    // Solution space is 1-dimensional -> only need to find a single global scale.

    // Calculate the unscaled camera coordinates of the control points.
    Matrix3Xd v = ReshapeToMatrix3xN(sol_basis.col(0));

    // Calculate distances between the unscaled control points.
    VectorXd distances_cam = ComputeDistances(v);

    // Calculate the unknown scale from the distances between control points in the world frame
    // and the corresponding distances in the camera frame.
    double scale = distances.dot(distances_cam) / distances_cam.dot(distances_cam);

    // Find sign of the scale such that a majority of control points are in front of the camera.
    int num_points_behind = (v.row(2).array() < 0).count();
    if (num_points_behind > 1) {
      scale = -scale;
    }

    weights(0) = scale;

    // Planar and non-planar for 2-dimensional solution space
  } else if (sol_dims == 2) {
    // Reshape the two basis vectors of the solution space into two 3xN matrices.
    // The sought matrix ctl_points_cam is a linear combination of these matrices.
    // ctl_points_cam = w1*v1 + w2*v2, where the weights w1 and w2 are still unknown.
    Matrix3Xd v1 = ReshapeToMatrix3xN(sol_basis.col(0));
    Matrix3Xd v2 = ReshapeToMatrix3xN(sol_basis.col(1));

    // Looking for the linear combination preserving the distances between control points.
    // This is a non-linear problem in the linear combination weights (w1,w2).
    // Strategy:
    //   (Step1) Solve for (w11,w22,w12) = (w1*w1, w2*w2, w1*w2) (linearized problem)
    //   (Step2) Extract (w1,w2) from (w11,w22,w12).
    // Both steps are over-determined problems.

    // Construct coefficient matrix of linearized problem (see Eq.13 in the paper).
    // This matrix is 6x4 in the non-planar case.
    const int num_points = v1.cols();
    const int num_equations = num_points * (num_points - 1) / 2;
    int k = 0;
    MatrixXd coeff_matrix(num_equations, 3);
    for (int i = 0; i < num_points; i++) {
      for (int j = i + 1; j < num_points; j++) {
        Vector3d d1 = v1.col(i) - v1.col(j);
        Vector3d d2 = v2.col(i) - v2.col(j);
        coeff_matrix.row(k++) << d1.dot(d1), d2.dot(d2), 2 * d1.dot(d2);
      }
    }

    // Step1: Compute least squares solution to the linearized problem.
    VectorXd squared_distances = distances.cwiseProduct(distances);
    VectorXd lin_sol =
        coeff_matrix.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(squared_distances);

    // Solution of the linearized problem.
    const double w11 = lin_sol(0);
    const double w22 = lin_sol(1);
    const double w12 = lin_sol(2);

    // Step2: Calculate weights w1 and w2 (denoted beta_i in the EPnP paper)
    // from the non-linear constraints (w11,w22,w12) = (w1*w1, w2*w2, w1*w2).
    // There are always two solutions, (w1,w2) and (-w1,-w2) -> sign is solved further below.
    // Over-constrained quadratic problem - least-squares solution has no simple closed form.
    // Empirical tests show that solution order largely affects pose accuracy.
    // Strategy: we exactly match the larger magnitude wij and hope for less overall alg. error.
    double w1, w2;
    w1 = w2 = 0;
    if (w11 >= w22 && w11 > 0) {
      w1 = std::sqrt(w11);
      if (w22 > std::abs(w12)) {
        w2 = std::sqrt(w22);
      } else {
        w2 = w12 / w1;
      }
    } else if (w11 <= w22 && w22 > 0) {
      w2 = std::sqrt(w22);
      if (w11 > std::abs(w12)) {
        w1 = std::sqrt(w11);
      } else {
        w1 = w12 / w2;
      }
    }

    // There are always two solutions: (w1,w2) and (-w1,-w2) but
    // one corresponds to control points behind the camera.
    // Find sign such that a majority all control points are in front of the camera.
    VectorXd ctl_point_depths = w1 * v1.row(2) + w2 * v2.row(2);
    int num_points_behind = (ctl_point_depths.array() < 0).count();
    if (num_points_behind > int{num_points / 2}) {
      weights << -w1, -w2, 0, 0;
    } else {
      weights << w1, w2, 0, 0;
    }

    // Non-planar case for 2-dimensional solution space
  } else if (sol_dims == 3 && sol_basis.rows() == 12) {
    // Reshape the 3 basis vectors of the solution space into three 3xN matrices.
    // The sought matrix ctl_points_cam is a linear combination of these matrices.
    // ctl_points_cam = w1*v1 + w2*v2 + w3*v3, where the weights (w1,w2,w3) are still unknown.
    Matrix3Xd v1 = ReshapeToMatrix3xN(sol_basis.col(0));
    Matrix3Xd v2 = ReshapeToMatrix3xN(sol_basis.col(1));
    Matrix3Xd v3 = ReshapeToMatrix3xN(sol_basis.col(2));

    // Looking for the linear combination preserving the distances between control points.
    // This is a non-linear problem in the linear combination weights (w1,w2,w3).
    // Strategy:
    //   (Step1) Solve for (w11,w22,w33,w12,w13,w23) = (w1*w1, w2*w2, w3*w3, w1*w2, w1*w3, w2*w3)
    //           (linearized problem).
    //   (Step2) Extract (w1,w2,w3) from (w11,w22,w33,w12,w13,w23).

    // Step1: Construct coefficient matrix of linearized problem (Case N=3 in the paper).
    // This matrix is 6x6 in the non-planar case.
    const int num_points = v1.cols();
    const int num_equations = num_points * (num_points - 1) / 2;
    int k = 0;
    MatrixXd coeff_matrix(num_equations, 6);
    for (int i = 0; i < num_points; i++) {
      for (int j = i + 1; j < num_points; j++) {
        Vector3d d1 = v1.col(i) - v1.col(j);
        Vector3d d2 = v2.col(i) - v2.col(j);
        Vector3d d3 = v3.col(i) - v3.col(j);
        coeff_matrix.row(k++) << d1.dot(d1), d2.dot(d2), d3.dot(d3), 2 * d1.dot(d2), 2 * d1.dot(d3),
            2 * d2.dot(d3);
      }
    }

    // Step1: Compute least squares solution to the linearized problem.
    VectorXd squared_distances = distances.cwiseProduct(distances);
    VectorXd lin_sol =
        coeff_matrix.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(squared_distances);

    // Solution of the linearized problem.
    const double w11 = lin_sol(0);
    const double w22 = lin_sol(1);
    const double w33 = lin_sol(2);
    const double w12 = lin_sol(3);
    const double w13 = lin_sol(4);
    const double w23 = lin_sol(5);

    // Calculate weights w1, w2, w3 from the non-linear constraints
    // (w11, w22, w33, w12, w13, w23) = (w1*w1, w2*w2, w3*w3, w1*w2, w1*w3, w2*w3).
    // There are always two solutions (w1,w2,w3) and (-w1,-w2,-w3) -> sign solved further below.
    // Over-constrained quadratic problem - leats-squares solution has no simple closed form.
    // Empirical tests show that solution order largely affects pose accuracy.
    // Strategy: we exactly match the larger magnitude wij and hope for less overall alg. error.
    double w1 = 0, w2 = 0, w3 = 0;
    if (w11 > w22 && w11 > w33 && w11 > 0) {
      w1 = std::sqrt(w11);
      if (w22 > std::abs(w12)) {
        w2 = std::sqrt(w22);
      } else {
        w2 = w12 / w1;
      }
      if (w33 > std::abs(w13)) {
        w3 = std::sqrt(w33);
      } else {
        w3 = w13 / w1;
      }
    } else if (w22 > w11 && w22 > w33 && w22 > 0) {
      w2 = std::sqrt(w22);
      if (w11 > std::abs(w12)) {
        w1 = std::sqrt(w11);
      } else {
        w1 = w12 / w2;
      }
      if (w33 > std::abs(w23)) {
        w3 = std::sqrt(w33);
      } else {
        w3 = w23 / w2;
      }
    } else if (w33 > w11 && w33 > w22 && w33 > 0) {
      w3 = std::sqrt(w33);
      if (w11 > std::abs(w13)) {
        w1 = std::sqrt(w11);
      } else {
        w1 = w13 / w3;
      }
      if (w22 > std::abs(w23)) {
        w2 = std::sqrt(w22);
      } else {
        w2 = w23 / w3;
      }
    }

    // There are always two solutions: (w1,w2,w3) and (-w1,-w2,-w3) but
    // one corresponds to control points behind the camera.
    // Find sign such that a majority all control points are in front of the camera.
    VectorXd ctl_point_depths = w1 * v1.row(2) + w2 * v2.row(2) + w3 * v3.row(2);
    int num_points_behind = (ctl_point_depths.array() < 0).count();
    if (num_points_behind > int{num_points / 2}) {
      weights << -w1, -w2, -w3, 0;
    } else {
      weights << w1, w2, w3, 0;
    }
  }

  // TODO: planar N=3 case using relinearization
  // TODO: non-planar N=4 (affine camera) case using relinearization
  // planar N=4 case is underconstrained = ambiguous (3 equations but 4 unknowns)

  return weights;
}

}  // namespace epnp
}  // namespace pnp
}  // namespace isaac


text:
#include <assert.h>
#include <iostream>
#include <vector>
#include <boost/shared_ptr.hpp>
#include "stardate.h"
#include "lin_time.h"
#include "vsfilesystem.h"
#include "cmd/script/mission.h"

using std::cerr;
using std::cout;
using std::endl;
using std::vector;

class Faction;
extern vector< boost::shared_ptr<Faction> >factions;

StarDate::StarDate()
{
    initial_time = mission->getGametime();
    initial_star_time = NULL;
}

void StarDate::Init( double time )
{
    if (initial_star_time != NULL)
        delete[] initial_star_time;
    initial_time = mission->getGametime();
    initial_star_time = new double[factions.size()];
    for (unsigned int i = 0; i < factions.size(); i++)
        initial_star_time[i] = time;
}

//Get the current StarDate time in seconds
double StarDate::GetCurrentStarTime( int faction )
{
    //Get the number of seconds elapsed since the server start
    double time_since_server_started = mission->getGametime()-initial_time;
    //Add them to the current date
    if (initial_star_time == NULL)
        return time_since_server_started;
    else
        return initial_star_time[faction]+time_since_server_started;
}

//Needed to calculate relative message and mission times
//into stardate
double StarDate::GetElapsedStarTime(int faction)
{
    if (initial_star_time == NULL)
        return initial_time;
    else
        return initial_star_time[faction] - initial_time;
}

/*
 *********************************************************************************
 **** Trek Stardate stuff                                                      ***
 *********************************************************************************
 */

void StarDate::InitTrek( string date )
{
    if (initial_star_time != NULL)
        //we must be reinitializing;
        delete[] initial_star_time;
    initial_star_time = 0;
    initial_time = mission->getGametime();
    initial_star_time = new double[factions.size()];
    double init_time = this->ConvertTrekDate( date );
    VSFileSystem::vs_dprintf( 3, "Initializing stardate from a Trek date for %d factions", factions.size() );
    for (unsigned int i = 0; i < factions.size(); i++)
        initial_star_time[i] = init_time;
}

/// The stardate format is ddd.hhmm:sec with
/// ddd, hh, mm, sec = star day, hour, minute, seconds
/// The second equals our terran second. The minute however has 480 seconds (HOURS_DIV=8)
/// The hour has 60 minutes and the day has 100 hours

//Convert a StarDate time into a Stardate string
string StarDate::ConvertFullTrekDate( double date )
{
    unsigned int days, hours, minutes, seconds;
    char cdate[32];

    // Get the number of days dividing seconds number by the number of seconds in a day: 100*60*60*8 = 2880000
    days    = (unsigned int) date/2880000;
    // Modulo gives us the number of stardate seconds elapsed in the current day
    date    = (unsigned int) date%2880000;
    // Get the hours elapsed in the day by dividing by number of seconds in a stardate hour: 60*60*8 = 28800
    hours   = (unsigned int) date/28800;
    // Modulo gives us the number of seconds elapsed in that hour
    date    = (unsigned int) date%28800;
    //Get the number of minutes elapsed in that hour by dividing by the number of seconds in a minute: 60*8 = 480
    minutes = (unsigned int) date/480;
    //The remaining seconds in the date
    //The seconds are counted from 0 to 480 before the minute effectively is incremented
    seconds = (unsigned int) date%480;

    //The hour/minute part is displayed like HHMM divided by HOURS_DIV which is 8 for now
    unsigned int hhmm  = (hours*100+minutes);

    sprintf( cdate, "%d.%.4d:%.3d", days, hhmm, seconds );
    return string( cdate );
}

string StarDate::ConvertTrekDate( double date )
{
    unsigned int days, hours, minutes;
    char cdate[32];

    days    = (unsigned int) date/2880000;
    date    = (unsigned int) date%2880000;
    hours   = (unsigned int) date/28800;
    date    = (unsigned int) date%28800;
    minutes = (unsigned int) date/480;

    unsigned int hhmm = (hours*100+minutes);

    sprintf( cdate, "%d.%.4d", days, hhmm );
    return string( cdate );
}

//Convert a StarDate into a number of seconds
double StarDate::ConvertTrekDate( string date )
{
    unsigned int days, hours, minutes, tmphrs, seconds, nb, pos;
    double res;
    //Replace the dot with 'a' so sscanf won't take it for a decimal symbol
    pos = date.find( "." );
    date.replace( pos, 1, "a" );
    if ( ( nb = sscanf( date.c_str(), "%da%4d:%3d", &days, &tmphrs, &seconds ) ) != 3 )
        VSFileSystem::vs_dprintf( 3, "!!! ERROR reading date\n");

    //Extract number of hours
    hours   = tmphrs/100;
    //Extract number of minutes
    minutes = tmphrs%100;

    res     = days*2880000+hours*28800+minutes*480+seconds;
    std::string formatted = ConvertFullTrekDate(res);
    VSFileSystem::vs_dprintf( 3, "Converted date to %ld, which stardate is %s\n", long(res), formatted.c_str() );
    return res;
}

//Get the current StarDate in a string
string StarDate::GetFullTrekDate( int faction )
{
    return ConvertFullTrekDate( this->GetCurrentStarTime( faction ) );
}

//Get the current StarDate in a string - short format
string StarDate::GetTrekDate( int faction )
{
    return ConvertTrekDate( this->GetCurrentStarTime( faction ) );
}

//Convert the string xxxx.y date format into a float representing the same data xxxx.y
float StarDate::GetFloatFromTrekDate( int faction )
{
    float  float_date;
    string cur_date = this->GetFullTrekDate( faction );
    sscanf( cur_date.c_str(), "%f", &float_date );

    return float_date;
}

/*
 *********************************************************************************
 **** DAN.A Stardate stuff                                                     ***
 *********************************************************************************
 */

void InitSDate( string date )
{}

string GetSDate( int faction = 0 )
{
    return string( "" );
}

string GetFullSDate( int faction = 0 )
{
    return string( "" );
}

/*
 *********************************************************************************
 **** String stardate formats conversions                                      ***
 *********************************************************************************
 */

string SDateFromTrekDate( string trekdate )
{
    return string( "" );
}

string TrekDateFromSDate( string sdate )
{
    return string( "" );
}