file simple/simple/random_tools.hpp

[No description available]

Classes

Name
classCholesky
classRan_old
classRan
classExponDev
classNormalDev
classBasicDevs
classMultiNormalDev
classAdvanceDevs
classRandomPlane
classRandomBasis
classTransformRandomBasis
classMultiNormDev

Functions

Name
template <typename T >
T *
matrix(const int xN)
template <typename T >
T **
matrix(const int xN, const int yN)
template <typename T >
T ***
matrix(const int xN, const int yN, const int zN)
template <typename T >
T *
matrix(const int xN, T in)
template <typename T >
T **
matrix(const int xN, const int yN, T in)
template <typename T >
T ***
matrix(const int xN, const int yN, const int zN, T in)
template <typename T >
T *
matrix(const int xN, T * in)
template <typename T >
T **
matrix(const int xN, const int yN, T ** in)
template <typename T >
T ***
matrix(const int xN, const int yN, const int zN, T *** in)
template <typename T >
void
del(T * temp)
template <typename T >
void
del(T ** temp, int xN)
template <typename T >
void
del(T *** temp, int xN, int yN)
template <class T >
const T
SQ(const T a)
template <class T >
const T
SQR(const T a)

Functions Documentation

function matrix

template <typename T >
T * matrix(
    const int xN
)

function matrix

template <typename T >
T ** matrix(
    const int xN,
    const int yN
)

function matrix

template <typename T >
T *** matrix(
    const int xN,
    const int yN,
    const int zN
)

function matrix

template <typename T >
T * matrix(
    const int xN,
    T in
)

function matrix

template <typename T >
T ** matrix(
    const int xN,
    const int yN,
    T in
)

function matrix

template <typename T >
T *** matrix(
    const int xN,
    const int yN,
    const int zN,
    T in
)

function matrix

template <typename T >
T * matrix(
    const int xN,
    T * in
)

function matrix

template <typename T >
T ** matrix(
    const int xN,
    const int yN,
    T ** in
)

function matrix

template <typename T >
T *** matrix(
    const int xN,
    const int yN,
    const int zN,
    T *** in
)

function del

template <typename T >
void del(
    T * temp
)

function del

template <typename T >
void del(
    T ** temp,
    int xN
)

function del

template <typename T >
void del(
    T *** temp,
    int xN,
    int yN
)

function SQ

template <class T >
inline const T SQ(
    const T a
)

function SQR

template <class T >
inline const T SQR(
    const T a
)

Source code

#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

Updated on 2024-07-18 at 13:53:32 +0000