random_tools.hpp

Go to the documentation of this file.

#ifndef RANDOM_TOOLS_HPP
#define RANDOM_TOOLS_HPP
#include <iostream>
#include <cstdio>

#include "gambit/Utils/threadsafe_rng.hpp"

template <typename T>
T *matrix(const int xN)
{
     T *temp = new T[xN];
     
     return temp;
}

template <typename T>
T **matrix(const int xN, const int yN)
{
     T **temp = new T*[xN];
     
     for (int i = 0; i < xN; i++)
          temp[i] = new T[yN];
     
     return temp;
}

template <typename T>
T ***matrix(const int xN, const int yN, const int zN)
{
     T ***temp = new T**[xN];
     
     for (int i = 0; i < xN; i++)
     {
          temp[i] = new T*[yN];
          for (int j = 0; j < yN; j++)
               temp[i][j] = new T[zN];
     }
     return temp;
}

template <typename T>
T *matrix(const int xN, T in)
{
        T *temp = new T[xN];

        for (int i = 0; i < xN; i++)
                temp[i] = in;

        return temp;
}

template <typename T>
T **matrix(const int xN, const int yN, T in)
{
        T **temp = new T*[xN];
     
        for (int i = 0; i < xN; i++)
        {
                temp[i] = new T[yN];
                for (int j = 0; j < yN; j++)
                        temp[i][j] = in;
        }
     
        return temp;
}

template <typename T>
T ***matrix(const int xN, const int yN, const int zN, T in)
{
        T ***temp = new T**[xN];
     
        for (int i = 0; i < xN; i++)
        {
                temp[i] = new T*[yN];
                for (int j = 0; j < yN; j++)
                {
                        temp[i][j] = new T[zN];
                        for (int k = 0; k < zN; k++)
                                temp[i][j][k] = in;
                }
        }
        return temp;
}

template <typename T>
                T *matrix(const int xN, T *in)
{
        T *temp = new T[xN];
        for (int i = 0; i < xN; i++)
                temp[i] = in[i];
     
        return temp;
}

template <typename T>
                T **matrix(const int xN, const int yN, T **in)
{
        T **temp = new T*[xN];
     
        for (int i = 0; i < xN; i++)
        {
                temp[i] = new T[yN];
                for (int j = 0; j < yN; j++)
                        temp[i][j] = in[i][j];
        }
     
        return temp;
}

template <typename T>
                T ***matrix(const int xN, const int yN, const int zN, T ***in)
{
        T ***temp = new T**[xN];
     
        for (int i = 0; i < xN; i++)
        {
                temp[i] = new T*[yN];
                for (int j = 0; j < yN; j++)
                {
                        temp[i][j] = new T[zN];
                        for (int k = 0; k < zN; k++)
                                temp[i][j][k] = in[i][j][k];
                }
        }
        return temp;
}


template <typename T>
void del(T *temp)
{
     delete[] temp;
     temp = NULL;
}

template <typename T>
void del(T **temp, int xN)
{
     for (int i = 0; i < xN; i++)
          delete[] temp[i];
     
     delete[] temp;
     temp = NULL;
}

template <typename T>
void del(T ***temp, int xN, int yN)
{
     for (int i = 0; i < xN; i++)
     {
          for (int j = 0; j < yN; j++)
               delete[] temp[i][j];
          delete[] temp[i];
     }
     
     delete[] temp;
     temp = NULL;
}

template<class T>
inline const T SQ(const T a) {return a*a;}

template<class T>
inline const T SQR(const T a) {return a*a;}

class Cholesky
{
        private:
                double **el;
                
        protected:
                int num;
        
                
        public:
                Cholesky(const int nin) : num(nin)
                {
                        el = matrix <double> (num, num);
                }
                
