Commit 35b0aa2b authored by Joseph Pallamidessi's avatar Joseph Pallamidessi

Retabbing, and other minor aesthetic change

parent 541b0b9d
......@@ -7,74 +7,74 @@
//Functions for cma
long Caleatoire::alea_Start(long unsigned inseed)
{
long tmp;
int i;
this->flgstored = 0;
this->startseed = inseed;
if (inseed < 1)
inseed = 1;
this->aktseed = inseed;
for (i = 39; i >= 0; --i)
{
tmp = this->aktseed/127773;
this->aktseed = 16807 * (this->aktseed - tmp * 127773)
- 2836 * tmp;
if (this->aktseed < 0) this->aktseed += 2147483647;
if (i < 32)
this->rgalea[i] = this->aktseed;
}
this->aktalea = this->rgalea[0];
return inseed;
long tmp;
int i;
this->flgstored = 0;
this->startseed = inseed;
if (inseed < 1)
inseed = 1;
this->aktseed = inseed;
for (i = 39; i >= 0; --i)
{
tmp = this->aktseed/127773;
this->aktseed = 16807 * (this->aktseed - tmp * 127773)
- 2836 * tmp;
if (this->aktseed < 0) this->aktseed += 2147483647;
if (i < 32)
this->rgalea[i] = this->aktseed;
}
this->aktalea = this->rgalea[0];
return inseed;
}
long Caleatoire::alea_init(long unsigned inseed)
{
clock_t cloc = clock();
this->flgstored = 0;
if (inseed < 1) {
while ((long) (cloc - clock()) == 0)
; /* TODO: remove this for time critical applications? */
inseed = (long)abs((long)(100*time(NULL)+clock()));
}
return this->alea_Start(inseed);
clock_t cloc = clock();
this->flgstored = 0;
if (inseed < 1) {
while ((long) (cloc - clock()) == 0)
; /* TODO: remove this for time critical applications? */
inseed = (long)abs((long)(100*time(NULL)+clock()));
}
return this->alea_Start(inseed);
}
double Caleatoire::alea_Uniform()
{
long tmp;
tmp = this->aktseed/127773;
this->aktseed = 16807 * (this->aktseed - tmp * 127773)
- 2836 * tmp;
if (this->aktseed < 0)
this->aktseed += 2147483647;
tmp = this->aktalea / 67108865;
this->aktalea = this->rgalea[tmp];
this->rgalea[tmp] = this->aktseed;
return (double)(this->aktalea)/(2.147483647e9);
long tmp;
tmp = this->aktseed/127773;
this->aktseed = 16807 * (this->aktseed - tmp * 127773)
- 2836 * tmp;
if (this->aktseed < 0)
this->aktseed += 2147483647;
tmp = this->aktalea / 67108865;
this->aktalea = this->rgalea[tmp];
this->rgalea[tmp] = this->aktseed;
return (double)(this->aktalea)/(2.147483647e9);
}
double Caleatoire::alea_Gauss()
{
double x1, x2, rquad, fac;
if (this->flgstored)
{
this->flgstored = 0;
return this->hold;
}
do
{
x1 = 2.0 * this->alea_Uniform() - 1.0;
x2 = 2.0 * this->alea_Uniform() - 1.0;
rquad = x1*x1 + x2*x2;
} while(rquad >= 1 || rquad <= 0);
fac = sqrt(-2.0*log(rquad)/rquad);
this->flgstored = 1;
this->hold = fac * x1;
return fac * x2;
double x1, x2, rquad, fac;
if (this->flgstored)
{
this->flgstored = 0;
return this->hold;
}
do
{
x1 = 2.0 * this->alea_Uniform() - 1.0;
x2 = 2.0 * this->alea_Uniform() - 1.0;
rquad = x1*x1 + x2*x2;
} while(rquad >= 1 || rquad <= 0);
fac = sqrt(-2.0*log(rquad)/rquad);
this->flgstored = 1;
this->hold = fac * x1;
return fac * x2;
}
void CCmaes::Cmaes_init_param(int lambda, int mu){
......@@ -100,7 +100,7 @@ for (i=0; i<this->mu; ++i)
this->weights[i] /= s1;
if(this->mu < 1 || this->mu > this->lambda || (this->mu==this->lambda && this->weights[0]==this->weights[this->mu-1])){
printf("readpara_SetWeights(): invalid setting of mu or lambda\n");
exit(0);
exit(0);
}
/*supplement defaults*/
......@@ -126,314 +126,314 @@ this->updateCmode.modulo *= this->facupdateCmode;
while ((int) (cloc - clock()) == 0)
; /* TODO: remove this for time critical applications!? */
this->seed = (unsigned int)abs((long)(100*time(NULL)+clock()));
this->seed = (unsigned int)abs((long)(100*time(NULL)+clock()));
}
void CCmaes::TestMinStdDevs()
/* increases sigma */
{
int i;
if (this->rgDiffMinChange == NULL)
return;
else{
for (i = 0; i < this->dim; ++i)
while (this->sigma * sqrt(this->C[i][i]) < this->rgDiffMinChange[i])
this->sigma *= exp(0.05+this->cs/this->damps);
}
int i;
if (this->rgDiffMinChange == NULL)
return;
else{
for (i = 0; i < this->dim; ++i)
while (this->sigma * sqrt(this->C[i][i]) < this->rgDiffMinChange[i])
this->sigma *= exp(0.05+this->cs/this->damps);
}
} /* cmaes_TestMinStdDevs() */
int Check_Eigen(int taille, double **C, double *diag, double **Q)
{
/* compute Q diag Q^T and Q Q^T to check */
int i, j, k, res = 0;
double cc, dd;
for (i=0; i < taille; ++i)
for (j=0; j < taille; ++j) {
for (cc=0.,dd=0., k=0; k < taille; ++k) {
cc += diag[k] * Q[i][k] * Q[j][k];
dd += Q[i][k] * Q[j][k];
}
/* check here, is the normalization the right one? */
if (fabs(cc - C[i>j?i:j][i>j?j:i])/sqrt(C[i][i]*C[j][j]) > 1e-10 && fabs(cc - C[i>j?i:j][i>j?j:i]) > 3e-14)
{
printf("cmaes_t:Eigen(): imprecise result detected \n");
++res;
}
if (fabs(dd - (i==j)) > 1e-10) {
printf("cmaes_t:Eigen(): imprecise result detected (Q not orthog.)\n");
++res;
}
}
return res;
/* compute Q diag Q^T and Q Q^T to check */
int i, j, k, res = 0;
double cc, dd;
for (i=0; i < taille; ++i)
for (j=0; j < taille; ++j) {
for (cc=0.,dd=0., k=0; k < taille; ++k) {
cc += diag[k] * Q[i][k] * Q[j][k];
dd += Q[i][k] * Q[j][k];
}
/* check here, is the normalization the right one? */
if (fabs(cc - C[i>j?i:j][i>j?j:i])/sqrt(C[i][i]*C[j][j]) > 1e-10 && fabs(cc - C[i>j?i:j][i>j?j:i]) > 3e-14)
{
printf("cmaes_t:Eigen(): imprecise result detected \n");
++res;
}
if (fabs(dd - (i==j)) > 1e-10) {
printf("cmaes_t:Eigen(): imprecise result detected (Q not orthog.)\n");
++res;
}
}
return res;
}
double myhypot(double a, double b)
/* sqrt(a^2 + b^2) numerically stable. */
{
double r = 0;
if (fabs(a) > fabs(b)) {
r = b/a;
r = fabs(a)*sqrt(1+r*r);
} else if (b != 0) {
r = a/b;
r = fabs(b)*sqrt(1+r*r);
}
return r;
double r = 0;
if (fabs(a) > fabs(b)) {
r = b/a;
r = fabs(a)*sqrt(1+r*r);
} else if (b != 0) {
r = a/b;
r = fabs(b)*sqrt(1+r*r);
}
return r;
}
void Householder2(int n, double **V, double *d, double *e) {
int i,j,k;
int i,j,k;
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
}
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
}
/* Householder reduction to tridiagonal form */
/* Householder reduction to tridiagonal form */
for (i = n-1; i > 0; i--) {
for (i = n-1; i > 0; i--) {
/* Scale to avoid under/overflow */
/* Scale to avoid under/overflow */
double scale = 0.0;
double h = 0.0;
for (k = 0; k < i; k++) {
scale = scale + fabs(d[k]);
}
if (scale == 0.0) {
e[i] = d[i-1];
for (j = 0; j < i; j++) {
d[j] = V[i-1][j];
V[i][j] = 0.0;
V[j][i] = 0.0;
}
} else {
double scale = 0.0;
double h = 0.0;
for (k = 0; k < i; k++) {
scale = scale + fabs(d[k]);
}
if (scale == 0.0) {
e[i] = d[i-1];
for (j = 0; j < i; j++) {
d[j] = V[i-1][j];
V[i][j] = 0.0;
V[j][i] = 0.0;
}
} else {
/* Generate Householder vector */
/* Generate Householder vector */
double f, g, hh;
double f, g, hh;
for (k = 0; k < i; k++) {
d[k] /= scale;
h += d[k] * d[k];
}
f = d[i-1];
g = sqrt(h);
if (f > 0) {
g = -g;
}
e[i] = scale * g;
h = h - f * g;
d[i-1] = f - g;
for (j = 0; j < i; j++) {
e[j] = 0.0;
}
for (k = 0; k < i; k++) {
d[k] /= scale;
h += d[k] * d[k];
}
f = d[i-1];
g = sqrt(h);
if (f > 0) {
g = -g;
}
e[i] = scale * g;
h = h - f * g;
d[i-1] = f - g;
for (j = 0; j < i; j++) {
e[j] = 0.0;
}
/* Apply similarity transformation to remaining columns */
/* Apply similarity transformation to remaining columns */
for (j = 0; j < i; j++) {
f = d[j];
V[j][i] = f;
g = e[j] + V[j][j] * f;
for (k = j+1; k <= i-1; k++) {
g += V[k][j] * d[k];
e[k] += V[k][j] * f;
}
e[j] = g;
}
f = 0.0;
for (j = 0; j < i; j++) {
e[j] /= h;
f += e[j] * d[j];
}
hh = f / (h + h);
for (j = 0; j < i; j++) {
e[j] -= hh * d[j];
}
for (j = 0; j < i; j++) {
f = d[j];
g = e[j];
for (k = j; k <= i-1; k++) {
V[k][j] -= (f * e[k] + g * d[k]);
}
d[j] = V[i-1][j];
V[i][j] = 0.0;
}
}
d[i] = h;
}
for (j = 0; j < i; j++) {
f = d[j];
V[j][i] = f;
g = e[j] + V[j][j] * f;
for (k = j+1; k <= i-1; k++) {
g += V[k][j] * d[k];
e[k] += V[k][j] * f;
}
e[j] = g;
}
f = 0.0;
for (j = 0; j < i; j++) {
e[j] /= h;
f += e[j] * d[j];
}
hh = f / (h + h);
for (j = 0; j < i; j++) {
e[j] -= hh * d[j];
}
for (j = 0; j < i; j++) {
f = d[j];
g = e[j];
for (k = j; k <= i-1; k++) {
V[k][j] -= (f * e[k] + g * d[k]);
}
d[j] = V[i-1][j];
V[i][j] = 0.0;
}
}
d[i] = h;
}
/* Accumulate transformations */
/* Accumulate transformations */
for (i = 0; i < n-1; i++) {
double h;
V[n-1][i] = V[i][i];
V[i][i] = 1.0;
h = d[i+1];
if (h != 0.0) {
for (k = 0; k <= i; k++) {
d[k] = V[k][i+1] / h;
}
for (j = 0; j <= i; j++) {
double g = 0.0;
for (k = 0; k <= i; k++) {
g += V[k][i+1] * V[k][j];
}
for (k = 0; k <= i; k++) {
V[k][j] -= g * d[k];
}
}
}
for (k = 0; k <= i; k++) {
V[k][i+1] = 0.0;
}
}
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
V[n-1][j] = 0.0;
}
V[n-1][n-1] = 1.0;
e[0] = 0.0;
for (i = 0; i < n-1; i++) {
double h;
V[n-1][i] = V[i][i];
V[i][i] = 1.0;
h = d[i+1];
if (h != 0.0) {
for (k = 0; k <= i; k++) {
d[k] = V[k][i+1] / h;
}
for (j = 0; j <= i; j++) {
double g = 0.0;
for (k = 0; k <= i; k++) {
g += V[k][i+1] * V[k][j];
}
for (k = 0; k <= i; k++) {
V[k][j] -= g * d[k];
}
}
}
for (k = 0; k <= i; k++) {
V[k][i+1] = 0.0;
}
}
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
V[n-1][j] = 0.0;
}
V[n-1][n-1] = 1.0;
e[0] = 0.0;
} /* Housholder() */
void QLalgo2 (int n, double *d, double *e, double **V) {
int i, k, l, m;
double f = 0.0;
double tst1 = 0.0;
double eps = 2.22e-16; /* Math.pow(2.0,-52.0); == 2.22e-16 */
int i, k, l, m;
double f = 0.0;
double tst1 = 0.0;
double eps = 2.22e-16; /* Math.pow(2.0,-52.0); == 2.22e-16 */
/* shift input e */
for (i = 1; i < n; i++) {
e[i-1] = e[i];
}
e[n-1] = 0.0; /* never changed again */
/* shift input e */
for (i = 1; i < n; i++) {
e[i-1] = e[i];
}
e[n-1] = 0.0; /* never changed again */
for (l = 0; l < n; l++) {
for (l = 0; l < n; l++) {
/* Find small subdiagonal element */
/* Find small subdiagonal element */
if (tst1 < fabs(d[l]) + fabs(e[l]))
tst1 = fabs(d[l]) + fabs(e[l]);
m = l;
while (m < n) {
if (fabs(e[m]) <= eps*tst1) {
/* if (fabs(e[m]) + fabs(d[m]+d[m+1]) == fabs(d[m]+d[m+1])) { */
break;
}
m++;
}
if (tst1 < fabs(d[l]) + fabs(e[l]))
tst1 = fabs(d[l]) + fabs(e[l]);
m = l;
while (m < n) {
if (fabs(e[m]) <= eps*tst1) {
/* if (fabs(e[m]) + fabs(d[m]+d[m+1]) == fabs(d[m]+d[m+1])) { */
break;
}
m++;
}
/* If m == l, d[l] is an eigenvalue, */
/* otherwise, iterate. */
/* If m == l, d[l] is an eigenvalue, */
/* otherwise, iterate. */
if (m > l) {
int iter = 0;
do { /* while (fabs(e[l]) > eps*tst1); */
double dl1, h;
double g = d[l];
double p = (d[l+1] - g) / (2.0 * e[l]);
double r = myhypot(p, 1.);
iter = iter + 1; /* Could check iteration count here */
if (m > l) {
int iter = 0;
do { /* while (fabs(e[l]) > eps*tst1); */
double dl1, h;
double g = d[l];
double p = (d[l+1] - g) / (2.0 * e[l]);
double r = myhypot(p, 1.);
iter = iter + 1; /* Could check iteration count here */
/* Compute implicit shift */
/* Compute implicit shift */
if (p < 0) {
r = -r;
}
d[l] = e[l] / (p + r);
d[l+1] = e[l] * (p + r);
dl1 = d[l+1];
h = g - d[l];
for (i = l+2; i < n; i++) {
d[i] -= h;
}
f = f + h;
if (p < 0) {
r = -r;
}
d[l] = e[l] / (p + r);
d[l+1] = e[l] * (p + r);
dl1 = d[l+1];
h = g - d[l];
for (i = l+2; i < n; i++) {
d[i] -= h;
}
f = f + h;
/* Implicit QL transformation. */
/* Implicit QL transformation. */
p = d[m];
{
double c = 1.0;
double c2 = c;
double c3 = c;
double el1 = e[l+1];
double s = 0.0;
double s2 = 0.0;
for (i = m-1; i >= l; i--) {
c3 = c2;
c2 = c;
s2 = s;
g = c * e[i];
h = c * p;
r = myhypot(p, e[i]);
e[i+1] = s * r;
s = e[i] / r;
c = p / r;
p = c * d[i] - s * g;
d[i+1] = h + s * (c * g + s * d[i]);
p = d[m];
{
double c = 1.0;
double c2 = c;
double c3 = c;
double el1 = e[l+1];
double s = 0.0;
double s2 = 0.0;
for (i = m-1; i >= l; i--) {
c3 = c2;
c2 = c;
s2 = s;
g = c * e[i];
h = c * p;
r = myhypot(p, e[i]);
e[i+1] = s * r;
s = e[i] / r;
c = p / r;
p = c * d[i] - s * g;
d[i+1] = h + s * (c * g + s * d[i]);
/* Accumulate transformation. */
/* Accumulate transformation. */
for (k = 0; k < n; k++) {
h = V[k][i+1];
V[k][i+1] = s * V[k][i] + c * h;
V[k][i] = c * V[k][i] - s * h;
}
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d[l] = c * p;
}
for (k = 0; k < n; k++) {
h = V[k][i+1];
V[k][i+1] = s * V[k][i] + c * h;
V[k][i] = c * V[k][i] - s * h;
}
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d[l] = c * p;
}
/* Check for convergence. */
/* Check for convergence. */
} while (fabs(e[l]) > eps*tst1);
}
d[l] = d[l] + f;
e[l] = 0.0;
}
} while (fabs(e[l]) > eps*tst1);
}
d[l] = d[l] + f;
e[l] = 0.0;
}
/* Sort eigenvalues and corresponding vectors. */
int j;
double p;
for (i = 0; i < n-1; i++) {
k = i;
p = d[i];
for (j = i+1; j < n; j++) {
if (d[j] < p) {
k = j;
p = d[j];
}
}
if (k != i) {
d[k] = d[i];
d[i] = p;
for (j = 0; j < n; j++) {
p = V[j][i];
V[j][i] = V[j][k];
V[j][k] = p;
}
}
}
int j;
double p;
for (i = 0; i < n-1; i++) {
k = i;
p = d[i];
for (j = i+1; j < n; j++) {
if (d[j] < p) {
k = j;
p = d[j];
}
}
if (k != i) {
d[k] = d[i];
d[i] = p;
for (j = 0; j < n; j++) {
p = V[j][i];
V[j][i] = V[j][k];
V[j][k] = p;
}
}
}
} /* QLalgo2 */
void Eigen( int taille, double **C, double *diag, double **Q, double *rgtmp)
{
int i, j;
if (rgtmp == NULL) /* was OK in former versions */
printf("cmaes_t:Eigen(): input parameter double *rgtmp must be non-NULL\n");
/* copy C to Q */
if (C != Q) {
for (i=0; i < taille; ++i)
for (j = 0; j <= i; ++j)
Q[i][j] = Q[j][i] = C[i][j];