example4.cpp

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

//
// This file contains some example code, which may be used to call 
// the Fast Fourier Transform Class.
// This code is not meant to be usefull in itself, but is provided as 
// an example of how the FFT class may be used.
//
// Copyright (C) 1997-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: Will Pitt, Mary Steven & Breen Sweeney, Mark Williams 
// 
///////////////////////////////////////////////////////////////////////////
// 
// Brief Description of Code:
// 
// This routine is a program which uses the Fast Fourier Transform class 
// (btl_fourier class) to carry out a One-dimensional, Two-dimensional, 
// and a Three-dimensional Fast Fourier Transform.
// This program is not meant to be used, and carries out no useful function
// itself. It is provided as an example of how the routines can be called
// from other C++ programs.
//
// For further details of the fourier_transform algorithms see the 
// documentation for the individual classes.
///////////////////////////////////////////////////////////////////////////

#include <string>
#include <complex>
#include <fstream>
using namespace std;

#include "btl_fourier.h"
#include "btl_numerics.h"
using namespace btl;


void TestResults(const string& testName, const complex_type& error)
{
    cout << " Ideally difference should be zero "
         << ">>>>> Difference = "
         << error
         << " : "
         << testName
         << " <<<<<\n"
         << endl;
}

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

template<class RandomAccessIterator>
void FillComplexContainer(RandomAccessIterator start,
                          const unsigned int size)
{
    // Set different random number seed each time this function is called.
    static unsigned int seed = 0;
    seed++;
    srand(seed);

    for (unsigned int i=0; i<size; i++)
    	start[i] = complex_type((rand()/10000.0),(rand()/10000.0));
}

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

template<class RandomAccessIterator >
complex<double> max_mismatch(RandomAccessIterator firsta,
                             RandomAccessIterator lasta,
                             RandomAccessIterator firstb,
                             complex<double> result)
{
    double maxError = 0.0, error = 0.0;
    for (; firsta<lasta; firsta++, firstb++) 
    {
         error = norm(*firsta - *firstb);
         if (error > maxError){
             result = *firsta - *firstb;
             maxError = error;
         }
    }           
    return result;
}

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

template <class RandomAccessIterator>
complex_type OneDimForwardBackwardTest(RandomAccessIterator testarr,
    	    	         	       RandomAccessIterator intermediate,
    	    	         	       RandomAccessIterator copy,
    	    	         	       const unsigned int size)
{    	    	     	 

    complex_type ZERO = 0.0;
    	// Transform in forward direction
    	//
    fourier_transform(testarr, intermediate, size, true);
   
    	// Transform back - overwriting the origin data
    	//
    fourier_transform(intermediate, testarr, size, false);
    
    	// Compare results with original
    	//
    return( max_mismatch(testarr,(testarr+size),copy,ZERO) );
}

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

template <class RandomAccessIterator>
complex_type TwoDimForwardBackwardTest(RandomAccessIterator testarr,
    	    	         	       RandomAccessIterator intermediate,
    	    	         	       RandomAccessIterator copy,
    	    	         	       const unsigned int nrows,
    	    	         	       const unsigned int ncols)
{    	    	     	 

    complex_type ZERO = 0.0;    
    	// Transform in forward direction
    	//
    fourier_transform(testarr, intermediate, nrows, ncols, true);
   
    	// Transform back - overwriting the origin data
    	//
    fourier_transform(intermediate, testarr, nrows, ncols, false);
    
    	// Compare results with original
    	//
    unsigned int size = nrows*ncols;    	
    return( max_mismatch(testarr,(testarr+size),copy,ZERO) );
}

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