                Cholesky(double **a, const int nin) : num(nin)
                {
                        el = matrix <double> (num, num);
                        double sum = 0;
                        int i, j, k;
                        
                        for (i = 0; i < num; i++)
                                for (j = 0; j < num; j++)
                                        el[i][j] = a[i][j];

                        for (i = 0; i < num; i++)
                        {
                                for (j = i; j < num; j++)
                                {
                                        for(sum = el[i][j], k = i - 1; k >= 0; k--)
                                                sum -= el[i][k]*el[j][k];
                                        if(i ==j)
                                        {
                                                if(sum <= 0.0)
                                                {
                                                        std::cout << "Cholesky failed " << sum << std::endl;
                                                        getchar();
                                                }
                                                el[i][i] = sqrt(sum);
                                        }
                                        else
                                                el[j][i] = sum/el[i][i];
                                }
                        }
                        for (i = 0; i < num; i++)
                                for (j = 0; j < i; j++)
                                        el[j][i] = 0.0;
                }
                
                bool EnterMatM(double **a, const int min)
                {
                        double sum = 0;
                        int i, j, k;
                        //ofstream out1("cov.dat");
                        for (i = 0; i < num; i++)
                        {
                                for (j = 0; j < num; j++)
                                {
                                        el[i][j] = a[i][j];
                                        //out1 << el[i][j] << "   " << endl;
                                }
                                //out1 << endl;
                        }
                        
                        if (num >= min)
                        {
                                k = min -1;
                                int l = num - min + 1;
                                for (i = 0; i < l; i++)
                                {
                                        for (j = k+i; j < num; j++)
                                                el[i][j] = el[j][i] = 0.0;
                                }
                        }
                        
                        for (i = 0; i < num; i++)
                        {
                                for (j = i; j < num; j++)
                                {
                                        for(sum = el[i][j], k = i - 1; k >= 0; k--)
                                                sum -= el[i][k]*el[j][k];
                                        if(i ==j)
                                        {
                                                if(sum <= 0.0)
                                                        return true;
                                                el[i][i] = sqrt(sum);
                                        }
                                        else
                                                el[j][i] = sum/el[i][i];
                                }
                        }
                        for (i = 0; i < num; i++)
                                for (j = 0; j < i; j++)
                                        el[j][i] = 0.0;
                                
//                                              ofstream out2("cov2.dat");
//                      double **ct = matrix <double> (n, n, 0.0);
//                      for (i = 0; i < num; i++)
//                              for (j = 0; j < n; j++)
//                                      for (k = 0; k < num; k++)
//                                              ct[i][j] += el[i][k]*el[j][k];
//                      for (i = 0; i < num; i++)
//                      {
//                              for (j = 0; j < num; j++)
//                                      out2 << ct[i][j] << "   ";
//                              out2 << endl;
//                      }
//                      cout << "finished" << endl;
//                      getchar();
                                
                        return false;
                }
                
                bool EnterMat(double **a)
                {
                        double sum = 0;
                        int i, j, k;
                        
                        for (i = 0; i < num; i++)
                        {
                                for (j = 0; j < num; j++)
                                {
                                        el[i][j] = a[i][j];
                                }
                        }
                        
                        for (i = 0; i < num; i++)
                        {
                                for (j = i; j < num; j++)
                                {
                                        for(sum = el[i][j], k = i - 1; k >= 0; k--)
                                                sum -= el[i][k]*el[j][k];
                                        if(i ==j)
                                        {
                                                if(sum <= 0.0)
                                                        return false;
                                                el[i][i] = sqrt(sum);
                                        }
                                        else
                                                el[j][i] = sum/el[i][i];
                                }
                        }
                        for (i = 0; i < num; i++)
                                for (j = 0; j < i; j++)
                                        el[j][i] = 0.0;
                                
                        return true;
                }
                
                void EnterMat(double **a, int nin)
                {
                        del <double> (el, num);
                        num = nin;
                        el = matrix <double> (num, num);
                        EnterMat(a);
                }
                
