ownRandom.h

Go to the documentation of this file.

#ifndef _OWN_RANDOM_H_
#define _OWN_RANDOM_H_

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <vector>

class ownRandom{
                public:
                                float *normalRandomCache;
                                bool isOk;
                                int size;
                                
                                ownRandom(){
                                                isOk=false;
                                                size = 0;
                                                allocate(10000);
                                                srand( (unsigned)time( NULL ) ); 
                                }
                                ~ownRandom(){
                                
                                }
                                void allocate(int m_size){
                                                if(size != 0){ 
                                                                if(normalRandomCache)free(normalRandomCache);
                                                }
                                                 
                                                normalRandomCache = (float *) malloc(sizeof(float)*m_size);
                                                isOk=true;
                                                size = m_size;
                                }
                                
                                void fillCache(){
                                                if(size > 0 && normalRandomCache != NULL){
                                                                for(int i=0;i<size;i++){
                                                                                normalRandomCache[i] = normalRandom(); 
                                                                                isOk = true;     
                                                                }
                                                }else{
                                                                allocate(10000);
                                                                fillCache();
                                                }
                                }
                                float getCachedUniformRandom(){
                                                int rv;
                                                float scale;
                                                
                                                if(!isOk) fillCache();
                                                
                                                scale = (float) size / (float)RAND_MAX;
                                                rv = rand(); 
                                                rv = (int) ((float) rv * scale);
                                                if(rv >=size){
                                                                fprintf(stderr,"ERROR in getCachedUniformRandom(). rv=%d\n",rv);
                                                                rv = size-1; 
                                                }
                                                return normalRandomCache[rv];
                                }
                                
                                float uniformRandom(void){
                                                return ( ((float) rand() / (float) RAND_MAX));
                                }
                                
                                float normalRandom()
                                {
                                                static float gset;
                                                static int randomStored = 0;
                                                float fac,rsq,v1,v2;

                                                if(randomStored){
                                                                return gset;
                                                }
    
                                                do {
                                                                v1=2.0*uniformRandom()-1.0; //pick two uniform numbers in the square extending
                                                                                                                                                                                //from -1 to +1 in each direction, 
                                                                v2=2.0*uniformRandom()-1.0;
                                                                rsq=v1*v1+v2*v2;                //see if they are in the unit circle,
                                                } while (rsq >= 1.0 || rsq == 0.0); //if not try again
        
                                                fac=sqrt(-2.0*log(rsq)/rsq);
                                                gset=v1*fac; // can be used also as normal distributed random number!!
                                                return v2*fac;
                                }
                                
                                void covRandom(float& xRand, float& yRand, float cov[3]){
                                                float a,b,c;
                                                float eig1,eig2;
                                                float vx1,vy1,vx2,vy2;
                                                float d;
                                                float x,y;
                                                
                                                xRand = normalRandom();
                                                yRand = normalRandom();
                                                // calculate eigen values
                                                a = 1.0;
                                                b = -(cov[0]+cov[2]);
                                                c = cov[0]*cov[2] - cov[1]*cov[1];

                                                eig1 = (-b + sqrt(b*b-4*a*c))/(2*a);
                                                eig2 = (-b - sqrt(b*b-4*a*c))/(2*a);

                                                // calculate eigen vectors
                                                vx1 = cov[1];
                                                vy1 = eig1 - cov[0];
                                                d = sqrt(vx1*vx1+vy1+vy1);
                                                vx1 /= d;
                                                vy1 /= d;

                                                vx2 = cov[1];
                                                vy2 = eig2 - cov[0];
                                                d = sqrt(vx2*vx2+vy2+vy2);
                                                vx2 /= d;
                                                vy2 /= d;
        
                                                // Transform the random variable according to the covariance
                                                x = xRand * sqrt(eig1);
                                                y = yRand * sqrt(eig2);
        
                                                xRand = x * vx1 + y * vx2;
                                                yRand = x * vy1 + y * vy2;
                                }
};
#endif