btl_numerics.h

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

        for(;j<KK;j++)    ran_u[j-LL] = u[j];
        fill_buffer(BUFFER_SIZE);
}

T operator() () 
{
        if(current_position == BUFFER_SIZE) fill_buffer(BUFFER_SIZE);
        return buffer[current_position++];
}

}; // end of rng for floating point numbers


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

template <>
class random_uniform<long>
{
private:

// The lagged subtractive generator x(i) = (x(i-KK) - x(i-LL))mod 2**MM 

	static const long MM = 1L<<30;   // modulus = 2**30
	static const int KK = 100;       // the long lag
	static const int LL = 37;        // the short lag
                                         // lag pairs (127, 30) and (258, 83) can also be used
	static const int TT = 70;        // guarantees many seperate streams of numbers 
					 // ~ 2**TT streams when TT + MM = KK
	static const int DK = 200;       // KK+ KK
	static const int KL = 63;        // KK - LL

	long ran_x[KK];		         // current state of random number generator

        int current_position;
        static const int BUFFER_SIZE = 1009;
        long buffer[BUFFER_SIZE];

void fill_buffer(int n)
{
	register int i,j;
	for (j=0;j<KK;j++)     buffer[j] = ran_x[j];
	for (;j<n;j++)         buffer[j] = ((buffer[j-KK] - buffer[j-LL])&(MM-1));  // sum mod MM-1 
	for (i=0;i<LL;i++,j++)  ran_x[i] = ((buffer[j-KK] - buffer[j-LL])&(MM-1));
	for (;i<KK;i++,j++)     ran_x[i] = ((buffer[j-KK] - ran_x[i-LL])&(MM-1));
        current_position = 0;
}

public:

random_uniform(long seed)
{ 
// add debug for n > 1073741821 i.e. (MM-3)
    register int t,j;
    long  x[DK-1]; 
    register long ss = ((seed+2)&(MM-2));

    for (j=0;j<KK;j++)
    {
        x[j] = ss;
        ss *= 2;
        if(ss>=MM) ss -= (MM-2);                
    }
    for(;j<DK-1;j++) x[j] = 0;
    x[1]++;
    ss = seed&(MM-1);                  // ss forced < MM, both s and t are counters
    t = TT - 1;

    while (t)
    { 
        for (j=KK-1;j>0;j--)  x[j+j] = x[j];
        for (j=DK-2;j>KK-LL;j-=2) x[DK-1-j] = (x[j]&(MM-2));
        for (j=DK-2;j>=KK;j--)
        {
            if ((x[j]%2) == 1)
            {
                x[j-KL] = ((x[j-KL] - x[j])&(MM-1));
                x[j-KK] = ((x[j-KK] - x[j])&(MM-1));
            }
	}
        if ((ss%2) == 1)                               // if s is odd
        {
            for (j=KK;j>0;j--) x[j] = x[j-1];
            x[0] = x[KK];
            if ((x[KK]%2) == 1) x[LL] = ((x[LL] - x[KK])&(MM-1)) ; //???
        }

         if (ss) ss /= 2;             // successive integer divides of seed value by 2 
	 else t--;		      // when 0 is reached countdown TT further cycles
     }
 
     for(j=0;j<LL;j++) ran_x[j+KL] = x[j];
     for(;j<KK;j++)    ran_x[j-LL] = x[j];
     fill_buffer(BUFFER_SIZE);
}


long operator() () 
{
        if(current_position == BUFFER_SIZE) fill_buffer(BUFFER_SIZE);
        return buffer[current_position++];
}

}; // end of rng for long integers


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

template <>
class random_uniform<int>
{
private:

// The lagged subtractive generator x(i) = (x(i-KK) - x(i-LL))mod 2**MM 

	static const unsigned int MM = 1U<<15;   // modulus = 2**15
	static const int KK = 100;       // the long lag
	static const int LL = 37;        // the short lag
                                         // lag pairs (127, 30) and (258, 83) can also be used
	static const int TT = 70;        // guarantees many seperate streams of numbers 
					 // ~ 2**TT streams when TT + MM = KK
	static const int DK = 200;       // KK + KK
	static const int KL = 63;        // KK - LL

	int ran_x[KK];		         // current state of random number generator