                void ElMult (double *y, double *b)
                {
                        int i, j;
                        for(i = 0; i < num; i++)
                        {
                                b[i] = 0.0;
                                for (j = 0; j <= i; j++)
                                        b[i] += el[i][j]*y[j];
                        }
                }
                
                void ElMult (double *y)
                {
                        int i, j;
                        double *b = new double[num];
                        for(i = 0; i < num; i++)
                        {
                                b[i] = 0.0;
                                for (j = 0; j <= i; j++)
                                {
                                        b[i] += el[i][j]*y[j];
                                }
                        }
                        for (i = 0; i < num; i++)
                        {
                                y[i] = b[i];
                        }
                        delete[] b;
                }
                
                void Solve (double *b, double *x)
                {
                        int i, k;
                        double sum;
                        
                        for (i = 0; i < num; i++)
                        {
                                for (sum = b[i], k=i-1; k >=0; k--)
                                        sum -= el[i][k]*x[k];
                                x[i]=sum/el[i][i];
                        }
                        
                        for (i = num-1; i >=0; i--)
                        {
                                for (sum = x[i], k=i+1; k<num; k++)
                                        sum -= el[k][i]*x[k];
                                x[i]=sum/el[i][i];
                        }
                }
                
                double Square(double *y, double *y0)
                {
                        int i, j;
                        double sum;
                        double *x = new double[num];
                        
                        for (i = 0; i < num; i++)
                        {
                                for (sum = (y[i]-y0[i]), j=0; j < i; j++)
                                        sum -= el[i][j]*x[j];
                                x[i]=sum/el[i][i];
                        }
                        
                        sum = 0.0;
                        for (i = 0; i < num; i++)
                                sum += x[i]*x[i];
                        
                        delete[] x;
                        
                        return sum;
                }
                
                double Square(double *y, double *y0, int *map)
                {
                        int i, j;
                        double sum;
                        double *x = new double[num];
                        
                        for (i = 0; i < num; i++)
                        {
                                for (sum = (y[map[i]]-y0[i]), j=0; j < i; j++)
                                        sum -= el[i][j]*x[j];
                                x[i]=sum/el[i][i];
                        }
                        
                        sum = 0.0;
                        for (i = 0; i < num; i++)
                                sum += x[i]*x[i];
                        
                        delete[] x;
                        
                        return sum;
                }
                
                void Inverse(double **ainv)
                {
                        double sum;
                        
                        for (int i = 0; i < num; i++)
                                for (int j = 0; j <= i; j++)
                                {
                                        sum = (i == j ? 1.0 : 0.0);
                                        for (int k = i-1; k >= j; k--)
                                                sum -= el[i][k]*ainv[j][k];
                                        ainv[j][i] = sum/el[i][i];
                                }

                        for (int i = num - 1; i >= 0; i--)
                                for (int j = 0; j <= i; j++)
                                {
                                        sum = (i < j ? 0.0 : ainv[j][i]);
                                        for (int k = i + 1; k < num; k++)
                                                sum -= el[k][i]*ainv[j][k];
                                        ainv[i][j] = ainv[j][i] = sum/el[i][i];
                                }
                }

                double DetSqrt()
                {
                        double temp = 1.0;
                        for (int i = 0; i < num; i++)
                                temp *= el[i][i];
                        return temp;
                }
                ~Cholesky()
                {
                        del <double> (el, num);
                }
};

class Ran_old
{
        private:
                unsigned long long int u, v, w;
                
        public:
                Ran_old(unsigned long long int j) : v(4101842887655102017LL), w(1)
                {
                        u = j ^ v; int64();
                        v = u; int64(); int64();
                        w = v; int64(); int64();
                }
                inline unsigned long long int int64()
                {
                        u = u * 2862933555777941757LL + 7046029254386353087LL;
                        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
                        w = 4294957665U*(w & 0xffffffff) + (w >> 32);
                        unsigned long long int x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
                        return (x + v) ^ w;
                }
                inline double Doub(){return 5.42101086242752217E-20 * int64();}
                inline unsigned int int32(){return (unsigned int)int64();}
};

