From 00d1d5d6a770959776bd27d1dddd75864de1cc0b Mon Sep 17 00:00:00 2001 From: Guangyang Wen Date: Sun, 19 Nov 2017 14:03:54 +0100 Subject: [PATCH] refractor --- CMakeLists.txt | 1 + include/cma-es/cmaes.h | 606 ++++++++++++++++++++---------------- include/cma-es/parameters.h | 125 ++++---- include/cma-es/random.h | 16 +- include/cma-es/utils.h | 34 +- include/lexi.h | 60 +++- src/conversion.cpp | 4 +- test/example1.cpp | 40 ++- test/example2.cpp | 222 ++++++------- test/testlexi.cpp | 42 ++- 10 files changed, 650 insertions(+), 500 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 98a4009..b75a9f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ cmake_minimum_required (VERSION 3.1) project(color) set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wconversion -Wno-unused-private-field -Wno-unused-variable -Wno-unused-parameter") set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") diff --git a/include/cma-es/cmaes.h b/include/cma-es/cmaes.h index 4f8e8c9..95ae82f 100644 --- a/include/cma-es/cmaes.h +++ b/include/cma-es/cmaes.h @@ -109,13 +109,28 @@ #include #include #include +#include + +/** + * @class Individual + * + * @field x vector of parameters of native scalar type (float, double, etc.) + * @field value linearly ordered object of the the object function value, which + * supports all comparison operators, additive(+/-) group operations and + * initialization through 0 or std::numeric_limit::max() + */ +template class Individual { +public: + std::vector x; + Ordered value; +}; /** * @class CMAES * Evolution Strategies with Covariance Matrix Adaptation. The public interface * of the optimization algorithm. */ -template class CMAES { +template class CMAES { public: /** * Keys for get(). @@ -190,65 +205,67 @@ template class CMAES { //! Implementation version. std::string version; //!< Random number generator. - Random rand; + Random rand; //!< CMA-ES parameters. - Parameters params; + Parameters params; //! Step size. - T sigma; + Scalar sigma; //! Mean x vector, "parent". - T *xmean; + Scalar *xmean; //! Best sample ever. - T *xBestEver; + Individual xBestEver; + Scalar evalsBestEver; //! x-vectors, lambda offspring. - T **population; + Individual *population; //! Sorting index of sample population. - int *index; - //! History of function values. - T *funcValueHistory; + size_t *index; + //! History of last historySize function values. + Ordered *funcValueHistory; + size_t historySize; - T chiN; + Scalar chiN; //! Lower triangular matrix: i>=j for C[i][j]. - T **C; + Scalar **C; //! Matrix with normalize eigenvectors in columns. - T **B; + Scalar **B; //! Axis lengths. - T *rgD; + Scalar *rgD; //! Anisotropic evolution path (for covariance). - T *pc; + Scalar *pc; //! Isotropic evolution path (for step length). - T *ps; + Scalar *ps; //! Last mean. - T *xold; + Scalar *xold; //! Output vector. - T *output; + // Scalar *output; //! B*D*z. - T *BDz; + Scalar *BDz; //! Temporary (random) vector used in different places. - T *tempRandom; + Scalar *tempRandom; //! Objective function values of the population. - T *functionValues; + Ordered *functionValues; //!< Public objective function value array returned by init(). - T *publicFitness; + Ordered *publicFitness; //! Generation number. - T gen; + size_t gen; //! Algorithm state. enum { INITIALIZED, SAMPLED, UPDATED } state; // repeatedly used for output - T maxdiagC; - T mindiagC; - T maxEW; - T minEW; + Scalar maxdiagC; + Scalar mindiagC; + Scalar maxEW; + Scalar minEW; bool eigensysIsUptodate; bool doCheckEigen; //!< control via signals.par - T genOfEigensysUpdate; + Scalar genOfEigensysUpdate; - T dMaxSignifKond; - T dLastMinEWgroesserNull; + Scalar dMaxSignifKond; + Scalar dLastMinEWgroesserNull; bool isResumeDone; @@ -272,13 +289,13 @@ template class CMAES { * @param diag (output) N eigenvalues. * @param Q (output) Columns are normalized eigenvectors. */ - void eigen(T *diag, T **Q, T *rgtmp) { + void eigen(Scalar *diag, Scalar **Q, Scalar *rgtmp) { assert(rgtmp && "eigen(): input parameter rgtmp must be non-NULL"); if (C != Q) // copy C to Q { - for (int i = 0; i < params.N; ++i) - for (int j = 0; j <= i; ++j) + for (size_t i = 0; i < params.N; ++i) + for (size_t j = 0; j <= i; ++j) Q[i][j] = Q[j][i] = C[i][j]; } @@ -291,12 +308,12 @@ template class CMAES { * operations writes to error file. * @return number of detected inaccuracies */ - int checkEigen(T *diag, T **Q) { - // compute Q diag Q^T and Q Q^T to check + int checkEigen(Scalar *diag, Scalar **Q) { + // compute Q diag Q^Scalar and Q Q^Scalar to check int res = 0; for (int i = 0; i < params.N; ++i) for (int j = 0; j < params.N; ++j) { - T cc = 0., dd = 0.; + Scalar cc = 0., dd = 0.; for (int k = 0; k < params.N; ++k) { cc += diag[k] * Q[i][k] * Q[j][k]; dd += Q[i][k] * Q[j][k]; @@ -304,9 +321,9 @@ template class CMAES { // check here, is the normalization the right one? const bool cond1 = fabs(cc - C[i > j ? i : j][i > j ? j : i]) / sqrt(C[i][i] * C[j][j]) > - T(1e-10); + Scalar(1e-10); const bool cond2 = - fabs(cc - C[i > j ? i : j][i > j ? j : i]) > T(3e-14); + fabs(cc - C[i > j ? i : j][i > j ? j : i]) > Scalar(3e-14); if (cond1 && cond2) { std::stringstream s; s << i << " " << j << ": " << cc << " " @@ -317,7 +334,7 @@ template class CMAES { << std::endl; ++res; } - if (std::fabs(dd - (i == j)) > T(1e-10)) { + if (std::fabs(dd - (i == j)) > Scalar(1e-10)) { std::stringstream s; s << i << " " << j << " " << dd; if (params.logWarnings) @@ -339,26 +356,26 @@ template class CMAES { * @param V input: matrix output of Householder. output: basis of * eigenvectors, according to d */ - void ql(T *d, T *e, T **V) { + void ql(Scalar *d, Scalar *e, Scalar **V) { const int n = params.N; - T f(0); - T tst1(0); - const T eps(2.22e-16); // 2.0^-52.0 = 2.22e-16 + Scalar f(0); + Scalar tst1(0); + const Scalar eps(2.22e-16); // 2.0^-52.0 = 2.22e-16 // shift input e - T *ep1 = e; - for (T *ep2 = e + 1, *const end = e + n; ep2 != end; ep1++, ep2++) + Scalar *ep1 = e; + for (Scalar *ep2 = e + 1, *const end = e + n; ep2 != end; ep1++, ep2++) *ep1 = *ep2; - *ep1 = T(0); // never changed again + *ep1 = Scalar(0); // never changed again for (int l = 0; l < n; l++) { // find small subdiagonal element - T &el = e[l]; - T &dl = d[l]; - const T smallSDElement = std::fabs(dl) + std::fabs(el); + Scalar &el = e[l]; + Scalar &dl = d[l]; + const Scalar smallSDElement = std::fabs(dl) + std::fabs(el); if (tst1 < smallSDElement) tst1 = smallSDElement; - const T epsTst1 = eps * tst1; + const Scalar epsTst1 = eps * tst1; int m = l; while (m < n) { if (std::fabs(e[m]) <= epsTst1) @@ -369,18 +386,18 @@ template class CMAES { // if m == l, d[l] is an eigenvalue, otherwise, iterate. if (m > l) { do { - T h, g = dl; - T &dl1r = d[l + 1]; - T p = (dl1r - g) / (T(2) * el); - T r = myhypot(p, T(1)); + Scalar h, g = dl; + Scalar &dl1r = d[l + 1]; + Scalar p = (dl1r - g) / (Scalar(2) * el); + Scalar r = myhypot(p, Scalar(1)); // compute implicit shift if (p < 0) r = -r; - const T pr = p + r; + const Scalar pr = p + r; dl = el / pr; h = g - dl; - const T dl1 = el * pr; + const Scalar dl1 = el * pr; dl1r = dl1; for (int i = l + 2; i < n; i++) d[i] -= h; @@ -388,32 +405,32 @@ template class CMAES { // implicit QL transformation. p = d[m]; - T c(1); - T c2(1); - T c3(1); - const T el1 = e[l + 1]; - T s(0); - T s2(0); + Scalar c(1); + Scalar c2(1); + Scalar c3(1); + const Scalar el1 = e[l + 1]; + Scalar s(0); + Scalar s2(0); for (int i = m - 1; i >= l; i--) { c3 = c2; c2 = c; s2 = s; - const T &ei = e[i]; + const Scalar &ei = e[i]; g = c * ei; h = c * p; r = myhypot(p, ei); e[i + 1] = s * r; s = ei / r; c = p / r; - const T &di = d[i]; + const Scalar &di = d[i]; p = c * di - s * g; d[i + 1] = h + s * (c * g + s * di); // accumulate transformation. for (int k = 0; k < n; k++) { - T &Vki1 = V[k][i + 1]; + Scalar &Vki1 = V[k][i + 1]; h = Vki1; - T &Vki = V[k][i]; + Scalar &Vki = V[k][i]; Vki1 = s * Vki + c * h; Vki *= c; Vki -= s * h; @@ -437,7 +454,7 @@ template class CMAES { * @param d output: diagonal * @param e output: [0..n-1], off diagonal (elements 1..n-1) */ - void householder(T **V, T *d, T *e) { + void householder(Scalar **V, Scalar *d, Scalar *e) { const int n = params.N; for (int j = 0; j < n; j++) { @@ -448,9 +465,9 @@ template class CMAES { for (int i = n - 1; i > 0; i--) { // scale to avoid under/overflow - T scale = 0.0; - T h = 0.0; - for (T *pd = d, *const dend = d + i; pd != dend; pd++) { + Scalar scale = 0.0; + Scalar h = 0.0; + for (Scalar *pd = d, *const dend = d + i; pd != dend; pd++) { scale += std::fabs(*pd); } if (scale == 0.0) { @@ -462,26 +479,26 @@ template class CMAES { } } else { // generate Householder vector - for (T *pd = d, *const dend = d + i; pd != dend; pd++) { + for (Scalar *pd = d, *const dend = d + i; pd != dend; pd++) { *pd /= scale; h += *pd * *pd; } - T &dim1 = d[i - 1]; - T f = dim1; - T g = f > 0 ? -std::sqrt(h) : std::sqrt(h); + Scalar &dim1 = d[i - 1]; + Scalar f = dim1; + Scalar g = f > 0 ? -std::sqrt(h) : std::sqrt(h); e[i] = scale * g; h = h - f * g; dim1 = f - g; - memset((void *)e, 0, (size_t)i * sizeof(T)); + memset((void *)e, 0, (size_t)i * sizeof(Scalar)); // apply similarity transformation to remaining columns for (int j = 0; j < i; j++) { f = d[j]; V[j][i] = f; - T &ej = e[j]; + Scalar &ej = e[j]; g = ej + V[j][j] * f; for (int k = j + 1; k <= i - 1; k++) { - T &Vkj = V[k][j]; + Scalar &Vkj = V[k][j]; g += Vkj * d[k]; e[k] += Vkj * f; } @@ -489,16 +506,16 @@ template class CMAES { } f = 0.0; for (int j = 0; j < i; j++) { - T &ej = e[j]; + Scalar &ej = e[j]; ej /= h; f += ej * d[j]; } - T hh = f / (h + h); + Scalar hh = f / (h + h); for (int j = 0; j < i; j++) { e[j] -= hh * d[j]; } for (int j = 0; j < i; j++) { - T &dj = d[j]; + Scalar &dj = d[j]; f = dj; g = e[j]; for (int k = j; k <= i - 1; k++) { @@ -514,8 +531,8 @@ template class CMAES { // accumulate transformations const int nm1 = n - 1; for (int i = 0; i < nm1; i++) { - T h; - T &Vii = V[i][i]; + Scalar h; + Scalar &Vii = V[i][i]; V[n - 1][i] = Vii; Vii = 1.0; h = d[i + 1]; @@ -524,9 +541,9 @@ template class CMAES { d[k] = V[k][i + 1] / h; } for (int j = 0; j <= i; j++) { - T g = 0.0; + Scalar g = 0.0; for (int k = 0; k <= i; k++) { - T *Vk = V[k]; + Scalar *Vk = V[k]; g += Vk[i + 1] * Vk[j]; } for (int k = 0; k <= i; k++) { @@ -539,7 +556,7 @@ template class CMAES { } } for (int j = 0; j < n; j++) { - T &Vnm1j = V[n - 1][j]; + Scalar &Vnm1j = V[n - 1][j]; d[j] = Vnm1j; Vnm1j = 0.0; } @@ -550,8 +567,9 @@ template class CMAES { /** * Dirty index sort. */ - void sortIndex(const T *rgFunVal, int *iindex, int n) { - int i, j; + template + void sortIndex(const LessThanComparable *rgFunVal, size_t *iindex, size_t n) { + size_t i, j; for (i = 1, iindex[0] = 0; i < n; ++i) { for (j = i; j > 0; --j) { if (rgFunVal[iindex[j - 1]] < rgFunVal[i]) @@ -566,27 +584,29 @@ template class CMAES { const int N = params.N; bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; - if (params.ccov != T(0)) { + if (params.ccov != Scalar(0)) { // definitions for speeding up inner-most loop - const T mucovinv = T(1) / params.mucov; - const T commonFactor = params.ccov * (diag ? (N + T(1.5)) / T(3) : T(1)); - const T ccov1 = std::min(commonFactor * mucovinv, T(1)); - const T ccovmu = std::min(commonFactor * (T(1) - mucovinv), T(1) - ccov1); - const T sigmasquare = sigma * sigma; - const T onemccov1ccovmu = T(1) - ccov1 - ccovmu; - const T longFactor = - (T(1) - hsig) * params.ccumcov * (T(2) - params.ccumcov); + const Scalar mucovinv = Scalar(1) / params.mucov; + const Scalar commonFactor = + params.ccov * (diag ? (N + Scalar(1.5)) / Scalar(3) : Scalar(1)); + const Scalar ccov1 = std::min(commonFactor * mucovinv, Scalar(1)); + const Scalar ccovmu = + std::min(commonFactor * (Scalar(1) - mucovinv), Scalar(1) - ccov1); + const Scalar sigmasquare = sigma * sigma; + const Scalar onemccov1ccovmu = Scalar(1) - ccov1 - ccovmu; + const Scalar longFactor = + (Scalar(1) - hsig) * params.ccumcov * (Scalar(2) - params.ccumcov); eigensysIsUptodate = false; // update covariance matrix - for (int i = 0; i < N; ++i) - for (int j = diag ? i : 0; j <= i; ++j) { - T &Cij = C[i][j]; + for (size_t i = 0; i < N; ++i) + for (size_t j = diag ? i : 0; j <= i; ++j) { + Scalar &Cij = C[i][j]; Cij = onemccov1ccovmu * Cij + ccov1 * (pc[i] * pc[j] + longFactor * Cij); for (int k = 0; k < params.mu; ++k) { // additional rank mu update - const T *rgrgxindexk = population[index[k]]; + const std::vector &rgrgxindexk = population[index[k]].x; Cij += ccovmu * params.weights[k] * (rgrgxindexk[i] - xold[i]) * (rgrgxindexk[j] - xold[j]) / sigmasquare; } @@ -594,7 +614,7 @@ template class CMAES { // update maximal and minimal diagonal value maxdiagC = mindiagC = C[0][0]; for (int i = 1; i < N; ++i) { - const T &Cii = C[i][i]; + const Scalar &Cii = C[i][i]; if (maxdiagC < Cii) maxdiagC = Cii; else if (mindiagC > Cii) @@ -613,7 +633,8 @@ template class CMAES { for (int i = 0; i < params.N; ++i) while (this->sigma * std::sqrt(this->C[i][i]) < this->params.rgDiffMinChange[i]) - this->sigma *= std::exp(T(0.05) + this->params.cs / this->params.damps); + this->sigma *= + std::exp(Scalar(0.05) + this->params.cs / this->params.damps); } /** @@ -621,11 +642,11 @@ template class CMAES { * @param x Search space vector. * @param eps Mutation factor. */ - void addMutation(T *x, T eps = 1.0) { + void addMutation(std::vector &x, Scalar eps = 1.0) { for (int i = 0; i < params.N; ++i) tempRandom[i] = rgD[i] * rand.gauss(); for (int i = 0; i < params.N; ++i) { - T sum = 0.0; + Scalar sum = 0.0; for (int j = 0; j < params.N; ++j) sum += B[i][j] * tempRandom[j]; x[i] = xmean[i] + eps * sigma * sum; @@ -678,52 +699,52 @@ template class CMAES { file << "function evaluations " << (long)countevals << std::endl; file << "elapsed (CPU) time [s] " << std::setprecision(2) << eigenTimings.totaltotaltime << std::endl; - file << "function value f(x)=" << population[index[0]][params.N] - << std::endl; + file << "function value f(x)=" << population[index[0]].value << std::endl; file << "maximal standard deviation " << sigma * std::sqrt(maxdiagC) << std::endl; file << "minimal standard deviation " << sigma * std::sqrt(mindiagC) << std::endl; file << "sigma " << sigma << std::endl; file << "axisratio " - << (maxElement(rgD, params.N) / minElement(rgD, params.N)) + << (maxElement(rgD, (size_t)params.N) / + minElement(rgD, (size_t)params.N)) << std::endl; - file << "xbestever found after " << std::setprecision(0) - << xBestEver[params.N + 1] << "evaluations, function value " - << xBestEver[params.N] << std::endl; - for (int i = 0; i < params.N; ++i) - file << " " << std::setw(12) << xBestEver[i] + file << "xbestever found after " << std::setprecision(0) << evalsBestEver + << "evaluations, function value " << xBestEver.value << std::endl; + for (size_t i = 0; i < params.N; ++i) + file << " " << std::setw(12) << xBestEver.x[i] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "xbest (of last generation, function value " - << population[index[0]][params.N] << ")" << std::endl; - for (int i = 0; i < params.N; ++i) - file << " " << std::setw(12) << population[index[0]][i] + << population[index[0]].value << ")" << std::endl; + for (size_t i = 0; i < params.N; ++i) + file << " " << std::setw(12) << population[index[0]].x[i] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "xmean" << std::endl; - for (int i = 0; i < params.N; ++i) + for (size_t i = 0; i < params.N; ++i) file << " " << std::setw(12) << xmean[i] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "Standard deviation of coordinate axes (sigma*sqrt(diag(C)))" << std::endl; - for (int i = 0; i < params.N; ++i) + for (size_t i = 0; i < params.N; ++i) file << " " << std::setw(12) << sigma * std::sqrt(C[i][i]) << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "Main axis lengths of mutation ellipsoid (sigma*diag(D))" << std::endl; - for (int i = 0; i < params.N; ++i) + for (size_t i = 0; i < params.N; ++i) tempRandom[i] = rgD[i]; std::sort(tempRandom, tempRandom + params.N); - for (int i = 0; i < params.N; ++i) - file << " " << std::setw(12) << sigma * tempRandom[params.N - 1 - i] + for (size_t i = 0; i < params.N; ++i) + file << " " << std::setw(12) + << sigma * tempRandom[(size_t)params.N - 1 - i] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "Longest axis (b_i where d_ii=max(diag(D))" << std::endl; - int k = maxIndex(rgD, params.N); - for (int i = 0; i < params.N; ++i) + size_t k = maxIndex(rgD, (size_t)params.N); + for (size_t i = 0; i < params.N; ++i) file << " " << std::setw(12) << B[i][k] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << "Shortest axis (b_i where d_ii=max(diag(D))" << std::endl; - k = minIndex(rgD, params.N); - for (int i = 0; i < params.N; ++i) + k = minIndex(rgD, (size_t)params.N); + for (size_t i = 0; i < params.N; ++i) file << " " << std::setw(12) << B[i][k] << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); file << std::endl; @@ -752,7 +773,7 @@ template class CMAES { file << std::endl; } if (key & WKFBestEver) { - file << xBestEver[params.N] << std::endl; + file << xBestEver.value << std::endl; } if (key & WKCGeneration) { file << gen << std::endl; @@ -764,8 +785,8 @@ template class CMAES { file << params.lambda << std::endl; } if (key & WKB) { - int *iindex = new int[params.N]; - sortIndex(rgD, iindex, params.N); + size_t *iindex = new size_t[params.N]; + sortIndex(rgD, iindex, (size_t)params.N); for (int i = 0; i < params.N; ++i) for (int j = 0; j < params.N; ++j) { file << B[j][iindex[params.N - 1 - i]]; @@ -779,8 +800,8 @@ template class CMAES { file << std::endl; } if (key & WKXBest) { - for (int i = 0; i < params.N; ++i) - file << (i == 0 ? "" : "\t") << population[index[0]][i]; + for (size_t i = 0; i < params.N; ++i) + file << (i == 0 ? "" : "\t") << population[index[0]].x[i]; file << std::endl; } if (key & WKClock) { @@ -795,10 +816,10 @@ template class CMAES { } public: - T countevals; //!< objective function evaluations + Scalar countevals; //!< objective function evaluations Timing eigenTimings; - CMAES() : version("1.0alpha") {} + CMAES() : version("1.1alpha") {} /** * Releases the dynamically allocated memory, including that of the return @@ -809,24 +830,24 @@ template class CMAES { delete[] ps; delete[] tempRandom; delete[] BDz; - delete[]-- xmean; - delete[]-- xold; - delete[]-- xBestEver; - delete[]-- output; + delete[] xmean; + delete[] xold; + // delete[]-- xBestEver; + // delete[]-- output; delete[] rgD; for (int i = 0; i < params.N; ++i) { delete[] C[i]; delete[] B[i]; } - for (int i = 0; i < params.lambda; ++i) - delete[]-- population[i]; + // for (int i = 0; i < params.lambda; ++i) + // delete[]-- population[i]; delete[] population; delete[] C; delete[] B; delete[] index; delete[] publicFitness; delete[] functionValues; - delete[]-- funcValueHistory; + delete[] funcValueHistory; } /** @@ -836,101 +857,110 @@ template class CMAES { * pass them to updateDistribution(). Not that after the desctructor * was called, the array is deleted. */ - T *init(const Parameters ¶meters) { + Ordered *init(const Parameters ¶meters) { params = parameters; stopMessage = ""; - T trace(0); + Scalar trace(0); for (int i = 0; i < params.N; ++i) trace += params.rgInitialStds[i] * params.rgInitialStds[i]; sigma = std::sqrt(trace / params.N); - chiN = std::sqrt((T)params.N) * (T(1) - T(1) / (T(4) * params.N) + - T(1) / (T(21) * params.N * params.N)); + chiN = std::sqrt((Scalar)params.N) * + (Scalar(1) - Scalar(1) / (Scalar(4) * params.N) + + Scalar(1) / (Scalar(21) * params.N * params.N)); eigensysIsUptodate = true; doCheckEigen = false; genOfEigensysUpdate = 0; isResumeDone = false; - T dtest; - for (dtest = T(1); dtest && dtest < T(1.1) * dtest; dtest *= T(2)) - if (dtest == dtest + T(1)) + Scalar dtest; + for (dtest = Scalar(1); dtest != 0 && dtest < Scalar(1.1) * dtest; + dtest *= Scalar(2)) + if (dtest == dtest + Scalar(1)) break; - dMaxSignifKond = dtest / T(1000); // not sure whether this is really save, - // 100 does not work well enough + dMaxSignifKond = + dtest / Scalar(1000); // not sure whether this is really save, + // 100 does not work well enough gen = 0; countevals = 0; state = INITIALIZED; - dLastMinEWgroesserNull = T(1); + dLastMinEWgroesserNull = Scalar(1); printtime = writetime = firstwritetime = firstprinttime = 0; - pc = new T[params.N]; - ps = new T[params.N]; - tempRandom = new T[params.N + 1]; - BDz = new T[params.N]; - xmean = new T[params.N + 2]; - xmean[0] = params.N; - ++xmean; - xold = new T[params.N + 2]; - xold[0] = params.N; - ++xold; - xBestEver = new T[params.N + 3]; - xBestEver[0] = params.N; - ++xBestEver; - xBestEver[params.N] = std::numeric_limits::max(); - output = new T[params.N + 2]; - output[0] = params.N; - ++output; - rgD = new T[params.N]; - C = new T *[params.N]; - B = new T *[params.N]; - publicFitness = new T[params.lambda]; - functionValues = new T[params.lambda]; + pc = new Scalar[params.N]; + ps = new Scalar[params.N]; + tempRandom = new Scalar[params.N + 1]; + BDz = new Scalar[params.N]; + xmean = new Scalar[params.N]; + // xmean[0] = params.N; + // ++xmean; + xold = new Scalar[params.N]; + // xold[0] = params.N; + // ++xold; + xBestEver.x.resize((size_t)params.N); + // xBestEver = new Scalar[params.N + 3]; + // xBestEver[0] = params.N; + // ++xBestEver; + // xBestEver[params.N] = std::numeric_limits::max(); + xBestEver.value = std::numeric_limits::max(); + // output = new Scalar[params.N + 2]; + // output[0] = params.N; + // ++output; + rgD = new Scalar[params.N]; + C = new Scalar *[params.N]; + B = new Scalar *[params.N]; + publicFitness = new Ordered[params.lambda]; + functionValues = new Ordered[params.lambda]; // functionValues[0] = params.lambda; // ++functionValues; - const int historySize = 10 + (int)ceil(3. * 10. * params.N / params.lambda); - funcValueHistory = new T[historySize + 1]; - funcValueHistory[0] = (T)historySize; - funcValueHistory++; + historySize = + 10 + (size_t)ceil(3. * 10. * (double)params.N / params.lambda); + funcValueHistory = new Ordered[historySize]; + // funcValueHistory[0] = (Scalar)historySize; + // funcValueHistory++; for (int i = 0; i < params.N; ++i) { - C[i] = new T[i + 1]; - B[i] = new T[params.N]; + C[i] = new Scalar[i + 1]; + B[i] = new Scalar[params.N]; } - index = new int[params.lambda]; - for (int i = 0; i < params.lambda; ++i) + index = new size_t[params.lambda]; + for (size_t i = 0; i < params.lambda; ++i) index[i] = i; - population = new T *[params.lambda]; + population = new Individual[params.lambda]; for (int i = 0; i < params.lambda; ++i) { - population[i] = new T[params.N + 2]; - population[i][0] = params.N; - population[i]++; - for (int j = 0; j < params.N; j++) - population[i][j] = 0.0; + population[i].x.resize((size_t)params.N); + // population[i] = new Scalar[params.N + 2]; + // population[i][0] = params.N; + // population[i]++; + // for (int j = 0; j < params.N; j++) + // population[i][j] = 0.0; + population[i].value = std::numeric_limits::max(); } // initialize newed space for (int i = 0; i < params.lambda; i++) { - functionValues[i] = std::numeric_limits::max(); + functionValues[i] = std::numeric_limits::max(); } for (int i = 0; i < historySize; i++) { - funcValueHistory[i] = std::numeric_limits::max(); + funcValueHistory[i] = std::numeric_limits::max(); } for (int i = 0; i < params.N; ++i) for (int j = 0; j < i; ++j) C[i][j] = B[i][j] = B[j][i] = 0.; for (int i = 0; i < params.N; ++i) { - B[i][i] = T(1); - C[i][i] = rgD[i] = params.rgInitialStds[i] * std::sqrt(params.N / trace); + B[i][i] = Scalar(1); + C[i][i] = rgD[i] = + params.rgInitialStds[i] * std::sqrt((double)params.N / trace); C[i][i] *= C[i][i]; - pc[i] = ps[i] = T(0); + pc[i] = ps[i] = Scalar(0); } - minEW = minElement(rgD, params.N); + minEW = minElement(rgD, (size_t)params.N); minEW = minEW * minEW; - maxEW = maxElement(rgD, params.N); + maxEW = maxElement(rgD, (size_t)params.N); maxEW = maxEW * maxEW; maxdiagC = C[0][0]; @@ -1110,7 +1140,7 @@ template class CMAES { * @return A pointer to a "population" of lambda N-dimensional multivariate * normally distributed samples. */ - T *const *samplePopulation() { + Individual *samplePopulation() { bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; // calculate eigensystem @@ -1120,8 +1150,8 @@ template class CMAES { else { for (int i = 0; i < params.N; ++i) rgD[i] = std::sqrt(C[i][i]); - minEW = square(minElement(rgD, params.N)); - maxEW = square(maxElement(rgD, params.N)); + minEW = square(minElement(rgD, (size_t)params.N)); + maxEW = square(maxElement(rgD, (size_t)params.N)); eigensysIsUptodate = true; eigenTimings.start(); } @@ -1131,16 +1161,16 @@ template class CMAES { for (int iNk = 0; iNk < params.lambda; ++iNk) { // generate scaled random vector D*z - T *rgrgxink = population[iNk]; - for (int i = 0; i < params.N; ++i) + std::vector &rgrgxink = population[iNk].x; + for (size_t i = 0; i < params.N; ++i) if (diag) rgrgxink[i] = xmean[i] + sigma * rgD[i] * rand.gauss(); else tempRandom[i] = rgD[i] * rand.gauss(); if (!diag) - for (int i = 0; i < params.N; ++i) // add mutation sigma*B*(D*z) + for (size_t i = 0; i < params.N; ++i) // add mutation sigma*B*(D*z) { - T sum = 0.0; + Scalar sum = 0.0; for (int j = 0; j < params.N; ++j) sum += B[i][j] * tempRandom[j]; rgrgxink[i] = xmean[i] + sigma * sum; @@ -1163,11 +1193,10 @@ template class CMAES { * must hold. * @return A pointer to the resampled "population". */ - T *const *reSampleSingle(int i) { - T *x; + Individual *reSampleSingle(int i) { assert(i >= 0 && i < params.lambda && "reSampleSingle(): index must be between 0 and sp.lambda"); - x = population[i]; + std::vector &x = population[i].x; addMutation(x); return population; } @@ -1186,9 +1215,9 @@ template class CMAES { * @return A pointer to the resampled solution vector, equals input x for * x != NULL on input. */ - T *sampleSingleInto(T *x) { - if (!x) - x = new T[params.N]; + std::vector &sampleSingleInto(std::vector &x) { + // if (!x) + // x = new Scalar[params.N]; addMutation(x); return x; } @@ -1202,7 +1231,7 @@ template class CMAES { * sampled a new value. * @return A pointer to the resampled "population" member. */ - T const *reSampleSingleOld(T *x) { + Scalar const *reSampleSingleOld(Scalar *x) { assert(x && "reSampleSingleOld(): Missing input x"); addMutation(x); return x; @@ -1221,9 +1250,9 @@ template class CMAES { * @return A pointer to the perturbed solution vector, equals input x for * x != NULL. */ - T *perturbSolutionInto(T *x, T const *pxmean, T eps) { + Scalar *perturbSolutionInto(Scalar *x, Scalar const *pxmean, Scalar eps) { if (!x) - x = new T[params.N]; + x = new Scalar[params.N]; assert(pxmean && "perturbSolutionInto(): pxmean was not given"); addMutation(x, eps); return x; @@ -1236,7 +1265,7 @@ template class CMAES { * @param fitnessValues An array of \f$\lambda\f$ function values. * @return Mean value of the new distribution. */ - T *updateDistribution(const T *fitnessValues) { + Scalar *updateDistribution(const Ordered *fitnessValues) { const int N = params.N; bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; @@ -1252,16 +1281,27 @@ template class CMAES { params.logStream << "updateDistribution(): unexpected state" << std::endl; // assign function values - for (int i = 0; i < params.lambda; ++i) - population[i][N] = functionValues[i] = fitnessValues[i]; + for (int i = 0; i < params.lambda; ++i) { + // std::cout << i << ": " << fitnessValues[i] << std::endl; + population[i].value = functionValues[i] = fitnessValues[i]; + } // Generate index - sortIndex(fitnessValues, index, params.lambda); + // std::cout << "["; + // for (size_t i = 0; i < params.lambda; ++i) { + // std::cout << index[i] << ","; + // } + // std::cout << "] => ["; + sortIndex(fitnessValues, index, (size_t)params.lambda); + // for (size_t i = 0; i < params.lambda; ++i) { + // std::cout << index[i] << ","; + // } + // std::cout << "]" << std::endl; // Test if function values are identical, escape flat fitness if (fitnessValues[index[0]] == fitnessValues[index[(int)params.lambda / 2]]) { - sigma *= std::exp(T(0.2) + params.cs / params.damps); + sigma *= std::exp(Scalar(0.2) + params.cs / params.damps); if (params.logWarnings) { params.logStream << "Warning: sigma increased due to equal function values" @@ -1271,30 +1311,37 @@ template class CMAES { } // update function value history - for (int i = (int)*(funcValueHistory - 1) - 1; i > 0; --i) + for (size_t i = historySize - 1; i > 0; --i) funcValueHistory[i] = funcValueHistory[i - 1]; funcValueHistory[0] = fitnessValues[index[0]]; // update xbestever - if (xBestEver[N] > population[index[0]][N] || gen == 1) - for (int i = 0; i <= N; ++i) { - xBestEver[i] = population[index[0]][i]; - xBestEver[N + 1] = countevals; - } + // std::cout << population[index[0]].value << std::endl; + if (xBestEver.value > population[index[0]].value || gen == 1) { + // std::cout << xBestEver.value << " -> " << population[index[0]].value << + // std::endl; + xBestEver = population[index[0]]; + evalsBestEver = countevals; + // for (int i = 0; i < N; ++i) { + // // xBestEver.x[i] = population[index[0]].x[i]; /* TODO: use + // std::vector copy assign */ + // // xBestEver.value = population[index[0]].value; + // } + } - const T sqrtmueffdivsigma = std::sqrt(params.mueff) / sigma; + const Scalar sqrtmueffdivsigma = std::sqrt(params.mueff) / sigma; // calculate xmean and rgBDz~N(0,C) - for (int i = 0; i < N; ++i) { + for (size_t i = 0; i < N; ++i) { xold[i] = xmean[i]; xmean[i] = 0.; for (int iNk = 0; iNk < params.mu; ++iNk) - xmean[i] += params.weights[iNk] * population[index[iNk]][i]; + xmean[i] += params.weights[iNk] * population[index[iNk]].x[i]; BDz[i] = sqrtmueffdivsigma * (xmean[i] - xold[i]); } // calculate z := D^(-1)* B^(-1)* rgBDz into rgdTmp for (int i = 0; i < N; ++i) { - T sum; + Scalar sum; if (diag) sum = BDz[i]; else { @@ -1306,15 +1353,15 @@ template class CMAES { } // cumulation for sigma (ps) using B*z - const T sqrtFactor = std::sqrt(params.cs * (T(2) - params.cs)); - const T invps = T(1) - params.cs; + const Scalar sqrtFactor = std::sqrt(params.cs * (Scalar(2) - params.cs)); + const Scalar invps = Scalar(1) - params.cs; for (int i = 0; i < N; ++i) { - T sum; + Scalar sum; if (diag) sum = tempRandom[i]; else { - sum = T(0); - T *Bi = B[i]; + sum = Scalar(0); + Scalar *Bi = B[i]; for (int j = 0; j < N; ++j) sum += Bi[j] * tempRandom[j]; } @@ -1322,20 +1369,21 @@ template class CMAES { } // calculate norm(ps)^2 - T psxps(0); + Scalar psxps(0); for (int i = 0; i < N; ++i) { - const T &rgpsi = ps[i]; + const Scalar &rgpsi = ps[i]; psxps += rgpsi * rgpsi; } // cumulation for covariance matrix (pc) using B*D*z~N(0,C) int hsig = std::sqrt(psxps) / - std::sqrt(T(1) - std::pow(T(1) - params.cs, T(2) * gen)) / + std::sqrt(Scalar(1) - std::pow(Scalar(1) - params.cs, + Scalar(2) * (Scalar)gen)) / chiN < - T(1.4) + T(2) / (N + 1); - const T ccumcovinv = 1. - params.ccumcov; - const T hsigFactor = - hsig * std::sqrt(params.ccumcov * (T(2) - params.ccumcov)); + Scalar(1.4) + Scalar(2) / (N + 1); + const Scalar ccumcovinv = 1. - params.ccumcov; + const Scalar hsigFactor = + hsig * std::sqrt(params.ccumcov * (Scalar(2) - params.ccumcov)); for (int i = 0; i < N; ++i) pc[i] = ccumcovinv * pc[i] + hsigFactor * BDz[i]; @@ -1343,8 +1391,8 @@ template class CMAES { adaptC2(hsig); // update of sigma - sigma *= - std::exp(((std::sqrt(psxps) / chiN) - T(1)) * params.cs / params.damps); + sigma *= std::exp(((std::sqrt(psxps) / chiN) - Scalar(1)) * params.cs / + params.damps); state = UPDATED; return xmean; @@ -1355,16 +1403,13 @@ template class CMAES { * @param key Key of the requested scalar. * @return The desired value. */ - T get(GetScalar key) { + Scalar get(GetScalar key) { switch (key) { case AxisRatio: - return maxElement(rgD, params.N) / minElement(rgD, params.N); + return maxElement(rgD, (size_t)params.N) / + minElement(rgD, (size_t)params.N); case Eval: return countevals; - case Fitness: - return functionValues[index[0]]; - case FBestEver: - return xBestEver[params.N]; case Generation: return gen; case MaxEval: @@ -1390,6 +1435,22 @@ template class CMAES { } } + /** + * Request a parameter of object function type from CMA-ES. + * @param key Key of the requested value. + * @return The desired value. + */ + Ordered getValue(GetScalar key) { + switch (key) { + case Fitness: + return functionValues[index[0]]; + case FBestEver: + return xBestEver.value; + default: + throw std::runtime_error("getValue(): No match found for key"); + } + } + /** * Request a vector parameter from CMA-ES. * @param key Key of the requested vector. @@ -1397,26 +1458,28 @@ template class CMAES { * overwritten during the next call to any member functions other * than get(). */ - const T *getPtr(GetVector key) { + const std::vector getPtr(GetVector key) { switch (key) { case DiagC: { - for (int i = 0; i < params.N; ++i) - output[i] = C[i][i]; - return output; + std::vector o((size_t)params.N); + for (size_t i = 0; i < params.N; ++i) + o[i] = C[i][i]; + return o; } case DiagD: - return rgD; + return std::vector(rgD, rgD + params.N); case StdDev: { - for (int i = 0; i < params.N; ++i) - output[i] = sigma * std::sqrt(C[i][i]); - return output; + std::vector o((size_t)params.N); + for (size_t i = 0; i < params.N; ++i) + o[i] = sigma * std::sqrt(C[i][i]); + return o; } case XBestEver: - return xBestEver; + return xBestEver.x; case XBest: - return population[index[0]]; + return population[index[0]].x; case XMean: - return xmean; + return std::vector(xmean, xmean + params.N); default: throw std::runtime_error("getPtr(): No match found for key"); } @@ -1429,7 +1492,7 @@ template class CMAES { * writing access to its elements. The memory must be explicitly * released using delete[]. */ - T *getNew(GetVector key) { return getInto(key, 0); } + Scalar *getNew(GetVector key) { return getInto(key, 0); } /** * Request a vector parameter from CMA-ES. @@ -1438,11 +1501,11 @@ template class CMAES { * written into. For mem == NULL new memory is allocated as with * calling getNew() and must be released by the user at some point. */ - T *getInto(GetVector key, T *res) { - T const *res0 = getPtr(key); + Scalar *getInto(GetVector key, Scalar *res) { + const std::vector res0 = getPtr(key); if (!res) - res = new T[params.N]; - for (int i = 0; i < params.N; ++i) + res = new Scalar[params.N]; + for (size_t i = 0; i < params.N; ++i) res[i] = res0[i]; return res; } @@ -1456,7 +1519,8 @@ template class CMAES { * @return Does any stop criterion match? */ bool testForTermination() { - T range, fac; + Ordered range; + Scalar fac; int iAchse, iKoo; int diag = params.diagonalCov == 1 || params.diagonalCov >= gen; int N = params.N; @@ -1476,11 +1540,11 @@ template class CMAES { // TolFun range = std::max(maxElement(funcValueHistory, - (int)std::min(gen, *(funcValueHistory - 1))), - maxElement(functionValues, params.lambda)) - + std::min((size_t)gen, historySize)), + maxElement(functionValues, (size_t)params.lambda)) - std::min(minElement(funcValueHistory, - (int)std::min(gen, *(funcValueHistory - 1))), - minElement(functionValues, params.lambda)); + std::min((size_t)gen, historySize)), + minElement(functionValues, (size_t)params.lambda)); if (gen > 0 && range <= params.stopTolFun) { message << "TolFun: function value differences " << range @@ -1488,9 +1552,9 @@ template class CMAES { } // TolFunHist - if (gen > *(funcValueHistory - 1)) { - range = maxElement(funcValueHistory, (int)*(funcValueHistory - 1)) - - minElement(funcValueHistory, (int)*(funcValueHistory - 1)); + if (gen > historySize) { + range = maxElement(funcValueHistory, historySize) - + minElement(funcValueHistory, historySize); if (range <= params.stopTolFunHist) message << "TolFunHist: history of function value changes " << range << " stopTolFunHist=" << params.stopTolFunHist << std::endl; @@ -1547,7 +1611,7 @@ template class CMAES { // Component of xmean is not changed anymore for (iKoo = 0; iKoo < N; ++iKoo) { if (xmean[iKoo] == - xmean[iKoo] + sigma * std::sqrt(C[iKoo][iKoo]) / T(5)) { + xmean[iKoo] + sigma * std::sqrt(C[iKoo][iKoo]) / Scalar(5)) { message << "NoEffectCoordinate: standard deviation 0.2*" << (sigma * std::sqrt(C[iKoo][iKoo])) << " in coordinate " << iKoo << " without effect" << std::endl; @@ -1595,7 +1659,7 @@ template class CMAES { /** * Conducts the eigendecomposition of C into B and D such that - * \f$C = B \cdot D \cdot D \cdot B^T\f$ and \f$B \cdot B^T = I\f$ + * \f$C = B \cdot D \cdot D \cdot B^Scalar\f$ and \f$B \cdot B^Scalar = I\f$ * and D diagonal and positive. * @param force For force == true the eigendecomposion is conducted even if * eigenvector and values seem to be up to date. @@ -1623,8 +1687,8 @@ template class CMAES { // find largest and smallest eigenvalue, they are supposed to be sorted // anyway - minEW = minElement(rgD, params.N); - maxEW = maxElement(rgD, params.N); + minEW = minElement(rgD, (size_t)params.N); + maxEW = maxElement(rgD, (size_t)params.N); if (doCheckEigen) // needs O(n^3)! writes, in case, error message in error // file @@ -1643,7 +1707,7 @@ template class CMAES { * @param newxmean new mean, if it is NULL, it will be set to the current mean * @return new mean */ - T const *setMean(const T *newxmean) { + Scalar const *setMean(const Scalar *newxmean) { assert(state != SAMPLED && "setMean: mean cannot be set inbetween the calls" "of samplePopulation and updateDistribution"); diff --git a/include/cma-es/parameters.h b/include/cma-es/parameters.h index 501aa57..f7fae99 100644 --- a/include/cma-es/parameters.h +++ b/include/cma-es/parameters.h @@ -7,48 +7,48 @@ #include #include -template class CMAES; +template class CMAES; /** * @class Parameters * Holds all parameters that can be adjusted by the user. */ -template class Parameters { - friend class CMAES; +template class Parameters { + friend class CMAES; public: /* Input parameters. */ //! Problem dimension, must stay constant. int N; //! Initial search space vector. - T *xstart; + Scalar *xstart; //! A typical value for a search space vector. - T *typicalX; + Scalar *typicalX; //! Indicates that the typical x is the initial point. bool typicalXcase; //! Initial standard deviations. - T *rgInitialStds; - T *rgDiffMinChange; + Scalar *rgInitialStds; + Scalar *rgDiffMinChange; /* Termination parameters. */ //! Maximal number of objective function evaluations. - T stopMaxFunEvals; - T facmaxeval; + Scalar stopMaxFunEvals; + Scalar facmaxeval; //! Maximal number of iterations. - T stopMaxIter; + Scalar stopMaxIter; //! Minimal fitness value. Only activated if flg is true. struct { bool flg; - T val; + Scalar val; } stStopFitness; //! Minimal value difference. - T stopTolFun; + Comparable stopTolFun; //! Minimal history value difference. - T stopTolFunHist; + Comparable stopTolFunHist; //! Minimal search space step size. - T stopTolX; + Scalar stopTolX; //! Defines the maximal condition number. - T stopTolUpXFactor; + Scalar stopTolUpXFactor; /* internal evolution strategy parameters */ /** @@ -60,37 +60,37 @@ template class Parameters { * Number of individuals used to recompute the mean. */ int mu; - T mucov; + Scalar mucov; /** * Variance effective selection mass, should be lambda/4. */ - T mueff; + Scalar mueff; /** * Weights used to recombinate the mean sum up to one. */ - T *weights; + Scalar *weights; /** * Damping parameter for step-size adaption, d = inifinity or 0 means adaption * is turned off, usually close to one. */ - T damps; + Scalar damps; /** * cs^-1 (approx. n/3) is the backward time horizon for the evolution path * ps and larger than one. */ - T cs; - T ccumcov; + Scalar cs; + Scalar ccumcov; /** * ccov^-1 (approx. n/4) is the backward time horizon for the evolution path * pc and larger than one. */ - T ccov; - T diagonalCov; + Scalar ccov; + Scalar diagonalCov; struct { - T modulo; - T maxtime; + Scalar modulo; + Scalar maxtime; } updateCmode; - T facupdateCmode; + Scalar facupdateCmode; /** * Determines the method used to initialize the weights. @@ -120,7 +120,7 @@ template class Parameters { weightMode(UNINITIALIZED_WEIGHTS), resumefile(""), logWarnings(false), logStream(std::cerr) { stStopFitness.flg = false; - stStopFitness.val = -std::numeric_limits::max(); + stStopFitness.val = -std::numeric_limits::max(); updateCmode.modulo = -1; updateCmode.maxtime = -1; } @@ -150,18 +150,19 @@ template class Parameters { * available, must be defined here or you have to set the * member manually. * @param inxstart Initial point in search space \f$x_0\f$, default (NULL) is - * \f$(0.5,\ldots,0.5)^T + N(0, initialStdDev^2) \in R^N\f$. - * This must be an array of size \f$N\f$. + * \f$(0.5,\ldots,0.5)^Scalar + N(0, initialStdDev^2) \in + * R^N\f$. This must be an array of size \f$N\f$. * @param inrgsigma Coordinatewise initial standard deviation of the sample * distribution (\f$\sigma \cdot \sqrt{C_{ii}} = * initialStdDev[i]\f$). The expected initial distance * between initialX and the optimum per coordinate should be * roughly initialStdDev. The entries should not differ by * several orders of magnitude. Default (NULL) is - * \f$(0.3,\ldots,0.3)^T \in R^N\f$. This must be an array of - * size \f$N\f$. + * \f$(0.3,\ldots,0.3)^Scalar \in R^N\f$. This must be an + * array of size \f$N\f$. */ - void init(int dimension = 0, const T *inxstart = 0, const T *inrgsigma = 0) { + void init(int dimension = 0, const Scalar *inxstart = 0, + const Scalar *inrgsigma = 0) { if (logWarnings) { if (!(xstart || inxstart || typicalX)) logStream << "Warning: initialX undefined. typicalX = 0.5...0.5." @@ -182,7 +183,7 @@ template class Parameters { diagonalCov = 0; // default is 0, but this might change in future if (!xstart) { - xstart = new T[N]; + xstart = new Scalar[N]; if (inxstart) { for (int i = 0; i < N; ++i) xstart[i] = inxstart[i]; @@ -198,13 +199,13 @@ template class Parameters { } if (!rgInitialStds) { - rgInitialStds = new T[N]; + rgInitialStds = new Scalar[N]; if (inrgsigma) { for (int i = 0; i < N; ++i) rgInitialStds[i] = inrgsigma[i]; } else { for (int i = 0; i < N; ++i) - rgInitialStds[i] = T(0.3); + rgInitialStds[i] = Scalar(0.3); } } @@ -218,7 +219,7 @@ template class Parameters { if (xstart) delete[] xstart; if (p.xstart) { - xstart = new T[N]; + xstart = new Scalar[N]; for (int i = 0; i < N; i++) xstart[i] = p.xstart[i]; } @@ -226,7 +227,7 @@ template class Parameters { if (typicalX) delete[] typicalX; if (p.typicalX) { - typicalX = new T[N]; + typicalX = new Scalar[N]; for (int i = 0; i < N; i++) typicalX[i] = p.typicalX[i]; } @@ -236,7 +237,7 @@ template class Parameters { if (rgInitialStds) delete[] rgInitialStds; if (p.rgInitialStds) { - rgInitialStds = new T[N]; + rgInitialStds = new Scalar[N]; for (int i = 0; i < N; i++) rgInitialStds[i] = p.rgInitialStds[i]; } @@ -244,7 +245,7 @@ template class Parameters { if (rgDiffMinChange) delete[] rgDiffMinChange; if (p.rgDiffMinChange) { - rgDiffMinChange = new T[N]; + rgDiffMinChange = new Scalar[N]; for (int i = 0; i < N; i++) rgDiffMinChange[i] = p.rgDiffMinChange[i]; } @@ -269,7 +270,7 @@ template class Parameters { if (weights) delete[] weights; if (p.weights) { - weights = new T[mu]; + weights = new Scalar[mu]; for (int i = 0; i < mu; i++) weights[i] = p.weights[i]; } @@ -294,25 +295,26 @@ template class Parameters { * Supplements default parameter values. */ void supplementDefaults() { + Scalar n = (Scalar)N; if (lambda < 2) - lambda = 4 + (int)(3.0 * log((double)N)); + lambda = 4 + (int)(3.0 * log(n)); if (mu <= 0) mu = lambda / 2; if (!weights) setWeights(weightMode); if (cs > 0) - cs *= (mueff + 2.) / (N + mueff + 3.); + cs *= (mueff + 2.) / (n + mueff + 3.); if (cs <= 0 || cs >= 1) - cs = (mueff + 2.) / (N + mueff + 3.); + cs = (mueff + 2.) / (n + mueff + 3.); if (ccumcov <= 0 || ccumcov > 1) - ccumcov = 4. / (N + 4); + ccumcov = 4. / (n + 4); if (mucov < 1) mucov = mueff; - T t1 = 2. / ((N + 1.4142) * (N + 1.4142)); - T t2 = (2. * mueff - 1.) / ((N + 2.) * (N + 2.) + mueff); + Scalar t1 = 2. / ((n + 1.4142) * (n + 1.4142)); + Scalar t2 = (2. * mueff - 1.) / ((n + 2.) * (n + 2.) + mueff); t2 = (t2 > 1) ? 1 : t2; t2 = (1. / mucov) * t1 + (1. - 1. / mucov) * t2; if (ccov >= 0) @@ -321,31 +323,32 @@ template class Parameters { ccov = t2; if (diagonalCov < 0) - diagonalCov = 2 + 100. * N / sqrt((double)lambda); + diagonalCov = 2 + 100. * n / sqrt((double)lambda); if (stopMaxFunEvals <= 0) - stopMaxFunEvals = facmaxeval * 900 * (N + 3) * (N + 3); + stopMaxFunEvals = facmaxeval * 900 * (n + 3) * (n + 3); else stopMaxFunEvals *= facmaxeval; if (stopMaxIter <= 0) stopMaxIter = ceil((double)(stopMaxFunEvals / lambda)); - if (damps < T(0)) - damps = T(1); + if (damps < Scalar(0)) + damps = Scalar(1); damps = damps * - (T(1) + - T(2) * std::max(T(0), std::sqrt((mueff - T(1)) / (N + T(1))) - - T(1))) * - (T)std::max( - T(0.3), - T(1) - // modify for short runs - (T)N / (T(1e-6) + std::min(stopMaxIter, - stopMaxFunEvals / lambda))) + + (Scalar(1) + + Scalar(2) * std::max(Scalar(0), std::sqrt((mueff - Scalar(1)) / + (n + Scalar(1))) - + Scalar(1))) * + (Scalar)std::max( + Scalar(0.3), + Scalar(1) - // modify for short runs + n / (Scalar(1e-6) + + std::min(stopMaxIter, stopMaxFunEvals / lambda))) + cs; if (updateCmode.modulo < 0) - updateCmode.modulo = 1. / ccov / (double)N / 10.; + updateCmode.modulo = 1. / ccov / n / 10.; updateCmode.modulo *= facupdateCmode; if (updateCmode.maxtime < 0) updateCmode.maxtime = 0.20; // maximal 20% of CPU-time @@ -357,7 +360,7 @@ template class Parameters { void setWeights(Weights mode) { if (weights) delete[] weights; - weights = new T[mu]; + weights = new Scalar[mu]; switch (mode) { case LINEAR_WEIGHTS: for (int i = 0; i < mu; ++i) @@ -375,7 +378,7 @@ template class Parameters { } // normalize weights vector and set mueff - T s1 = 0, s2 = 0; + Scalar s1 = 0, s2 = 0; for (int i = 0; i < mu; ++i) { s1 += weights[i]; s2 += weights[i] * weights[i]; diff --git a/include/cma-es/random.h b/include/cma-es/random.h index 50dcd34..2ad5dad 100644 --- a/include/cma-es/random.h +++ b/include/cma-es/random.h @@ -9,10 +9,10 @@ */ template class Random { // variables for uniform() - long int startseed; - long int aktseed; - long int aktrand; - long int rgrand[32]; + long startseed; + long aktseed; + long aktrand; + long rgrand[32]; // variables for gauss() bool stored; T hold; @@ -21,18 +21,18 @@ template class Random { /** * @param seed use clock if 0 */ - Random(long unsigned seed = 0) : hold(0.0) { + Random(long seed = 0) : hold(0.0) { stored = false; if (seed < 1) { - long int t = 100 * time(0) + clock(); - seed = (long unsigned)(t < 0 ? -t : t); + long t = 100 * time(0) + (long)clock(); + seed = (long)(t < 0 ? -t : t); } start(seed); } /** * @param seed 0 == 1 */ - void start(long unsigned seed) { + void start(long seed) { stored = false; startseed = seed; if (seed < 1) diff --git a/include/cma-es/utils.h b/include/cma-es/utils.h index 94d349b..e6e1887 100644 --- a/include/cma-es/utils.h +++ b/include/cma-es/utils.h @@ -9,33 +9,37 @@ #include #include -template T square(T d) { return d * d; } +template inline Scalar square(Scalar d) { return d * d; } -template T maxElement(const T *rgd, int len) { +template +inline Scalar maxElement(const Scalar *rgd, size_t len) { return *std::max_element(rgd, rgd + len); } -template T minElement(const T *rgd, int len) { +template +inline Scalar minElement(const Scalar *rgd, size_t len) { return *std::min_element(rgd, rgd + len); } -template int maxIndex(const T *rgd, int len) { - return std::max_element(rgd, rgd + len) - rgd; +template +inline size_t maxIndex(const Scalar *rgd, size_t len) { + return (size_t)(std::max_element(rgd, rgd + len) - rgd); } -template int minIndex(const T *rgd, int len) { - return std::min_element(rgd, rgd + len) - rgd; +template +inline size_t minIndex(const Scalar *rgd, size_t len) { + return (size_t)(std::min_element(rgd, rgd + len) - rgd); } /** sqrt(a^2 + b^2) numerically stable. */ -template T myhypot(T a, T b) { - const T fabsa = std::fabs(a), fabsb = std::fabs(b); +template Scalar myhypot(Scalar a, Scalar b) { + const Scalar fabsa = std::fabs(a), fabsb = std::fabs(b); if (fabsa > fabsb) { - const T r = b / a; - return fabsa * std::sqrt(T(1) + r * r); - } else if (b != T(0)) { - const T r = a / b; - return fabsb * std::sqrt(T(1) + r * r); + const Scalar r = b / a; + return fabsa * std::sqrt(Scalar(1) + r * r); + } else if (b != Scalar(0)) { + const Scalar r = a / b; + return fabsb * std::sqrt(Scalar(1) + r * r); } else - return T(0); + return Scalar(0); } diff --git a/include/lexi.h b/include/lexi.h index 2536141..695e7c8 100644 --- a/include/lexi.h +++ b/include/lexi.h @@ -1,11 +1,25 @@ +#include +#include #include -template class LexiProduct { - using Product = std::vector; +/* +support all compaison +support minus +LexiProduct(0) <= any lexi <= LexiProduct(std::numeric_limits::max()) +*/ +template class LexiProduct { + using Product = std::vector; public: LexiProduct() : LexiProduct(0) {} - LexiProduct(size_t size) : prd(size) {} + LexiProduct(Scalar val) : prd{val} {} + LexiProduct(const Product &prd2) : prd(prd2) {} + LexiProduct(Product &&prd2) : prd(std::move(prd2)) {} + + LexiProduct &operator=(const Scalar value) { + prd = std::vector{value}; + return *this; + } LexiProduct &operator=(const Product &prd2) { prd = prd2; return *this; @@ -15,37 +29,57 @@ template class LexiProduct { return *this; } - bool operator==(const LexiProduct &lexi2) { + bool operator==(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) == 0; } - bool operator!=(const LexiProduct &lexi2) { + bool operator!=(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) != 0; } - bool operator<(const LexiProduct &lexi2) { + bool operator<(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) < 0; } - bool operator<=(const LexiProduct &lexi2) { + bool operator<=(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) <= 0; } - bool operator>(const LexiProduct &lexi2) { + bool operator>(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) > 0; } - bool operator>=(const LexiProduct &lexi2) { + bool operator>=(const LexiProduct &lexi2) const { return lexicmp(*this, lexi2) >= 0; } + LexiProduct operator-(const LexiProduct &lexi2) const { + const size_t n = std::min(this->prd.size(), lexi2.prd.size()); + Product d(n); + for (size_t i = 0; i < n; ++i) { + d[i] = this->prd[i] - lexi2.prd[i]; + } + return LexiProduct(std::move(d)); + } + private: Product prd; - int lexicmp(const LexiProduct &lhs, const LexiProduct &rhs) { - const int lencmp = lhs.prd.size() - rhs.prd.size(); + Scalar lexicmp(const LexiProduct &lhs, + const LexiProduct &rhs) const { + const Scalar lencmp = (Scalar)(lhs.prd.size() - rhs.prd.size()); const size_t n = (lencmp < 0 ? lhs : rhs).prd.size(); for (size_t i = 0; i < n; ++i) { - const T d = lhs.prd[i] - rhs.prd[i]; + const Scalar d = lhs.prd[i] - rhs.prd[i]; if (d != 0) { return d; } } return lencmp; } -}; + + friend std::ostream &operator<<(std::ostream &os, const LexiProduct &lexi) { + const size_t n = lexi.prd.size(); + os << "["; + for (size_t i = 0; i < n - 1; ++i) { + os << lexi.prd[i] << ","; + } + os << lexi.prd[n - 1] << "]"; + return os; + } +}; \ No newline at end of file diff --git a/src/conversion.cpp b/src/conversion.cpp index bbb54e8..03d94be 100644 --- a/src/conversion.cpp +++ b/src/conversion.cpp @@ -1,9 +1,11 @@ -#define _USE_MATH_DEFINES #include #include #include namespace color { +double sgn(double val); +constexpr inline double pow3(double x); + double sgn(double val) { return (0. < val) - (val < 0.); } constexpr inline double pow3(double x) { return x * x * x; } diff --git a/test/example1.cpp b/test/example1.cpp index a1253c8..9c91442 100644 --- a/test/example1.cpp +++ b/test/example1.cpp @@ -6,24 +6,31 @@ #include #include +#include #include +#include + +using Y = LexiProduct; +// using Y = double; /** * The objective (fitness) function to be minized. "cigtab" function. */ -double fitfun(double const *x, int N) { +Y fitfun(const std::vector &x, int N) { double sum = 1e4 * x[0] * x[0] + 1e-4 * x[1] * x[1]; - for (int i = 2; i < N; ++i) + for (size_t i = 2; i < N; ++i) sum += x[i] * x[i]; - return sum; + return std::vector{0., sum}; } /** * The optimization loop. */ int main(int, char **) { - CMAES evo; - double *arFunvals, *const *pop, *xfinal; + CMAES evo; + Y *arFunvals; + Individual *pop; + double *xfinal; // Initialize everything const int dim = 22; @@ -33,7 +40,9 @@ int main(int, char **) { double stddev[dim]; for (int i = 0; i < dim; i++) stddev[i] = 0.3; - Parameters parameters; + Parameters parameters; + parameters.stopTolFun = std::vector{0, 1e-12}; + parameters.stopTolFunHist = std::vector{0, 1e-13}; // TODO Adjust parameters here parameters.init(dim, xstart, stddev); arFunvals = evo.init(parameters); @@ -60,19 +69,26 @@ int main(int, char **) { */ // evaluate the new search points using fitfun from above - for (int i = 0; i < evo.get(CMAES::Lambda); ++i) - arFunvals[i] = fitfun(pop[i], (int)evo.get(CMAES::Dimension)); + for (int i = 0; i < evo.get(CMAES::Lambda); ++i) + arFunvals[i] = + fitfun(pop[i].x, (int)evo.get(CMAES::Dimension)); // update the search distribution used for sampleDistribution() - evo.updateDistribution(arFunvals); + double const *xmean = evo.updateDistribution(arFunvals); + // for (size_t i = 0; i < (int)evo.get(CMAES::Dimension); ++i) { + // std::cout << xmean[i] << ","; + // } + // std::cout << std::endl; } std::cout << "Stop:" << std::endl << evo.getStopMessage(); - evo.writeToFile(CMAES::WKResume, + evo.writeToFile(CMAES::WKResume, "resumeevo1.dat"); // write resumable state of CMA-ES // get best estimator for the optimum, xmean - xfinal = - evo.getNew(CMAES::XMean); // "XBestEver" might be used as well + xfinal = evo.getNew(CMAES::XMean); // "XBestEver" + // might be + // used as + // well // do something with final solution and finally release memory delete[] xfinal; diff --git a/test/example2.cpp b/test/example2.cpp index 6ad1915..a9a7308 100644 --- a/test/example2.cpp +++ b/test/example2.cpp @@ -13,31 +13,34 @@ #include double **OrthogonalBasis(int DIM); -double f_rosenbrock(double const *x); -double f_rand(double const *x); -double f_constant(double const *x); -double f_kugelmin1(double const *x); -double f_sphere(double const *x); -double f_stepsphere(double const *x); -double f_cigar(double const *x); -double f_cigtab(double const *x); -double f_tablet(double const *x); -double f_elli(double const *x); -double f_ellirot(double const *x); -double f_elli100(double const *x); -double f_ellinumtest(double const *x); -double f_parabR(double const *x); -double f_sharpR(double const *x); -double f_diffpow(double const *x); -double f_diffpowrot(double const *x); -double f_gleichsys5(double const *x); - -double *optimize(double (*pFun)(double const *), int number_of_restarts, +double f_rosenbrock(const std::vector &x); +double f_rand(const std::vector &x); +double f_constant(const std::vector &x); +double f_kugelmin1(const std::vector &x); +double f_sphere(const std::vector &x); +double f_stepsphere(const std::vector &x); +double f_cigar(const std::vector &x); +double f_cigtab(const std::vector &x); +double f_tablet(const std::vector &x); +double f_elli(const std::vector &x); +double f_ellirot(const std::vector &x); +double f_elli100(const std::vector &x); +double f_ellinumtest(const std::vector &x); +double f_parabR(const std::vector &x); +double f_sharpR(const std::vector &x); +double f_diffpow(const std::vector &x); +double f_diffpowrot(const std::vector &x); +double f_gleichsys5(const std::vector &x); + +double *optimize(double (*pFun)(const std::vector &), + int number_of_restarts, double increment_factor_for_population_size); int main(int, char **) { - typedef double (*pfun_t)(double const *); - pfun_t rgpFun[99]; // array (range) of pointer to objective function + typedef double (*pfun_t)(const std::vector &); + size_t maxTest = 99; + std::vector rgpFun( + maxTest); // vector (range) of pointer to objective function double incpopsize = 2; double *x; @@ -61,11 +64,15 @@ int main(int, char **) { rgpFun[22] = f_ellirot; rgpFun[23] = f_diffpowrot; - int nb = 5; int nbrestarts = 0; - - // Optimize function - x = optimize(rgpFun[nb], nbrestarts, incpopsize); + for (size_t nb = 0; nb < maxTest; ++nb) { + // Optimize function + pfun_t f = rgpFun[nb]; + if (f == nullptr) { + continue; + } + x = optimize(f, nbrestarts, incpopsize); + } // here we could utilize the solution x, and finally free memory delete[] x; @@ -77,11 +84,11 @@ int main(int, char **) { * Somewhat extended interface for optimizing pFun with CMAES implementing a * restart procedure with increasing population size. */ -double *optimize(double (*pFun)(double const *), int nrestarts, +double *optimize(double (*pFun)(const std::vector &), int nrestarts, double incpopsize) { - CMAES evo; // the optimizer - double *const *pop; // sampled population - double *fitvals; // objective function values of sampled population + CMAES evo; // the optimizer + Individual *pop; // sampled population + double *fitvals; // objective function values of sampled population double fbestever = 0, *xbestever = NULL; // store best solution double fmean; int irun, @@ -100,7 +107,7 @@ double *optimize(double (*pFun)(double const *), int nrestarts, double stddev[dim]; for (int i = 0; i < dim; i++) stddev[i] = 0.3; - Parameters parameters; + Parameters parameters; // You can resume a previous run by specifying a file that contains the // resume data: // parameters.resumefile = "resumeevo2.dat"; @@ -128,7 +135,7 @@ double *optimize(double (*pFun)(double const *), int nrestarts, */ // Compute fitness value for each candidate solution - for (int i = 0; i < evo.get(CMAES::PopSize); ++i) { + for (int i = 0; i < evo.get(CMAES::PopSize); ++i) { /* You may resample the solution i until it lies within the feasible domain here, e.g. until it satisfies given box constraints (variable boundaries). The function @@ -143,7 +150,7 @@ double *optimize(double (*pFun)(double const *), int nrestarts, /* while (!is_feasible(pop[i])) evo.reSampleSingle(i); */ - fitvals[i] = (*pFun)(pop[i]); + fitvals[i] = (*pFun)(pop[i].x); } // update search distribution @@ -152,36 +159,41 @@ double *optimize(double (*pFun)(double const *), int nrestarts, fflush(stdout); } - lambda = (int)(incpopsize * - evo.get(CMAES::Lambda)); // needed for the restart - countevals = (int)evo.get(CMAES::Eval); // dito + lambda = + (int)(incpopsize * + evo.get(CMAES::Lambda)); // needed for the restart + countevals = (int)evo.get(CMAES::Eval); // dito // print some "final" output - std::cout << (int)evo.get(CMAES::Generation) << " generations, " - << (int)evo.get(CMAES::Eval) << " fevals (" - << evo.eigenTimings.totaltime - << " sec): f(x)=" << evo.get(CMAES::Fitness) << std::endl + std::cout << (int)evo.get(CMAES::Generation) + << " generations, " << (int)evo.get(CMAES::Eval) + << " fevals (" << evo.eigenTimings.totaltime + << " sec): f(x)=" << evo.getValue(CMAES::Fitness) + << std::endl << " (axis-ratio=" - << evo.get(CMAES::MaxAxisLength) / - evo.get(CMAES::MinAxisLength) - << ", max/min-stddev=" << evo.get(CMAES::MaxStdDev) << "/" - << evo.get(CMAES::MinStdDev) << ")" << std::endl + << evo.get(CMAES::MaxAxisLength) / + evo.get(CMAES::MinAxisLength) + << ", max/min-stddev=" + << evo.get(CMAES::MaxStdDev) << "/" + << evo.get(CMAES::MinStdDev) << ")" << std::endl << "Stop (run " << (irun + 1) << "):" << std::endl << evo.getStopMessage(); // write resume data - evo.writeToFile(CMAES::WKResume, "resumeevo2.dat"); + evo.writeToFile(CMAES::WKResume, "resumeevo2.dat"); // keep best ever solution - if (irun == 0 || evo.get(CMAES::FBestEver) < fbestever) { - fbestever = evo.get(CMAES::FBestEver); - xbestever = evo.getInto(CMAES::XBestEver, + if (irun == 0 || + evo.getValue(CMAES::FBestEver) < fbestever) { + fbestever = evo.getValue(CMAES::FBestEver); + xbestever = evo.getInto(CMAES::XBestEver, xbestever); // allocates memory if needed } // best estimator for the optimum is xmean, therefore check - if ((fmean = (*pFun)(evo.getPtr(CMAES::XMean))) < fbestever) { + const std::vector xmean = evo.getPtr(CMAES::XMean); + if ((fmean = (*pFun)(xmean)) < fbestever) { fbestever = fmean; - xbestever = evo.getInto(CMAES::XMean, xbestever); + xbestever = evo.getInto(CMAES::XMean, xbestever); } // abandon restarts if target fitness value was achieved or MaxFunEvals @@ -200,39 +212,39 @@ double *optimize(double (*pFun)(double const *), int nrestarts, return xbestever; // was dynamically allocated, should be deleted in the end } -double f_rand(double const *) { +double f_rand(const std::vector &) { double d = (double)rand() / RAND_MAX; while (d == 0.) d = (double)rand() / RAND_MAX; return d; } -double f_constant(double const *) { return 1; } +double f_constant(const std::vector &) { return 1; } static double SQR(double d) { return d * d; } -double f_stepsphere(double const *x) { - int i; +double f_stepsphere(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 0; i < DIM; ++i) sum += floor(x[i] * x[i]); return sum; } -double f_sphere(double const *x) { - int i; +double f_sphere(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 0; i < DIM; ++i) sum += x[i] * x[i]; return sum; } -double f_cigar(double const *x) { - int i; +double f_cigar(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 1; i < DIM; ++i) sum += x[i] * x[i]; @@ -241,10 +253,10 @@ double f_cigar(double const *x) { return sum; } -double f_cigtab(double const *x) { - int i; +double f_cigtab(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); sum = x[0] * x[0] + 1e8 * x[DIM - 1] * x[DIM - 1]; for (i = 1; i < DIM - 1; ++i) @@ -252,10 +264,10 @@ double f_cigtab(double const *x) { return sum; } -double f_tablet(double const *x) { - int i; +double f_tablet(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); sum = 1e6 * x[0] * x[0]; for (i = 1; i < DIM; ++i) @@ -264,8 +276,8 @@ double f_tablet(double const *x) { } /* a hack, memory is never released */ -double **OrthogonalBasis(int DIM) { - static int b_dim; +double **OrthogonalBasis(size_t DIM) { + static size_t b_dim; static double **b; double sp; int i, j, k; @@ -283,13 +295,13 @@ double **OrthogonalBasis(int DIM) { // Otherwise initialize basis b // allocate b - b = (double **)calloc((unsigned)DIM, sizeof(double *)); + b = (double **)calloc(DIM, sizeof(double *)); if (!b) { printf("calloc failed in function OrthogonalBasis in file example2.cpp"); exit(0); } for (i = 0; i < DIM; ++i) { - b[i] = (double *)calloc((unsigned)DIM, sizeof(double)); + b[i] = (double *)calloc(DIM, sizeof(double)); if (!b[i]) { printf("calloc failed in function Orthogonalbasis in file example2.cpp"); exit(0); @@ -319,10 +331,10 @@ double **OrthogonalBasis(int DIM) { return b; } -double f_ellirot(double const *x) { - int i, k; +double f_ellirot(const std::vector &x) { + size_t i, k; double sum = 0., y; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); double **B = OrthogonalBasis(DIM); if (DIM == 1) @@ -335,10 +347,10 @@ double f_ellirot(double const *x) { return sum; } -double f_elli(double const *x) { - int i; +double f_elli(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); if (DIM == 1) return x[0] * x[0]; @@ -347,10 +359,10 @@ double f_elli(double const *x) { return sum; } -double f_elli100(double const *x) { - int i; +double f_elli100(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); if (DIM == 1) return x[0] * x[0]; @@ -359,10 +371,10 @@ double f_elli100(double const *x) { return sum; } -double f_diffpow(double const *x) { - int i; +double f_diffpow(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); if (DIM == 1) return x[0] * x[0]; @@ -371,10 +383,10 @@ double f_diffpow(double const *x) { return sum; } -double f_diffpowrot(double const *x) { - int i, k; +double f_diffpowrot(const std::vector &x) { + size_t i, k; double sum = 0., y; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); double **B = OrthogonalBasis(DIM); if (DIM == 1) @@ -387,10 +399,10 @@ double f_diffpowrot(double const *x) { return sum; } -double f_kugelmin1(double const *x) { - int i; +double f_kugelmin1(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 1; i < DIM; ++i) sum += x[i] * x[i]; @@ -400,39 +412,39 @@ double f_kugelmin1(double const *x) { /** * Rosenbrock's Function, generalized. */ -double f_rosenbrock(double const *x) { +double f_rosenbrock(const std::vector &x) { double qualitaet; - int i; - int DIM = (int)(x[-1]); + size_t i; + size_t DIM = x.size(); qualitaet = 0.0; - for (i = DIM - 2; i >= 0; --i) - qualitaet += 100. * SQR(SQR(x[i]) - x[i + 1]) + SQR(1. - x[i]); + for (i = DIM - 1; i > 0; --i) + qualitaet += 100. * SQR(SQR(x[i - 1]) - x[i]) + SQR(1. - x[i - 1]); return (qualitaet); } -double f_parabR(double const *x) { - int i; +double f_parabR(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 1; i < DIM; ++i) sum += x[i] * x[i]; return -x[0] + 100. * sum; } -double f_sharpR(double const *x) { - int i; +double f_sharpR(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); for (i = 1; i < DIM; ++i) sum += x[i] * x[i]; return -x[0] + 100 * sqrt(sum); } -double f_ellinumtest(double const *x) { - int i; +double f_ellinumtest(const std::vector &x) { + size_t i; double sum = 0.; - int DIM = (int)(x[-1]); + size_t DIM = x.size(); static double maxVerhaeltnis = 0.; if (maxVerhaeltnis == 0.) { for (maxVerhaeltnis = 1.; @@ -460,7 +472,7 @@ double f_ellinumtest(double const *x) { * This is the reason why the square of the left side is in the quality * function. */ -double f_gleichsys5(double const *x) { +double f_gleichsys5(const std::vector &x) { double qualitaet = 0.0; static double koeff[5][6] = {/* c_1, c_2, c_3, c_4, c_5, c_0 */ @@ -469,7 +481,7 @@ double f_gleichsys5(double const *x) { {27, 1413, 191, 1032, 118, -94}, {199, 5402, 1032, 29203, 2331, 78172}, {21, 684, 118, 2331, 199, 5648}}; - int i, j; + size_t i, j; double sum; for (i = 0; i < 5; ++i) { diff --git a/test/testlexi.cpp b/test/testlexi.cpp index a5577c6..7000dbc 100644 --- a/test/testlexi.cpp +++ b/test/testlexi.cpp @@ -1,26 +1,40 @@ #include -#include #include +#include + +void check(bool b, const std::string &msg) { + if (!b) { + std::cout << msg << std::endl; + } +} int main() { - LexiProduct p1, p2, p2a, p3, p4; + LexiProduct p1, p2, dp, p2a, p3, p4; p1 = std::vector{11, 3, 4}; p2 = std::vector{11, 3, 7}; + dp = std::vector{11, 3, 7}; p2a = std::vector{11, 3, 7}; p3 = std::vector{11, 3, 7, 1}; p4 = std::vector{11, 2, 7, 1}; - assert(p2 == p2a); - assert(p1 != p2); - assert(p1 < p2); - assert(p1 <= p2); - assert(p2 > p1); - assert(p2 >= p1); - assert(p2 >= p2a); - assert(p2 <= p2a); - assert(p2 < p3); - assert(p3 > p2); - assert(p2 > p4); + check(p2 == p2a, "p2 == p2a not passed"); + check(p1 != p2, "p1 != p2 not passed"); + check(p1 < p2, "p1 < p2 not passed"); + check(p1 <= p2, "p1 <= p2 not passed"); + check(p2 > p1, "p2 > p1 not passed"); + check(p2 >= p1, "p2 >= p1 not passed"); + check(p2 >= p2a, "p2 >= p2a not passed"); + check(p2 <= p2a, "p2 <= p2a not passed"); + check(p2 < p3, "p2 < p3 not passed"); + check(p3 > p2, "p3 > p2 not passed"); + check(p2 > p4, "p2 > p4 not passed"); + check(p1 - p2 < + LexiProduct(std::vector{std::numeric_limits::max()}), + "p1 - p2 < [max] not passed"); + + LexiProduct f1, f2; + f1 = std::vector{3.5}; + f2 = std::vector{3.8}; - std::cout << "testlexi: All tests passed" << std::endl; + check(f1 < f2, "f1 < f2 not passed"); }; \ No newline at end of file