file simple/simple/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;
//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