class Ran
{
public:
        Ran(unsigned long long int){}
        inline double Doub(){return Gambit::Random::draw();}
};

class ExponDev : public Ran
{
        private:
                double beta;
        public:
                ExponDev(double din, unsigned long long ii) : Ran(ii), beta(din){}
                double Dev()
                {
                        double u;
                        do
                        {
                                u = Doub();
                        }
                        while(u == 0);
                        
                        return -log(u)/beta;
                }
};

class NormalDev : public Ran
{
        private:
                double mu, sig;
                
        public:
                NormalDev(double mmu, double ssig, unsigned long long i) : Ran(i), mu(mmu), sig(ssig){}
                double  Dev()
                {
                        double u, v, x, y, q;
                        do
                        {
                                u = Doub();
                                v = 1.7156*(Doub() - 0.5);
                                x = u - 0.449871;
                                y = fabs(v) + 0.386595;
                                q = x*x + y*(0.19600*y-0.25472*x);
                        }
                        while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
                        
                        return mu + sig*v/u;
                }
};

class BasicDevs : public Ran
{
        private:
                
        public:
                BasicDevs(unsigned long long i) : Ran(i) {}

                double  Dev()
                {
                        double u, v, x, y, q;
                        do
                        {
                                u = Doub();
                                v = 1.7156*(Doub() - 0.5);
                                x = u - 0.449871;
                                y = fabs(v) + 0.386595;
                                q = x*x + y*(0.19600*y-0.25472*x);
                        }
                        while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
                        
                        return v/u;
                }
                
                double ExpDev()
                {
                        double u;
                        do
                        {
                                u = Doub();
                        }
                        while(u == 0);
                        
                        return -log(u);
                }
};

class MultiNormalDev : public Ran, public Cholesky
{
        private:
                int mm;
                double f;
                double *spt;
                
        public:
                MultiNormalDev (double **vvar, double fin, unsigned long long int j, int nin) : Ran(j), Cholesky(vvar, nin), mm(nin), f(fin)
                {
                        spt = matrix <double> (mm);
                }
                void Dev(double *pt, double *mean)
                {
                        double u, v, x, y, q;
                        for (int i = 0; i < mm; i++)
                        {
                                do
                                {
                                        u = Doub();
                                        v = 1.7156*(Doub() - 0.5);
                                        x = u - 0.449871;
                                        y = fabs(v) + 0.386595;
                                        q = x*x + y*(0.19600*y-0.25472*x);
                                }
                                while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
                                spt[i] = v/u;
                        }
                        ElMult(spt, pt);
                        for (int i = 0; i < mm; i++)
                                pt[i] = f*pt[i] + mean[i];
                }
                void Dev(double **cvar, double *pt, double *mean)
                {
                        double u, v, x, y, q;
                        for (int i = 0; i < mm; i++)
                        {
                                do
                                {
                                        u = Doub();
                                        v = 1.7156*(Doub() - 0.5);
                                        x = u - 0.449871;
                                        y = fabs(v) + 0.386595;
                                        q = x*x + y*(0.19600*y-0.25472*x);
                                }
                                while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
                                spt[i] = v/u;
                        }
                        EnterMat(cvar);
                        ElMult(spt, pt);
                        for (int i = 0; i < mm; i++)
                                pt[i] = f*pt[i] + mean[i];
                }
                ~MultiNormalDev()
                {
                        del <double> (spt);
                }
};

class AdvanceDevs : public BasicDevs, public Cholesky
{
        private:
                double fac;
                
        public:
                AdvanceDevs(int nin, double din, unsigned long long iin) :  BasicDevs(iin), Cholesky(nin), fac(din)
                {
                }
                
                AdvanceDevs(double **vvar, const int nin, double din, unsigned long long iin) :  BasicDevs(iin), Cholesky(vvar, nin), fac(din)
                {
                }
                