template <class RandomAccessIterator>
complex_type ThreeDimForwardBackwardTest(RandomAccessIterator testarr,
    	    	         	         RandomAccessIterator intermediate,
    	    	         	         RandomAccessIterator copy,
    	    	         	         const unsigned int nrows,
    	    	         	         const unsigned int ncols,
    	    	         	         const unsigned int nlayers)
{    	    	     	 
    
    complex_type ZERO = 0.0;    
    	// Transform in forward direction
    	//
    fourier_transform(testarr, intermediate, nrows, ncols, nlayers, true);
   
    	// Transform back - overwriting the origin data
    	//
    fourier_transform(intermediate, testarr, nrows, ncols, nlayers, false);
    
    	// Compare results with original
    	//
    unsigned int size = nrows*ncols*nlayers;    	
    return( max_mismatch(testarr,(testarr+size),copy,ZERO)  );
}

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

void Check1DVectorInput(const unsigned int size)
{
   
    	// Create test data in array form
    	//
    vector<complex_type> testarr(size);
    FillComplexContainer(testarr.begin(),size);

    	// Copy original for comparison purposes
    	//
    vector<complex_type> copy = testarr;

    	// Create array to store intermediate results;
    	//
    vector<complex_type> intermediate(size);

    	// Do test
    	//
    complex_type error = OneDimForwardBackwardTest(testarr.begin(),
                                             	   intermediate.begin(),
                                             	   copy.begin(),
                                             	   size);
                                             
    TestResults("Results for vector of complex",error);
}

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

void Check2DVectorInput(const unsigned int nrows,
                            const unsigned int ncols)
{    
    	// Create test data in array form
    	//
    unsigned int size = nrows*ncols;
    vector<complex_type> testarr(size);
    FillComplexContainer(testarr.begin(),size);

    	// Copy original for comparison purposes
    	//
    vector<complex_type> copy = testarr;

    	// Create array to store intermediate results;
    	//
    vector<complex_type> intermediate(size);

    	// Do test
    	//
    complex_type error = TwoDimForwardBackwardTest(testarr.begin(),
                                             	   intermediate.begin(),
                                             	   copy.begin(),
                                             	   nrows,
                                             	   ncols);
                                             
    TestResults("Results for vector containing 2D data.",error);
}

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

void Check3DVectorInput(const unsigned int nrows,
                        const unsigned int ncols,
                        const unsigned int nlayers)
{
    	// Create test data in array form
    	//
    unsigned int size = nrows*ncols*nlayers;
    vector<complex_type> testarr(size);
    FillComplexContainer(testarr.begin(),size);

    	// Copy original for comparison purposes
    	//
    vector<complex_type> copy = testarr;

    	// Create array to store intermediate results;
    	//
    vector<complex_type> intermediate(size);

    	// Do test
    	//
    complex_type error = ThreeDimForwardBackwardTest(testarr.begin(),
                                             	     intermediate.begin(),
                                             	     copy.begin(),
                                             	     nrows,
                                             	     ncols,
                                             	     nlayers);
                                             
    TestResults("Results for vector containing 3D data.",error);
}


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

void Check1DArrayInput(const unsigned int size)
{
    	// Create test data in array form
    	//

    complex_type *testarr = new complex_type[size];
    FillComplexContainer(testarr,size);

    	// Copy original for comparison purposes
    	//
    complex_type *copy = new complex_type[size];
    for (unsigned int i=0; i<size; i++) copy[i] = testarr[i];

    	// Create array to store intermediate results;
    	//
    complex_type *intermediate = new complex_type[size];

    	// Do test
    	//
    complex_type error = OneDimForwardBackwardTest(testarr,
                                                   intermediate,
                                                   copy,
                                                   size);
                                                   
    TestResults("Results for array of complex",error);
     
    	// delete arrays
    	//
    delete[] testarr;
    delete[] copy;
    delete[] intermediate; 
}

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

void Check2DArrayInput(const unsigned int nrows, const unsigned int ncols)
{
    	// Create test data in array form
    	//
    unsigned int size = nrows*ncols;
    complex_type *testarr = new complex_type[size];
    FillComplexContainer(testarr,size);

    	// Copy original for comparison purposes
    	//
    complex_type *copy = new complex_type[size];
    for (unsigned int i=0; i<size; i++) copy[i] = testarr[i];

    	// Create array to store intermediate results;
    	//
    complex_type *intermediate = new complex_type[size];

    	// Do test
    	//
    complex_type error = TwoDimForwardBackwardTest(testarr,
                                                   intermediate,
                                                   copy,
                                                   nrows,
                                                   ncols);
                                                   
    TestResults("Results for array containing 2D data.",error);
     
    	// delete arrays
    	//
    delete[] testarr;
    delete[] copy;
    delete[] intermediate; 
}

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

