btl_vector_algorithms.h

利用这个模板可以分析基因表达数据

//..........................................................................................
    	    /**#: [Description="Returns the sum of the elements of a triple."]*/
 
 	template<class InputIterator, class T>
    	T sum(InputIterator firsta, T init)
	{
	    return init += *firsta + *(firsta+1) + *(firsta+2);
	}

//..........................................................................................
    	    /**#: [Description="Returns the sum of the vector elements, 
                   but accumulates the elements in order 
                   of ascending value in order to minimise rounding errors."]*/
 
	template<class InputIterator, class T>
    	T sum_precise(InputIterator first, InputIterator last, T init)
	{
            typename vector<T>::const_iterator i;
	    vector<T> v(first,last);
	    sort(v.begin(), v.end());
	    for (i=v.begin(); i!=v.end(); i++)
	    	init += *i;
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Returns the sum of the squares of the vector elements."]*/
  
	template<class InputIterator, class T>
    	T sum_of_squares(InputIterator first, InputIterator last, T init)
	{
            for ( ; first != last; first++)
	        init += *first * *first;
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Returns the sum of the squares of elements of a triple."]*/
 
 	template<class InputIterator, class T>
    	T sum_of_squares(InputIterator firsta, T init)
	{
	    return init += (*firsta * *firsta) + (*(firsta+1) * *(firsta+1)) 
			                        + (*(firsta+2) * *(firsta+2));
	}

//..........................................................................................
    	    /**#: [Description="Returns the sum of the squares of the vector elements, 
                   but accumulates the elements in order 
                   of ascending value in order to minimise rounding errors."]*/
 
	template<class InputIterator, class T>
    	T sum_of_squares_precise(InputIterator first, InputIterator last, T init)
	{
            typename vector<T>::const_iterator i;
	    vector<T> v(first,last);
	    sort(v.begin(), v.end());
	    for (i=v.begin(); i!=v.end(); i++)
	    	init += *i * *i;
	    return init;
	}
//..........................................................................................
    	    /**#: [Description="Returns the length of a vector."]*/
  
	template<class InputIterator1, class T>
    	T magnitude(InputIterator1 firsta, InputIterator1 lasta, T init )
	{
	    return sqrt(sum_of_squares(firsta,lasta,init));
	}

//..........................................................................................
    	    /**#: [Description="Returns the precise length of a vector."]*/
  
	template<class InputIterator1, class T>
    	T magnitude_precise(InputIterator1 firsta, InputIterator1 lasta, T init )
	{
	    return sqrt(sum_of_squares_precise(firsta,lasta,init));
	}

//..........................................................................................
    	    /**#: [Description="Returns the length of a triple."]*/
 
 	template<class InputIterator1, class T>
    	T magnitude(InputIterator1 firsta, T init)
	{
	    return sqrt(sum_of_squares(firsta,init));
	}

//} // end of namespace

//.............................................................................

    	    /**#: [Description="Rotates each triple in a container using a given 
	                        rotation matrix and origin"]*/

template <class ForwardIterator, class ConstForwardIterator1, class ConstForwardIterator2>
void rotate(ForwardIterator begin, ForwardIterator end, 
	    ConstForwardIterator1 rotmatrix, ConstForwardIterator2 origin)
{
	for(; begin!=end; begin+=3)
	{
	    *begin     -= *origin;
 	    *(begin+1) -= *(origin+1);
	    *(begin+2) -= *(origin+2);
	    BTL_REAL temp_x = *begin;
           BTL_REAL temp_y = *(begin+1);

           *begin     = (*begin * *rotmatrix)     + (*(begin+1) * *(rotmatrix+1))
                      + (*(begin+2) * *(rotmatrix+2));

           *(begin+1) = (temp_x * *(rotmatrix+3)) + (temp_y * *(rotmatrix+4))
                      + (*(begin+2) * *(rotmatrix+5));

           *(begin+2) = (temp_x * *(rotmatrix+6)) + (temp_y * *(rotmatrix+7))
                      + (*(begin+2) * *(rotmatrix+8));
	    *begin     += *origin;
 	    *(begin+1) += *(origin+1);
	    *(begin+2) += *(origin+2);
	}
        
}                    				       
//.............................................................................
    	    /**#: [Description="Rotates each triple in a container about a given 
	                        axis and origin by a given number of radians."] */

template <class ForwardIterator, class ConstForwardIterator1, class ConstForwardIterator2, class T>
void rotate(ForwardIterator first, ForwardIterator last, 
	    ConstForwardIterator1 axis, ConstForwardIterator2 origin, T angle)
	{
	// Constructor for 3x3 rotation matrix for a rotation about a given axis by a
	// given number of radians. Create the following matrix:

	//  [cosA+a1x.a1x(1-cosA)      -a3xsinA+a1a2(1-cosA)	a2xsinA+a1xa3x(1-cosA) ] 
	//  [a3xsinA+a1xa2(1-cosA)    cosA+a2xa2x(1-cosA)   	-a1xsinA+a2xa3x(1-cosA)] 
	//  [-a2xsinA+a3xa1x(1-cosA)   a1xsinA+(1-cosA)a2xa3x 	cosA+(1-cosA)a3xa3x   ]

	// where (a1x,a2x,a3x) is the unit vector in the direction of the axis and A is the angle
	// of rotation in radians.

            BTL_REAL recipmagofaxis = 1.0/sqrt(((*axis * *axis)+(*(axis+1) * *(axis+1))+(*(axis+2) * *(axis+2)))) ;
	    BTL_REAL a1x     = *axis * recipmagofaxis;
	    BTL_REAL a2x     = *(axis+1) * recipmagofaxis;
	    BTL_REAL a3x     = *(axis+2) * recipmagofaxis;
	    BTL_REAL cosA    = cos(angle);
	    BTL_REAL sinA    = sin(angle);
	    BTL_REAL a1mcosA = a1x * (1.0-cosA);    	
	    BTL_REAL a2mcosA = a2x * (1.0-cosA);    	
	    BTL_REAL a3mcosA = a3x * (1.0-cosA);    	
	    BTL_REAL a1sinA  = a1x * sinA;    	
	    BTL_REAL a2sinA  = a2x * sinA;    	
	    BTL_REAL a3sinA  = a3x * sinA;

            vector<BTL_REAL> rotation;
	    rotation.push_back(( cosA    + a1x * a1mcosA));
            rotation.push_back((-a3sinA  + a1x * a2mcosA));
            rotation.push_back(( a2sinA  + a1x * a3mcosA));
            rotation.push_back(( a3sinA  + a1x * a2mcosA));
            rotation.push_back(( cosA    + a2x * a2mcosA));
            rotation.push_back((-a1sinA  + a2x * a3mcosA));
            rotation.push_back((-a2sinA  + a1x * a3mcosA));
            rotation.push_back(( a1sinA  + a2x * a3mcosA));
            rotation.push_back(( cosA    + a3x * a3mcosA));

            rotate(first,last,rotation.begin(),origin);
	}

//.............................................................................
    	    /**#: [Description="Translates a each triple in a container 
	                        by a given vector"]*/

template <class ForwardIterator, class ConstForwardIterator>
void translate(ForwardIterator begin, ForwardIterator end, ConstForwardIterator trans)
{
	for(; begin!=end; begin+=3)
	{
	    *begin     += *trans;
 	    *(begin+1) += *(trans+1);
	    *(begin+2) += *(trans+2);
	}	
}


_BTL_END_NAMESPACE

//    const double DEGRAD = M_PI/180.0; 	// converts degrees to radians
//    double radAngle = angle * DEGRAD;

#endif