example3.cpp

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

//
// This file contains some example code, which may be used to call 
// the matrix, vector and statistics algorithms.
//
// This code is not necessarily meant to be useful in itself, but is provided as 
// an example of how the class may be used.
//
// Copyright (C) 1999 Software Engineering Group, Crystallography Department,
// Birkbeck College, Malet Street, London WC1E 7HX, U.K.
// ([e-mail] or [e-mail])
// 
// 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
// Handle) 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.
/////////////////////////////////////////////////////////////////////////////////////
//
// Author: Mark Williams 
// 
/////////////////////////////////////////////////////////////////////////////////////
// 
// Brief Description of Code:
// 
// Carry out simple numerical tasks.
//
// For further details see the documentation for the matrix, vector classes 
// and the algorithms header files.
/////////////////////////////////////////////////////////////////////////////////////
// 
#include <iostream>
using namespace std;

#include "btl_numeric_vector.h"
#include "btl_matrix.h"
#include "btl_vector_algorithms.h"
#include "btl_matrix_algorithms.h"
#include "btl_statistics.h"
using namespace btl;

int main()
{

	matrix<double> A(3,4,2.3);
	matrix<double> B(4,3,2.6);
	matrix<double> C(3,3,1.2);
	matrix<double> D(3,4);

	// member function - A,B,C are instances of the matrix class
	C = A * B;

// Vector stuff

	double temp = 0.0;
	temp = scalar_product(A.begin(),A.end(),B.begin(),temp);
	cout << "temp " << temp << endl;

	temp = 0.0;
	temp = scalar_product(A.begin(),B.begin(),temp);
	cout << "temp triple " << temp << endl;

	vector_product(A.begin(),B.begin(),C.begin());
	cout << "vector product " << C(1,1) << endl;

	triple_vector_product(A.begin(),B.begin(),C.begin(),D.begin());
	cout << "triple vector product " << D(1,1) << " " << D(1,2) << 
	" " << D(1,3) << endl;

	direct_product(A.begin(),A.end(),B.begin(),D.begin());
	cout << "direct " << D(1,1) << endl;

	D += temp;
	cout << "direct+temp " << D(1,1) << endl;

	temp = 0.0;
	temp = triple_scalar_product(A.begin(),B.begin(),D.begin(),temp);
	cout << "triple_scalar " << temp << endl;

	temp = 0.0;
	temp = separation_squared(A.begin(),A.end(),B.begin(),temp);
	cout << "separation squared " << temp << endl;

	temp = 0.0;
	temp = separation(A.begin(),A.end(),B.begin(),temp);
	cout << "separation " << temp << endl;

	temp = 0.0;
	temp = sum(A.begin(),A.end(),temp);
	cout << "sum "<< temp << endl;
	
	temp = 0.0;
	temp = sum_precise(A.begin(),A.end(),temp);
	cout << "sum precise "<< temp << endl;

	temp = 0.0;
	temp = sum_of_squares(A.begin(),A.end(),temp);
	cout << "sum squares "<< temp << endl;

	temp = 0.0;
	temp = sum_of_squares_precise(A.begin(),A.end(),temp);
	cout << "sum squares precise "<< temp << endl;

	temp = 0.0;
	temp = sum(A.begin(),temp);
	cout << "sum triple "<< temp << endl;

	temp = 0.0;
	temp = sum_of_squares(A.begin(),temp);
	cout << "sum squares triple "<< temp << endl;

	temp = 0.0;
	temp = magnitude(A.begin(),A.end(),temp);
	cout << "magnitude "<< temp << endl;

	temp = 0.0;
	temp = magnitude_precise(A.begin(),A.end(),temp);
	cout << "magnitude precise "<< temp << endl;

	temp = 0.0;
	temp = magnitude(A.begin(),temp);
	cout << "magnitude triple "<< temp << endl;

// Statistics stuff

	double mean1 = 0.0;
	mean1 = mean(A.begin(),A.end(),mean1);
	cout << "mean "<< mean1 << endl;

	temp = 0.0;
	temp = mean_absolute_deviation(A.begin(),A.end(),temp,mean1);
	cout << "mean abs dev "<< temp << endl;

	temp = 0.0;
	temp = variance(A.begin(),A.end(),temp,mean1);
	cout << "variance "<< temp << endl;

	double temp1=0.0;
	double temp2=0.0;
	double temp3=0.0;
	double temp4=0.0;
	double temp5=0.0;
	normal_statistics(A.begin(),A.end(),temp1,temp2,temp3,temp4,temp5);
	cout << "mean "<< temp1 << endl;
	cout << "mad "<< temp2 << endl;
	cout << "variance "<< temp3 << endl;
	cout << "skew "<< temp4 << endl;
	cout << "kurtosis "<< temp << endl;

// Matrix stuff

	vector<double> E(5);
	cout << "A " << A(1,1) << endl;
	column_means(A.begin(),A.end(),A.rows(),E.begin());
	cout << "A column mean " << E[0] << endl;

	copy_column(A.begin(),A.end(),A.rows(),E.begin(),1);		
	cout << "A copied to E " << A(1,1) << " " << E[1] << endl;
	
	temp=0.0;
	temp=determinant(C.begin(),C.end(),C.rows(),temp);		
	cout << "det C "<< temp << endl;
	cout << "matrix C "<< C << endl;

	matrix<double> F(4,4);
	matrix<double> G(4,4);
	matrix<double> H(4,4);
	matrix<double> I(4,4);
	matrix<double> J(4,4);
	matrix<double> K(3,3);
	int i,j;
	temp=0.0;
	for(i=1;i<5;i++)
	{
	    for(j=1;j<5;j++)
	    {
		temp += 1.0;
		F(i,j)=temp;
	    }
	}
	cout << "matrix F "<< F << endl;

	transpose(F.begin(),F.end(),F.rows(),G.begin());
	cout << "matrix G=FT "<< G << endl;

	matrix_matrixtranspose_product(F.begin(),F.end(),F.rows(),
					G.begin(),G.end(),G.rows(),
					H.begin());
	cout << "matrix FGT "<< H << endl;

	matrixtranspose_matrix_product(F.begin(),F.end(),F.rows(),
					G.begin(),G.end(),G.rows(),
					I.begin());
	cout << "matrix FTG "<< I << endl;

	matrix_product(G.begin(),G.end(),G.rows(),
			G.begin(),G.end(),G.rows(),
			J.begin());
	cout << "matrix GG "<< J << endl;

	matrix_minor(J.begin(),J.end(),J.rows(),K.begin(),1,3);
	cout << "matrix minor of GG "<< K << endl;

	matrix<double> L(4,4,1.3);
	temp=0.0;
	K(1,1)= 1.0;
	K(1,2)= 3.0;
	K(1,3)= 7.0;
	K(2,1)= 3.0;
	K(2,2)= 4.0;
	K(2,3)= 5.0;
	K(3,1)= 3.0;
	K(3,2)= 2.0;
	K(3,3)= 2.0;
	cout << "K "<< K << endl;
	temp=determinant(K.begin(),K.end(),K.rows(),temp);		
	cout << "det K "<< temp << endl;

	matrix<double> M(3,3,1.3);
	adjoint(K.begin(),K.end(),K.rows(),M.begin());		
	cout << "adjoint K "<< M << endl;
	inverse_square(K.begin(),K.end(),K.rows(),M.begin());		
	cout << "inverse K "<< M << endl;

	matrix<double> N(3,3,0.0);
	matrix_product(M.begin(),M.end(),M.rows(),
			K.begin(),K.end(),K.rows(),
			N.begin());
	cout << "invK.K "<< N << endl;

	matrix<double> U;
	numeric_vector<double> X;
	matrix<double> V;

	K(1,1)= 4.0; K(1,2)= 2.0; K(1,3)= 14.0;
	K(2,1)= 2.0; K(2,2)= 17.0; K(2,3)= -5.0;
	K(3,1)= 14.0; K(3,2)= -5.0; K(3,3)= 83.0;

	inverse_cholesky(K.begin(),K.end(),K.rows(),M.begin());		
	cout << "inverse K "<< M << endl;

	N *= 0.0;
	matrix_product(M.begin(),M.end(),M.rows(),
			K.begin(),K.end(),K.rows(),
			N.begin());
	cout << "invK.K "<< N << endl;

	K(1,1)= 0.0;
	K(1,2)= 1.0;
	K(1,3)= 0.0;
	K(2,1)= 1.0;
	K(2,2)= 0.0;
	K(2,3)= 0.0;
	K(3,1)= 0.0;
	K(3,2)= 0.0;
	K(3,3)= 0.0;
	numeric_vector<double> P(3);
	eigen_solution(K.begin(),K.end(),K.rows(),M.begin(),P.begin());
	cout << "eigenvectors " << M << endl;
	cout << "eigenvalues " << P << endl;


	

	assert(C(1,2)==C(2,1));
	cout << "multiply " << C(1,2) << endl;

    	return 0;
}