void Check3DArrayInput(const unsigned int nrows, 
    	    	       const unsigned int ncols,
    	    	       const unsigned int nlayers)
{
    	// Create test data in array form
    	//
    unsigned int size = nrows*ncols*nlayers;
    complex_type *testarr = new complex_type[size];
    FillComplexContainer(testarr,size);

    	// Copy original for comparison purposes
    	//
    complex_type *copy = new complex_type[size];
    for (unsigned int i=0; i<size; i++) copy[i] = testarr[i];

    	// Create array to store intermediate results;
    	//
    complex_type *intermediate = new complex_type[size];

    	// Do test
    	//
    complex_type error = ThreeDimForwardBackwardTest(testarr,
                                                     intermediate,
                                                     copy,
                                                     nrows,
                                                     ncols,
                                                     nlayers);
                                                   
    TestResults("Results for array containing 3D data.",error);
     
    	// delete arrays
    	//
    delete[] testarr;
    delete[] copy;
    delete[] intermediate; 
}

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

void Series_1Dtest()
{

    typedef vector<complex_type> Container;
    unsigned int dataSize = 1024;
    Container sin_freq_10(dataSize),
    	      sin_freq_20(dataSize),
    	      sin_freq_30(dataSize),
    	      sin_freq_sum(dataSize);

    double step=0.0, stepSize = 2.0 * M_PI / dataSize;
    for (unsigned int i=0; i<dataSize; i++, step += stepSize)
    //for (int i=0; i<dataSize; i++, step += stepSize)
    {
    	sin_freq_10[i] = complex_type(sin(10.0*step));
    	sin_freq_20[i] = complex_type(sin(20.0*step));
    	sin_freq_30[i] = complex_type(sin(30.0*step));
    	sin_freq_sum[i] = sin_freq_10[i] + sin_freq_20[i] + sin_freq_30[i];
    }
    
    ofstream fin;
    fin.open("fft1D.in");

    // check that file can be opened
    if (!fin)
    {
        cerr << "\n\nError creating fft1D.in" << endl;
        return;
    }
    
    
    // Please note the ii,iii etc variables are for the benefit of the SGI CC
    // compiler 
    
    for (unsigned int ii=0; ii<dataSize; ii++)          // DEBUG was i<datasize
   
    //for (int ii=0; i<dataSize; ii++)
    	
    fin << real(sin_freq_sum[ii]) << '\n';
        
    Container output_10(dataSize),
    	      output_20(dataSize),
    	      output_30(dataSize),
    	      output_sum(dataSize);
    	      
    fourier_transform(sin_freq_10.begin(),output_10.begin(), dataSize, true);
    fourier_transform(sin_freq_20.begin(),output_20.begin(), dataSize, true);
    fourier_transform(sin_freq_30.begin(),output_30.begin(), dataSize, true);
    fourier_transform(sin_freq_sum.begin(),output_sum.begin(), dataSize, true);
   
    vector<double> powerSpectrum_sum(dataSize);
    Container FFT_sum(dataSize);
    double theoretical_result, diff, maxdiff = 0.0;

    for (unsigned int iii=0; iii<dataSize; iii++)
    {
    	FFT_sum[iii] = output_10[iii]+output_20[iii]+output_30[iii];
    	theoretical_result  = norm(FFT_sum[iii]);
    	powerSpectrum_sum[iii] = norm(output_sum[iii]);
    	diff = fabs(theoretical_result - powerSpectrum_sum[iii]);
    	if (diff > maxdiff) maxdiff = diff;
    }

    complex_type ZERO = 0.0;
    cout << "1D Sin series test: max difference from ideal result = " << maxdiff
         << "\n";