                void MultiDev(double *pin, double *p0)
                {
                        int i;
                        double dist = 0.0;

                        double vec[num];
                        double norm = 0.0;
                        for (i = 0; i < num; i++)
                        {
                                vec[i] = Dev();
                                norm += vec[i]*vec[i];
                        }
                        norm = sqrt(norm);
                        
                        if(Doub() < 0.33)
                                dist = ExpDev();
                        else
                        {
                                double temp;
                                for (i = 0; i < 2; i++)
                                {
                                        temp = Dev();
                                        dist += temp*temp;
                                }
                                dist = sqrt(dist/2.0);
                        }

                        ElMult(vec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fac*dist*pin[i]/norm + p0[i];
                        }

                        return;
                }
                
                void MultiDev(double **cvar, double *pin, double *p0)
                {
                        int i;
                        double dist = 0.0;

                        double vec[num];
                        double norm = 0.0;
                        for (i = 0; i < num; i++)
                        {
                                vec[i] = Dev();
                                norm += vec[i]*vec[i];
                        }
                        norm = sqrt(norm);
                        
                        if(Doub() < 0.33)
                                dist = ExpDev();
                        else
                        {
                                double temp;
                                for (i = 0; i < 2; i++)
                                {
                                        temp = Dev();
                                        dist += temp*temp;
                                }
                                dist = sqrt(dist/2.0);
                        }

                        EnterMat(cvar);
                        ElMult(vec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fac*dist*pin[i]/norm + p0[i];
                        }

                        return;
                }
                
                void EllipseDev(double *pin, double *p0, double fin)
                {
                        int i;
                        double dist = 0.0;
                        
                        double vec[num];
                        double norm = 0.0;
                        for (i = 0; i < num; i++)
                        {
                                vec[i] = Dev();
                                norm += vec[i]*vec[i];
                        }
                        norm = sqrt(norm);
                        
                        dist = pow(Doub(), 1.0/num);
                        
                        ElMult(vec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fin*dist*pin[i] + p0[i];
                        }

                        return;
                }
                
                void EllipseDev(double **cvar, double *pin, double *p0, double fin)
                {
                        int i;
                        double dist = 0.0;

                        double vec[num];
                        double norm = 0.0;
                        for (i = 0; i < num; i++)
                        {
                                vec[i] = Dev();
                                norm += vec[i]*vec[i];
                        }
                        norm = sqrt(norm);
                        
                        dist = pow(Doub(), 1.0/num);
                        
                        EnterMat(cvar);
                        ElMult(vec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fin*dist*pin[i] + p0[i];
                        }

                        return;
                }
};

class RandomPlane : public AdvanceDevs
{
private:
        double **rotVec;
        double **currentVec;
        double **endVec;
        double **endEndVec;
        double **sEndVec;
        int proj, extra;
        double alim, alimt;
        double sqrtAlim, powAlim, ratioAlim, ratioAlimt;
        
public:
        RandomPlane(const int projin, const int nin, const double din, const double alim, const double alimt, unsigned long long iin) : AdvanceDevs(nin, din, iin), proj(projin), alim(alim), alimt(alimt)
        {
                rotVec = matrix <double> (nin, nin);
                RandRot();
                if (proj > num) proj = num;
                extra = num % proj;
                endVec = currentVec + num - extra;
                endEndVec = currentVec + num;
                sEndVec = currentVec + num - proj;
                sqrtAlim = sqrt(alim);
                powAlim = pow(alim, 0.5 - proj);
                double vol1 = 2.0*(sqrt(alim) - 1.0);
                double vol2 = (1.0 - pow(alim, 0.5 - proj))/(proj - 0.5);
                ratioAlim = vol1/(vol1 + vol2);
                ratioAlimt = (alimt - 1.0)/(2.0*alimt + proj - 2.0);
        }
        
