btl_vector_algorithms.h

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

//
// btl_vector_algorithms.h
//
// This file contains the generic numeric functions for manipulation of vectors 
//
// These classes are part of the Bioinformatics Template Library (BTL).
//
// Copyright (C) 1997, 1998 Birkbeck College, Malet Street, London, U.K.
// Copyright (C) 2004, 2005 University College, Gower Street, London, U.K.
//
// This library is free software; you can redistribute it and/or modify it
// under the terms of the GNU Library General Public License as published 
// by the Free Software Foundation; either version 2 of the License, or 
// (at your option) any later version.  This library is distributed in the
// hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
// implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE.  See the GNU Library General Public License for more details.
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
///////////////////////////////////////////////////////////////////////////

#if !defined (BTL_VECTORALGORITHMS_H)
#define BTL_VECTORALGORITHMS_H 1

/**#: [Description ="A collection of generic numerical algorithms for
      vector algebra with an emphasis on the manipulation of vectors in 3D."]
    [Summary = "generic vector algorithms"] 
    [Authors = "D.S.Moss, W.R.Pitt, I.Tickle, M.A.Williams"]
    [Files = "<A HREF=./btl/btl_vector_algorithms.h>btl_vector_algorithms.h</A>"]
    [Dependencies="<A HREF=#numeric_vector>btl_numeric_vector.h</A>,
                   <A HREF=#matrix>btl_matrix.h</A>"]
    [Prerequisites="None<P>"
    <P>
*/

#include <vector>
#include <iterator>

#include "BTL.h"
#include "btl_numeric_vector.h"
#include "btl_matrix.h"

_BTL_BEGIN_NAMESPACE

using namespace std;
    	
//..........................................................................................
// VECTOR ALGEBRA
//..........................................................................................
    	    /**#: [Description="Scalar/dot product of two vectors."]*/
 
	template<class InputIterator1, class InputIterator2, class T >
    	T scalar_product(InputIterator1 firsta, InputIterator1 lasta,
                            InputIterator2 firstb, T init)
	{
            for ( ; firsta != lasta; firsta++, firstb++)
	        init += (*firsta * *firstb);
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Scalar/dot product of two triples ."]*/
 
	template<class InputIterator1, class InputIterator2, class T >
    	T scalar_product(InputIterator1 firsta, InputIterator2 firstb, T init)
	{
	    init += (*firsta * *firstb) + (*(firsta+1) * *(firstb+1)) 
                                           + (*(firsta+2) * *(firstb+2));
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Direct product of two matrices or vectors ."]*/

	template<class InputIterator1, class InputIterator2, class OutputIterator>
    	OutputIterator direct_product(InputIterator1 firsta, InputIterator1 lasta, 
                                      InputIterator2 firstb, OutputIterator result)
	{
	    if(firsta == firstb)
		for (; firsta != lasta; firsta++, result++)
		    *result += *firsta * *firsta;
	    else 
		for (; firsta != lasta; firsta++, firstb++, result++)
		    *result += *firsta * *firstb;
	    return ++result;			 
	}
//..........................................................................................
    	    /**#: [Description="Vector/cross product of two triples."]*/
 
	template<class InputIterator1, class InputIterator2, class OutputIterator>
    	OutputIterator vector_product(InputIterator1 firsta, InputIterator2 firstb, 
                                      OutputIterator result)
	{
	    *result	+= (*(firsta+1) * *(firstb+2)) - (*(firsta+2) * *(firstb+1));
	    *(result+1) += (*(firsta+2) * *firstb)     - (*firsta     * *(firstb+2));
	    *(result+2) += (*firsta     * *(firstb+1)) - (*(firsta+1) * *firstb);
	    return (result+3);
	}

//..........................................................................................
    	    /**#: [Description="Vector product of three triples."]*/
 
	template<class InputIterator1, class InputIterator2, class InputIterator3, 
		 class OutputIterator>
    	OutputIterator triple_vector_product(InputIterator1 firsta, InputIterator2 firstb, 
                                             InputIterator3 firstc, OutputIterator result)
	{
            typename iterator_traits<OutputIterator>::value_type adotc = 0.0;
	    adotc = scalar_product(firsta,firstc,adotc);
            typename iterator_traits<OutputIterator>::value_type adotb = 0.0; 
	    adotb = scalar_product(firsta,firstb,adotb);
	    *result     += (*firstb     * adotc) - (*firstc     * adotb);
	    *(result+1) += (*(firstb+1) * adotc) - (*(firstc+1) * adotb);
	    *(result+2) += (*(firstb+2) * adotc) - (*(firstc+2) * adotb);
	    return (result+3);
	}

//..........................................................................................
    	    /**#: [Description="Scalar product of three triples."]*/
 
	template<class InputIterator1, class InputIterator2, class InputIterator3, class T>
    	T triple_scalar_product(InputIterator1 firsta, InputIterator2 firstb, 
                                  InputIterator3 firstc, T init)
	{
	    init += (*firsta * ((*(firstb+1) * *(firstc+2))-(*(firstb+2)* *(firstc+1))))
		  + (*(firsta+1) * ((*(firstb+2) * *firstc)-(*firstb * *(firstc+2))))
		  + (*(firsta+2) * ((*firstb * *(firstc+1))-(*(firstb+1)* *firstc)));
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Calculate the squared distance between the points 
                                 represented by two vectors."]*/
 
	template<class InputIterator1, class InputIterator2, class T>
    	T separation_squared(InputIterator1 firsta, InputIterator1 lasta,
                            InputIterator2 firstb, T init)
	{
            for ( ; firsta != lasta; firsta++, firstb++)
            {
	        init += (*firsta - *firstb) * (*firsta - *firstb);
	    }	
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Calculate the squared distance between points represented 
                                by two triples."]*/
 
	template<class InputIterator1, class InputIterator2, class T>
    	T separation_squared(InputIterator1 firsta, InputIterator2 firstb, T init)
	{
	    init  += ((*firsta - *firstb) * (*firsta - *firstb)) 
	           + ((*(firsta+1) - *(firstb+1)) * (*(firsta+1) - *(firstb+1)))
	           + ((*(firsta+2) - *(firstb+2)) * (*(firsta+2) - *(firstb+2)));
	    return init;
	}

//..........................................................................................
    	    /**#: [Description="Calculate the distance between the points represented 
                                by two vectors."]*/
 
	template<class InputIterator1, class InputIterator2, class T>
    	T separation(InputIterator1 firsta, InputIterator1 lasta,
                            InputIterator2 firstb, T init)
	{
	    return T(sqrt(separation_squared(firsta,lasta,firstb,init)));
	}

//..........................................................................................
    	    /**#: [Description="Calculate the distance between the points represented 
                                by two triples."]*/
 
	template<class InputIterator1, class InputIterator2, class T>
    	T separation(InputIterator1 firsta, InputIterator2 firstb, T init)
	{
	    return T(sqrt(separation_squared(firsta,firstb,init)));
	}


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