        int current_position;
        static const int BUFFER_SIZE = 1009;
        int buffer[BUFFER_SIZE];

void fill_buffer(int n)
{
	register int i,j;
	for (j=0;j<KK;j++)     buffer[j] = ran_x[j];
	for (;j<n;j++)         buffer[j] = ((buffer[j-KK] - buffer[j-LL])&(MM-1));  // sum mod MM-1 
	for (i=0;i<LL;i++,j++)  ran_x[i] = ((buffer[j-KK] - buffer[j-LL])&(MM-1));
	for (;i<KK;i++,j++)     ran_x[i] = ((buffer[j-KK] - ran_x[i-LL])&(MM-1));
        current_position = 0;
}

public:

random_uniform(int seed)
{ 
    register int t,j;
    int x[DK-1]; 
    register unsigned int ss = ((seed+2)&(MM-2));

    for (j=0;j<KK;j++)
    {
        x[j] = ss;
        ss *= 2;
        if(ss>=MM) ss -= (MM-2);                
    }
    for(;j<DK-1;j++) x[j] = 0;
    x[1]++;
    ss = seed&(MM-1);                  // ss forced < MM, both s and t are counters
    t = TT - 1;

    while (t)
    { 
        for (j=KK-1;j>0;j--)  x[j+j] = x[j];
        for (j=DK-2;j>KK-LL;j-=2) x[DK-1-j] = (x[j]&(MM-2));
        for (j=DK-2;j>=KK;j--)
        {
            if ((x[j]%2) == 1)
            {
                x[j-KL] = ((x[j-KL] - x[j])&(MM-1));
                x[j-KK] = ((x[j-KK] - x[j])&(MM-1));
            }
	}
        if ((ss%2) == 1)                               // if s is odd
        {
            for (j=KK;j>0;j--) x[j] = x[j-1];
            x[0] = x[KK];
            if ((x[KK]%2) == 1) x[LL] = ((x[LL] - x[KK])&(MM-1)) ; //???
        }

         if (ss) ss /= 2;             // successive integer divides of seed value by 2 
	 else t--;		      // when 0 is reached countdown TT further cycles
     }
 
     for(j=0;j<LL;j++) ran_x[j+KL] = x[j];
     for(;j<KK;j++)    ran_x[j-LL] = x[j];
     fill_buffer(BUFFER_SIZE);
}


int operator() () 
{
        if(current_position == BUFFER_SIZE) fill_buffer(BUFFER_SIZE);
        return buffer[current_position++];
}

}; // end of rng for integers

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

	    	/**#: [Description="Portable random number generator that produces float
                                    or double precision random numbers from a normal distibution
                                    using the ratio method of Kinderman & Monahan (see Knuth TAOCP)."] */

template <class T = BTL_REAL>
class random_normal : public random_uniform<T>
{
private:

        double m, sd, mult1, mult2, temp;
	T u,v;

public:

random_normal(long seed) : random_uniform<T>(seed)
{
        m = 0.0; sd = 1.0;
        mult1 = -0.5*exp(1.0); mult2 = sqrt(8.0*exp(-1.0))*sd;
}
random_normal(long seed, T mean, T standard_deviation) : random_uniform<T>(seed) 
{
        m = mean; sd = standard_deviation;
        mult1 = -0.5*exp(1.0); mult2 = sqrt(8.0*exp(-1.0)) * sd;
}

T operator() () 
{
        u = random_uniform<T>::next_number();
        v = random_uniform<T>::next_number() - 0.5;
        while (v*v > mult1*u*u*log(u))
        {
            u = random_uniform<T>::next_number();
            v = random_uniform<T>::next_number() - 0.5;
        }
        return m + mult2*v/u;
}

}; // end of rng for normally distributed floating point numbers
 

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


	    	/**#: [Description="Portable random number generator that produces float
                                    or double precision random numbers from an exponential distibution."] */

template <class T = BTL_REAL>
class random_exponential : public random_uniform<T>
{
private:

	double m;

public:

random_exponential(long seed) : random_uniform<T>(seed)
{
	m = 1.0; 
}
random_exponential(long seed, T mean) : random_uniform<T>(seed)
{
	m = mean;
}

T operator() () 
{
	return m * log(random_uniform<T>::next_number()); 
}

}; // end of rng for exponentially distributed floating point numbers
 
//.........................................................................................

_BTL_END_NAMESPACE

#endif