        double WalkDev()
        {
                if (Doub() <= ratioAlim)
                        return SQ(1.0 + Doub()*(sqrtAlim - 1.0));
                else
                        return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5));
        }
        
        double TransDev()
        {
                if (Doub() < ratioAlimt)
                {
                        return -pow(Doub(), 1.0/(alimt + proj - 1.0));
                }
                else
                {
                        return -pow(Doub(), 1.0/(1.0 - alimt));
                }
        }
        
        double KWalkDev()
        {
                if (Doub() < 0.5)
                {
                        if (Doub() <= ratioAlim)
                        {
                                return SQ(1.0 + Doub()*(sqrtAlim - 1.0));
                        }
                        else
                        {
                                return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5));
                        }
                }
                else
                {
                        if (Doub() < ratioAlimt)
                        {
                                return -pow(Doub(), 1.0/(alimt + proj - 1.0));
                        }
                        else
                        {
                                return -pow(Doub(), 1.0/(1.0 - alimt));
                        }
                }
        }
        
        bool KWalkDev(double &Z)
        {
                if (Doub() < 0.5)
                {
                        Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim));
                        
                        return pow(Z, proj - 1) > Doub();
                }
                else
                {
                        if (Doub() < (alimt - 1.0)/(2.0*alimt))
                        {
                                Z = -pow(Doub(), 1.0/(alimt + 1.0));
                        }
                        else
                        {
                                Z = -pow(Doub(), 1.0/(1.0 - alimt));
                        }
                        
                        return pow(Z, proj - 2) > Doub();
                }
        }
        
        double WalkDev(double *ptrOut, double *ptr, double *ptr0)
        {
                double Z;
   
                Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim));
                
                Mult(ptrOut, ptr, ptr0, Z);
                
                return (proj - 1.0)*log(Z);
        }
        
        double TransDev(double *ptrOut, double *ptr, double *ptr0)
        {
                double Z;
                
                if (Doub() < (alimt - 1.0)/(2.0*alimt))
                {
                        Z = -pow(Doub(), 1.0/(alimt + 1.0));
                }
                else
                {
                        Z = -pow(Doub(), 1.0/(1.0 - alimt));
                }
                
                Mult(ptrOut, ptr, ptr0, Z);
                
                return (proj - 1.0)*log(-Z);
        }
        
        void Mult(double *ptrOut, double *ptr, double *ptr0, const double Z)
        {
                double z[proj];
                for (int i = 0; i < proj; i++)
                {
                        z[i] = 0.0;
                        for (int j = 0; j < num; j++)
                        {
                                z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j];
                        }
                }
                
                for (int j = 0; j < num; j++)
                {
                        ptrOut[j] = ptr[j];
                        for (int i = 0; i < proj; i++)
                        {
                                ptrOut[j] += z[i]*currentVec[i][j];
                        }
                }

                currentVec++;// += proj;
                if (currentVec == sEndVec)
                {
                        for (double **vec = currentVec, **vec2 = rotVec; vec != endEndVec; vec++, vec2++)
                        {
                                double *temp = *vec2;
                                *vec2 = *vec;
                                *vec = temp;
                        }
                        RandRot(extra);
                }
        }
        
        void Mult2(double *ptrOut, double *ptr, double *ptr0, const double Z)
        {
                RandRot(0, proj);
                
                double z[proj];
                for (int i = 0; i < proj; i++)
                {
                        z[i] = 0.0;
                        for (int j = 0; j < num; j++)
                        {
                                z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j];
                        }
                }
                
                for (int j = 0; j < num; j++)
                {
                        ptrOut[j] = ptr[j];
                        for (int i = 0; i < proj; i++)
                        {
                                ptrOut[j] += z[i]*currentVec[i][j];
                        }
                }
        }
        
        void RandRot(const int start = 0)
        {
                double temp;
                double vec[num];
                int i, j, k;
                        
                for (i = start; i < num; i++)
                {
                        for (j = 0; j < num; j++)
                                vec[j] = Dev();
                        for (j = 0; j < i; j++)
                        {
                                temp = 0.0;
                                for (k = 0; k < num; k++)
                                        temp += vec[k]*rotVec[j][k];
                                for (k = 0; k < num; k++)
                                        vec[k] -= temp*rotVec[j][k];
                        }
                        temp = 0.0;
                        for (j = 0; j < num; j++)
                                temp += vec[j]*vec[j];
                        temp = sqrt(temp);
                        for (j = 0; j < num; j++)
                                rotVec[i][j] = vec[j]/temp;
                }
                
                currentVec = rotVec;
        }
        
        void RandRot(const int start, const int end)
        {
                double temp;
                double vec[num];
                int i, j, k;
                        
                for (i = start; i < end; i++)
                {
                        for (j = 0; j < num; j++)
                                vec[j] = Dev();
                        for (j = 0; j < i; j++)
                        {
                                temp = 0.0;
                                for (k = 0; k < num; k++)
                                        temp += vec[k]*rotVec[j][k];
                                for (k = 0; k < num; k++)
                                        vec[k] -= temp*rotVec[j][k];
                        }
                        temp = 0.0;
                        for (j = 0; j < num; j++)
                                temp += vec[j]*vec[j];
                        temp = sqrt(temp);
                        for (j = 0; j < num; j++)
                                rotVec[i][j] = vec[j]/temp;
                }
                
                currentVec = rotVec;
        }
        
        int Dim() const {return proj;}
        
        ~RandomPlane()
        {
                del <double> (rotVec, num);
        }
};

