random_tools.hpp

[No description available]

Classes

	Name
class	Cholesky
class	Ran_old
class	Ran
class	ExponDev
class	NormalDev
class	BasicDevs
class	MultiNormalDev
class	AdvanceDevs
class	RandomPlane
class	RandomBasis
class	TransformRandomBasis
class	MultiNormDev

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;

            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;
        }

        bool EnterMat(const std::vector<std::vector<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)
        {
        }

        double MultiDevDist()
        {
            double dist = 0.0;

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

            return fac*dist;
        }

        double MultiDevPDF(double r, int dim)
        {
            return -(2.0 - dim)*std::log(r/fac) - r*r/fac/fac/2.0;
        }

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

            std::vector<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[0], 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;

            std::vector<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[0], 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;

            std::vector<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[0], 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;

            std::vector<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[0], 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 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), 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 - 2.0)*log(-Z);
    }

    void Mult2(double *ptrOut, double *ptr, double *ptr0, const double Z)
    {
        std::vector<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 Mult(double *ptrOut, double *ptr, double *ptr0, const double Z)
    {
        RandRot(0, proj);

        std::vector<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 HopBlow(double *ptrOut, double *ptrIn, double *ptr, double *ptr0)
    {
        double max = 0;

        RandRot(0, proj);

        std::vector<double> z(proj);
        for (int i = 0; i < proj; i++)
        {
            z[i] = 0.0;
            for (int j = 0; j < num; j++)
            {
                z[i] += (ptr[j] - ptr0[j])*currentVec[i][j];
            }
        }

        for (auto &&elem : z)
        {
            double pt = std::abs(elem);
            if (pt > max)
                max = pt;
        }

        double r = MultiDevDist();

        for (int j = 0; j < num; j++)
        {
            ptrOut[j] = ptrIn[j];
            for (int i = 0; i < proj; i++)
            {
                ptrOut[j] += r*max*currentVec[i][j]/3.0;
            }
        }
    }

    void RandRot(const int start = 0)
    {
        double temp;
        std::vector<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;
        std::vector<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;
            std::vector<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

file twalk/twalk/random_tools.hpp #

Classes #

Functions #

Functions Documentation #

function matrix #

function matrix #

function matrix #

function matrix #

function matrix #

function matrix #

function matrix #

function matrix #

function matrix #

function del #

function del #

function del #

function SQ #

function SQR #

Source code #

file twalk/twalk/random_tools.hpp

Classes

Functions

Functions Documentation

function matrix

function matrix

function matrix

function matrix

function matrix

function matrix

function matrix

function matrix

function matrix

function del

function del

function del

function SQ

function SQR

Source code