class RandomBasis : public BasicDevs
{
        private:
                double **rotVec;
                
        protected:
                int num;
                double **currentVec;
                double **endVec;

        public:
                RandomBasis(int nin, unsigned long long iin) : BasicDevs(iin), num(nin)
                {
                        rotVec = matrix <double> (nin, nin);
                        RandRot();
                        endVec = currentVec + num;
                }
                
                void ChangeDim(const int nin)
                {
                        del <double> (rotVec, num);
                        num = nin;
                        rotVec = matrix <double> (nin, nin);
                        RandRot();
                        endVec = currentVec + num;
                }
                        
                void RandRot()
                {
                        double temp;
                        double vec[num];
                        int i, j, k;
                                
                        for (i = 0; i < num; i++)
                        {
                                for (j = 0; j < num; j++)
                                        vec[j] = Dev();
                                for (j = 0; j < i; j++)
                                {
                                        temp = 0.0;
                                        for (k = 0; k < num; k++)
                                                temp += vec[k]*rotVec[j][k];
                                        for (k = 0; k < num; k++)
                                                vec[k] -= temp*rotVec[j][k];
                                }
                                temp = 0.0;
                                for (j = 0; j < num; j++)
                                        temp += vec[j]*vec[j];
                                temp = sqrt(temp);
                                for (j = 0; j < num; j++)
                                        rotVec[i][j] = vec[j]/temp;
                        }
                        
                        currentVec = rotVec;
                }
                        
                double RanMult(double **cin)
                {
                        int i, j;
                        double temp = 0.0;
                        
                        for (i = 0; i < num; i++)
                        {
                                for (j = 0; j < num; j++)
                                {
                                        temp += (*currentVec)[i]*cin[i][j]*(*currentVec)[j];
                                }
                        }
                        return temp;
                }
                        
                void RanMult(const double in, double *out)
                {
                        int i;
                        double *ptr = *currentVec;
                        
                        for (i = 0; i < num; i++)
                        {
                                *out++ = in*(*ptr++);
                                        
                        }
                        return;
                }
                        
                void RanMult(double *in, const double w, double *out)
                {
                        int i;
                        double *ptr = *currentVec;
                        
                        for (i = 0; i < num; i++)
                        {
                                *out++ = *in++ + w*(*ptr++);
                        }
                        return;
                }
                        
                double Mag(double *a, double *a0)
                {
                        double temp = 0.0;
                        int i;
                        double *ptr = *currentVec;
                        
                        for (i = 0; i < num; i++)
                                temp += (*a++-*a0++)*(*ptr++);
                                
                        return temp;
                }
                        
                void Adjust(double *a, const double lim, const int iin)
                {
                        double temp = (lim - a[iin])/(*currentVec)[iin];
                        double *ptr = *currentVec;
                        
                        for (int i = 0; i < num; i++)
                                a[i] += temp*(*ptr++);
                        return;
                }
                        
                virtual void operator ++ (int)
                {
                        if (++currentVec == endVec)
                        {
                                RandRot();
                        }
                }
                
                virtual ~RandomBasis()
                {
                        del <double> (rotVec, num);
                }
};

class TransformRandomBasis : public RandomBasis, public Cholesky
{
        private:
                int num;
                
        public:
                TransformRandomBasis(double **vvar, int nin, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(vvar, nin), num(nin)
                {
                        double **ptr = currentVec;
                        for (int i = 0; i < num; i++)
                        {
                                ElMult(*ptr++);
                        }
                }
        
                void operator ++ (int)
                {
                        if (++currentVec == endVec)
                        {
                                RandRot();
                                double **ptr = currentVec;
                                for (int i = 0; i < num; i++)
                                {
                                        ElMult(*ptr++);
                                }
                        }
                }
};

class MultiNormDev : public RandomBasis, public Cholesky
{
        private:
                double fac;
                int num;
                
        public:
                MultiNormDev(int nin, double din, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(nin), fac(din), num(nin)
                {
                }
                
                MultiNormDev(double **vvar, const int nin, double din, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(vvar, nin), fac(din), num(nin)
                {
                }
                
                void MultiDev(double *pin, double *p0)
                {
                        int i;
                        double dist = 0.0;

                        if(Doub() < 0.33)
                                dist = ExpDev();
                        else
                        {
                                double temp;
                                for (i = 0; i < 2; i++)
                                {
                                        temp = Dev();
                                        dist += temp*temp;
                                }
                                dist = sqrt(dist/2.0);
                        }

                        ElMult(*currentVec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fac*dist*pin[i] + p0[i];
                        }
                        (*this)++;

                        return;
                }
                
                void MultiDev(double **cvar, double *pin, double *p0)
                {
                        int i;
                        double dist = 0.0;

                        if(Doub() < 0.33)
                                dist = ExpDev();
                        else
                        {
                                double temp;
                                for (i = 0; i < 2; i++)
                                {
                                        temp = Dev();
                                        dist += temp*temp;
                                }
                                dist = sqrt(dist/2.0);
                        }

                        EnterMat(cvar);
                        ElMult(*currentVec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fac*dist*pin[i] + p0[i];
                        }
                        (*this)++;

                        return;
                }
                
                void MultiDevGauss(double **cvar, double *pin, double *p0)
                {
                        int i;
                        double dist = 0.0;

                        double temp;
                        for (i = 0; i < 2; i++)
                        {
                                temp = Dev();
                                dist += temp*temp;
                        }
                        dist = sqrt(dist/2.0);
                        

                        EnterMat(cvar);
                        ElMult(*currentVec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fac*dist*pin[i] + p0[i];
                        }
                        (*this)++;

                        return;
                }
                
                void EllipseDev(double *pin, double *p0, double fin)
                {
                        int i;
                        double dist = 0.0;

                        dist = pow(Doub(), 1.0/num);
                        
                        ElMult(*currentVec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fin*dist*pin[i] + p0[i];
                        }
                        (*this)++;

                        return;
                }
                
                void EllipseDev(double **cvar, double *pin, double *p0, double fin)
                {
                        int i;
                        double dist = 0.0;

                        dist = pow(Doub(), 1.0/num);
                        
                        EnterMat(cvar);
                        ElMult(*currentVec, pin);
                        for (i = 0; i < num; i++)
                        {
                                pin[i] = fin*dist*pin[i] + p0[i];
                        }
                        (*this)++;

                        return;
                }
};

#endif