From 39e6fbbba96bd9a61a32f951b85ba67a5e0ee1ff Mon Sep 17 00:00:00 2001 From: Guangyang Wen Date: Mon, 20 Nov 2017 17:12:09 +0100 Subject: [PATCH] perceptionL --- CMakeLists.txt | 8 + cmake/Modules/FindNumPy.cmake | 55 ++++++ include/cma-es/cmaes.h | 348 +++++++++++++++++++++------------- include/color.h | 2 - include/fitness.h | 17 ++ include/lexi.h | 15 +- src/conversion.cpp | 7 - src/fitness.cpp | 131 +++++++++++++ test/example1.cpp | 16 +- test/example2.cpp | 34 ++-- test/testCIEDE2000.cpp | 27 --- test/testconversion.cpp | 75 ++++++++ test/testfitness.cpp | 62 ++++++ test/testfitness2.cpp | 19 ++ 14 files changed, 619 insertions(+), 197 deletions(-) create mode 100644 cmake/Modules/FindNumPy.cmake create mode 100644 include/fitness.h create mode 100644 src/fitness.cpp create mode 100644 test/testconversion.cpp create mode 100644 test/testfitness.cpp create mode 100644 test/testfitness2.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b75a9f0..2b26fa5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,8 @@ cmake_minimum_required (VERSION 3.1) project(color) set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wconversion -Wno-unused-private-field -Wno-unused-variable -Wno-unused-parameter -fno-omit-frame-pointer -fsanitize=address") +set(CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_STATIC_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address") 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/") @@ -16,6 +18,12 @@ target_include_directories(color PUBLIC include PRIVATE ${PYTHON_INCLUDE_DIR} ${ add_executable(testCIEDE2000 test/testCIEDE2000.cpp) target_link_libraries(testCIEDE2000 color) +add_executable(testconversion test/testconversion.cpp) +target_link_libraries(testconversion color) +add_executable(testfitness test/testfitness.cpp) +target_link_libraries(testfitness color) +add_executable(testfitness2 test/testfitness2.cpp) +target_link_libraries(testfitness2 color) add_executable(evo1 test/example1.cpp) target_link_libraries(evo1 color) diff --git a/cmake/Modules/FindNumPy.cmake b/cmake/Modules/FindNumPy.cmake new file mode 100644 index 0000000..248f8c2 --- /dev/null +++ b/cmake/Modules/FindNumPy.cmake @@ -0,0 +1,55 @@ +# - Find the NumPy libraries +# This module finds if NumPy is installed, and sets the following variables +# indicating where it is. +# +# TODO: Update to provide the libraries and paths for linking npymath lib. +# +# NUMPY_FOUND - was NumPy found +# NUMPY_VERSION - the version of NumPy found as a string +# NUMPY_VERSION_MAJOR - the major version number of NumPy +# NUMPY_VERSION_MINOR - the minor version number of NumPy +# NUMPY_VERSION_PATCH - the patch version number of NumPy +# NUMPY_VERSION_DECIMAL - e.g. version 1.6.1 is 10601 +# NUMPY_INCLUDE_DIR - path to the NumPy include files + +unset(NUMPY_VERSION) +unset(NUMPY_INCLUDE_DIR) + +if(PYTHONINTERP_FOUND) + execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-c" + "import numpy as n; print(n.__version__); print(n.get_include());" + RESULT_VARIABLE __result + OUTPUT_VARIABLE __output + OUTPUT_STRIP_TRAILING_WHITESPACE) + + if(__result MATCHES 0) + string(REGEX REPLACE ";" "\\\\;" __values ${__output}) + string(REGEX REPLACE "\r?\n" ";" __values ${__values}) + list(GET __values 0 NUMPY_VERSION) + list(GET __values 1 NUMPY_INCLUDE_DIR) + + string(REGEX MATCH "^([0-9])+\\.([0-9])+\\.([0-9])+" __ver_check "${NUMPY_VERSION}") + if(NOT "${__ver_check}" STREQUAL "") + set(NUMPY_VERSION_MAJOR ${CMAKE_MATCH_1}) + set(NUMPY_VERSION_MINOR ${CMAKE_MATCH_2}) + set(NUMPY_VERSION_PATCH ${CMAKE_MATCH_3}) + math(EXPR NUMPY_VERSION_DECIMAL + "(${NUMPY_VERSION_MAJOR} * 10000) + (${NUMPY_VERSION_MINOR} * 100) + ${NUMPY_VERSION_PATCH}") + string(REGEX REPLACE "\\\\" "/" NUMPY_INCLUDE_DIR ${NUMPY_INCLUDE_DIR}) + else() + unset(NUMPY_VERSION) + unset(NUMPY_INCLUDE_DIR) + message(STATUS "Requested NumPy version and include path, but got instead:\n${__output}\n") + endif() + endif() +else() + message(STATUS "To find NumPy Python interpretator is required to be found.") +endif() + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(NumPy REQUIRED_VARS NUMPY_INCLUDE_DIR NUMPY_VERSION + VERSION_VAR NUMPY_VERSION) + +if(NUMPY_FOUND) + message(STATUS "NumPy ver. ${NUMPY_VERSION} found (include: ${NUMPY_INCLUDE_DIR})") +endif() diff --git a/include/cma-es/cmaes.h b/include/cma-es/cmaes.h index 95ae82f..db50821 100644 --- a/include/cma-es/cmaes.h +++ b/include/cma-es/cmaes.h @@ -106,9 +106,11 @@ #include #include #include +#include #include #include #include +#include #include /** @@ -166,7 +168,7 @@ template class CMAES { }; /** - * Keys for getPtr(). + * Keys for getVector(). */ enum GetVector { NoVector = 0, @@ -247,7 +249,7 @@ template class CMAES { //! Objective function values of the population. Ordered *functionValues; //!< Public objective function value array returned by init(). - Ordered *publicFitness; + // Ordered *publicFitness; //! Generation number. size_t gen; @@ -256,7 +258,9 @@ template class CMAES { // repeatedly used for output Scalar maxdiagC; + size_t maxdiagCi; Scalar mindiagC; + size_t mindiagCi; Scalar maxEW; Scalar minEW; @@ -274,14 +278,14 @@ template class CMAES { time_t firstwritetime; time_t firstprinttime; - std::string - stopMessage; //!< A message that contains all matched stop criteria. + //! stopCriteria[].first: criterion satisfied, .second: message + std::vector> stopCriteria; - std::string getTimeStr(void) { - time_t tm = time(0); - std::string timeStr(ctime(&tm)); - return timeStr.substr(0, 24); - } + // std::string getTimeStr(void) { + // time_t tm = time(0); + // std::string timeStr(ctime(&tm)); + // return timeStr.substr(0, 24); + // } /** * Calculating eigenvalues and vectors. @@ -376,11 +380,11 @@ template class CMAES { if (tst1 < smallSDElement) tst1 = smallSDElement; const Scalar epsTst1 = eps * tst1; - int m = l; - while (m < n) { + int m = l - 1; + while (m < n - 1) { + ++m; if (std::fabs(e[m]) <= epsTst1) break; - m++; } // if m == l, d[l] is an eigenvalue, otherwise, iterate. @@ -566,13 +570,14 @@ template class CMAES { /** * Dirty index sort. + * @param col should support [] indexing */ - template - void sortIndex(const LessThanComparable *rgFunVal, size_t *iindex, size_t n) { + template + void sortIndex(const Collection &col, 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]) + if (col[iindex[j - 1]] < col[i]) break; iindex[j] = iindex[j - 1]; // shift up } @@ -613,12 +618,15 @@ template class CMAES { } // update maximal and minimal diagonal value maxdiagC = mindiagC = C[0][0]; - for (int i = 1; i < N; ++i) { + for (size_t i = 1; i < N; ++i) { const Scalar &Cii = C[i][i]; - if (maxdiagC < Cii) + if (maxdiagC < Cii) { maxdiagC = Cii; - else if (mindiagC > Cii) + maxdiagCi = i; + } else if (mindiagC > Cii) { mindiagC = Cii; + mindiagCi = i; + } } } } @@ -841,26 +849,24 @@ template class CMAES { } // for (int i = 0; i < params.lambda; ++i) // delete[]-- population[i]; - delete[] population; delete[] C; delete[] B; delete[] index; - delete[] publicFitness; + // delete[] publicFitness; + delete[] population; delete[] functionValues; delete[] funcValueHistory; } /** - * Initializes the CMA-ES algorithm. + * Allocate ressources and initializes the CMA-ES algorithm. * @param parameters The CMA-ES parameters. - * @return Array of size lambda that can be used to assign fitness values and - * pass them to updateDistribution(). Not that after the desctructor - * was called, the array is deleted. */ - Ordered *init(const Parameters ¶meters) { + void init(const Parameters ¶meters) { params = parameters; - stopMessage = ""; + stopCriteria.resize(10); // There are ten stop criteria + // stopMessage = ""; Scalar trace(0); for (int i = 0; i < params.N; ++i) @@ -881,8 +887,8 @@ template class CMAES { if (dtest == dtest + Scalar(1)) break; dMaxSignifKond = - dtest / Scalar(1000); // not sure whether this is really save, - // 100 does not work well enough + dtest / Scalar(1e2); // not sure whether this is really save, + // 100 does not work well enough gen = 0; countevals = 0; @@ -912,7 +918,7 @@ template class CMAES { rgD = new Scalar[params.N]; C = new Scalar *[params.N]; B = new Scalar *[params.N]; - publicFitness = new Ordered[params.lambda]; + // publicFitness = new Ordered[params.lambda]; functionValues = new Ordered[params.lambda]; // functionValues[0] = params.lambda; // ++functionValues; @@ -982,20 +988,20 @@ template class CMAES { if (params.resumefile != "") resumeDistribution(params.resumefile); - return publicFitness; + // return publicFitness; } /** * Well, says hello. - * @return eg. "(5,10)-CMA-ES(mu_eff=3.4), Ver="1.0alpha", dimension=9" + * @return eg. "CMA-ES-1.1alpha lambda=263, mu=131, mu_eff=68.6358, + * dimension=14, diagonalIterations=0" */ std::string sayHello() { std::stringstream stream; - stream << "(" << params.mu << "," << params.lambda - << ")-CMA-ES(mu_eff=" << std::setprecision(1) << params.mueff - << "), Ver=\"" << version << "\", dimension=" << params.N - << ", diagonalIterations=" << (long)params.diagonalCov << " (" - << getTimeStr() << ")"; + stream << "CMA-ES-" << version << " lambda=" << params.lambda + << ", mu=" << params.mu << ", mu_eff=" << params.mueff + << ", dimension=" << params.N + << ", diagonalIterations=" << (long)params.diagonalCov; return stream.str(); } @@ -1135,8 +1141,6 @@ template class CMAES { } /** - * The search space vectors have a special form: they are arrays with N+1 - * entries. Entry number -1 is the dimension of the search space N. * @return A pointer to a "population" of lambda N-dimensional multivariate * normally distributed samples. */ @@ -1262,18 +1266,19 @@ template class CMAES { * Core procedure of the CMA-ES algorithm. Sets a new mean value and estimates * the new covariance matrix and a new step size for the normal search * distribution. - * @param fitnessValues An array of \f$\lambda\f$ function values. + * @param fitnessValues An iterable of \f$\lambda\f$ function values. * @return Mean value of the new distribution. */ - Scalar *updateDistribution(const Ordered *fitnessValues) { + template + Scalar *updateDistribution(const OrderedCollection &fitnessValues) { const int N = params.N; bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; assert(state != UPDATED && "updateDistribution(): You need to call " "samplePopulation() before update can take place."); - assert(fitnessValues && - "updateDistribution(): No fitness function value array input."); + // assert(fitnessValues && + // "updateDistribution(): No fitness function value array input."); if (state == SAMPLED) // function values are delivered here countevals += params.lambda; @@ -1281,7 +1286,7 @@ template class CMAES { params.logStream << "updateDistribution(): unexpected state" << std::endl; // assign function values - for (int i = 0; i < params.lambda; ++i) { + for (size_t i = 0; i < params.lambda; ++i) { // std::cout << i << ": " << fitnessValues[i] << std::endl; population[i].value = functionValues[i] = fitnessValues[i]; } @@ -1321,6 +1326,7 @@ template class CMAES { // std::cout << xBestEver.value << " -> " << population[index[0]].value << // std::endl; xBestEver = population[index[0]]; + // std::cout << xBestEver.value << std::endl; evalsBestEver = countevals; // for (int i = 0; i < N; ++i) { // // xBestEver.x[i] = population[index[0]].x[i]; /* TODO: use @@ -1458,7 +1464,7 @@ template class CMAES { * overwritten during the next call to any member functions other * than get(). */ - const std::vector getPtr(GetVector key) { + const std::vector getVector(GetVector key) { switch (key) { case DiagC: { std::vector o((size_t)params.N); @@ -1481,7 +1487,7 @@ template class CMAES { case XMean: return std::vector(xmean, xmean + params.N); default: - throw std::runtime_error("getPtr(): No match found for key"); + throw std::runtime_error("getVector(): No match found for key"); } } @@ -1492,7 +1498,7 @@ template class CMAES { * writing access to its elements. The memory must be explicitly * released using delete[]. */ - Scalar *getNew(GetVector key) { return getInto(key, 0); } + // Scalar *getNew(GetVector key) { return getInto(key, nullptr); } /** * Request a vector parameter from CMA-ES. @@ -1501,14 +1507,14 @@ 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. */ - Scalar *getInto(GetVector key, Scalar *res) { - const std::vector res0 = getPtr(key); - if (!res) - res = new Scalar[params.N]; - for (size_t i = 0; i < params.N; ++i) - res[i] = res0[i]; - return res; - } + // Scalar *getInto(GetVector key, Scalar *res) { + // const std::vector res0 = getVector(key); + // if (!res) + // res = new Scalar[params.N]; + // for (size_t i = 0; i < params.N; ++i) + // res[i] = res0[i]; + // return res; + // } /** * Some stopping criteria can be set in initials.par, with names starting @@ -1518,125 +1524,205 @@ template class CMAES { * that contains the matched stop criteria via getStopMessage(). * @return Does any stop criterion match? */ - bool testForTermination() { + bool testForTermination(bool makeMessage = false) { Ordered range; Scalar fac; int iAchse, iKoo; int diag = params.diagonalCov == 1 || params.diagonalCov >= gen; int N = params.N; - std::stringstream message; + // std::stringstream message; - if (stopMessage != "") { - message << stopMessage << std::endl; - } + // if (stopMessage != "") { + // message << stopMessage << std::endl; + // } - // function value reached - if ((gen > 1 || state > SAMPLED) && params.stStopFitness.flg && - functionValues[index[0]] <= params.stStopFitness.val) { - message << "Fitness: function value " << functionValues[index[0]] - << " <= stopFitness (" << params.stStopFitness.val << ")" - << std::endl; + // 0. function value reached + { + stopCriteria[0].first = + (gen > 1 || state > SAMPLED) && params.stStopFitness.flg && + functionValues[index[0]] <= params.stStopFitness.val; + if (stopCriteria[0].first || makeMessage) { + std::stringstream message; + message << "Fitness: function value " << functionValues[index[0]] + << " <= stopFitness (" << params.stStopFitness.val << ")"; + stopCriteria[0].second = message.str(); + } } - // TolFun - range = std::max(maxElement(funcValueHistory, - std::min((size_t)gen, historySize)), - maxElement(functionValues, (size_t)params.lambda)) - - std::min(minElement(funcValueHistory, - std::min((size_t)gen, historySize)), - minElement(functionValues, (size_t)params.lambda)); - - if (gen > 0 && range <= params.stopTolFun) { - message << "TolFun: function value differences " << range - << " < stopTolFun=" << params.stopTolFun << std::endl; + // 1. TolFun + { + range = std::max(maxElement(funcValueHistory, + std::min((size_t)gen, historySize)), + maxElement(functionValues, (size_t)params.lambda)) - + std::min(minElement(funcValueHistory, + std::min((size_t)gen, historySize)), + minElement(functionValues, (size_t)params.lambda)); + + stopCriteria[1].first = gen > 0 && range <= params.stopTolFun; + if (stopCriteria[1].first || makeMessage) { + std::stringstream message; + message << "TolFun: function value differences " << range + << " < stopTolFun=" << params.stopTolFun; + stopCriteria[1].second = message.str(); + } } - // TolFunHist - if (gen > historySize) { - range = maxElement(funcValueHistory, historySize) - - minElement(funcValueHistory, historySize); - if (range <= params.stopTolFunHist) + // 2. TolFunHist + { + if (gen > historySize) { + range = maxElement(funcValueHistory, historySize) - + minElement(funcValueHistory, historySize); + if (range <= params.stopTolFunHist) + stopCriteria[2].first = true; + } else { + stopCriteria[2].first = false; + } + if (stopCriteria[2].first || makeMessage) { + std::stringstream message; message << "TolFunHist: history of function value changes " << range - << " stopTolFunHist=" << params.stopTolFunHist << std::endl; + << " stopTolFunHist=" << params.stopTolFunHist; + stopCriteria[2].second = message.str(); + } } - // TolX - int cTemp = 0; - for (int i = 0; i < N; ++i) { - cTemp += (sigma * std::sqrt(C[i][i]) < params.stopTolX) ? 1 : 0; - cTemp += (sigma * pc[i] < params.stopTolX) ? 1 : 0; - } - if (cTemp == 2 * N) { - message << "TolX: object variable changes below " << params.stopTolX - << std::endl; + // 3. TolX + { + int cTemp = 0; + for (int i = 0; i < N; ++i) { + cTemp += (sigma * std::sqrt(C[i][i]) < params.stopTolX) ? 1 : 0; + cTemp += (sigma * pc[i] < params.stopTolX) ? 1 : 0; + } + stopCriteria[3].first = cTemp == 2 * N; + if (stopCriteria[3].first || makeMessage) { + std::stringstream message; + message << "TolX: object variable changes below " << params.stopTolX; + stopCriteria[3].second = message.str(); + } } - // TolUpX - for (int i = 0; i < N; ++i) { - if (sigma * std::sqrt(C[i][i]) > - params.stopTolUpXFactor * params.rgInitialStds[i]) { + // 4. TolUpX + { + stopCriteria[4].first = false; + for (int i = 0; i < N; ++i) { + if (sigma * std::sqrt(C[i][i]) > + params.stopTolUpXFactor * params.rgInitialStds[i]) { + stopCriteria[4].first = true; + break; + } + } + if (stopCriteria[4].first || makeMessage) { + std::stringstream message; message << "TolUpX: standard deviation increased by more than " << params.stopTolUpXFactor - << ", larger initial standard deviation recommended." - << std::endl; - break; + << ", larger initial standard deviation recommended."; + stopCriteria[4].second = message.str(); } } - // Condition of C greater than dMaxSignifKond - if (maxEW >= minEW * dMaxSignifKond) { - message << "ConditionNumber: maximal condition number " << dMaxSignifKond - << " reached. maxEW=" << maxEW << ",minEW=" << minEW - << ",maxdiagC=" << maxdiagC << ",mindiagC=" << mindiagC - << std::endl; + // 5. Condition of C greater than dMaxSignifKond + { + stopCriteria[5].first = maxEW >= minEW * dMaxSignifKond; + if (stopCriteria[5].first || makeMessage) { + std::stringstream message; + message << "ConditionNumber: maximal condition number " + << dMaxSignifKond << " reached. maxEW=" << maxEW + << ",minEW=" << minEW << ",maxdiagC=" << maxdiagC + << ",maxdiagCi=" << maxdiagCi << ",mindiagC=" << mindiagC + << ",mindiagCi=" << mindiagCi; + stopCriteria[5].second = message.str(); + } } - // Principal axis i has no effect on xmean, ie. x == x + 0.1* sigma* rgD[i]* - // B[i] - if (!diag) { - for (iAchse = 0; iAchse < N; ++iAchse) { - fac = 0.1 * sigma * rgD[iAchse]; - for (iKoo = 0; iKoo < N; ++iKoo) { - if (xmean[iKoo] != xmean[iKoo] + fac * B[iKoo][iAchse]) + // 6. Principal axis i has no effect on xmean, ie. x == x + 0.1* sigma* + // rgD[i]* B[i] + { + if (!diag) { + std::vector> no_effect; + for (iAchse = 0; iAchse < N; ++iAchse) { + fac = 0.1 * sigma * rgD[iAchse]; + for (iKoo = 0; iKoo < N; ++iKoo) { + if (xmean[iKoo] != xmean[iKoo] + fac * B[iKoo][iAchse]) + break; + } + if (iKoo == N) { + stopCriteria[6].first = true; + no_effect.emplace_back(fac / 0.1, iAchse); break; + } } - if (iKoo == N) { - message << "NoEffectAxis: standard deviation 0.1*" << (fac / 0.1) - << " in principal axis " << iAchse << " without effect" - << std::endl; - break; + if (stopCriteria[6].first || makeMessage) { + std::stringstream message; + message << "NoEffectAxis: "; + for (const auto &p : no_effect) + message << "[axis " << p.second << " std dev 0.1*" << p.first + << "] "; + stopCriteria[6].second = message.str(); } } } - // Component of xmean is not changed anymore - for (iKoo = 0; iKoo < N; ++iKoo) { - if (xmean[iKoo] == - xmean[iKoo] + sigma * std::sqrt(C[iKoo][iKoo]) / Scalar(5)) { + + // 7. Component of xmean is not changed anymore + { + for (iKoo = 0; iKoo < N; ++iKoo) { + if (xmean[iKoo] == + xmean[iKoo] + sigma * std::sqrt(C[iKoo][iKoo]) / Scalar(5)) { + stopCriteria[7].first = true; + break; + } + } + if (stopCriteria[7].first) { + std::stringstream message; message << "NoEffectCoordinate: standard deviation 0.2*" << (sigma * std::sqrt(C[iKoo][iKoo])) << " in coordinate " - << iKoo << " without effect" << std::endl; - break; + << iKoo << " without effect"; + stopCriteria[7].second = message.str(); + } else if (makeMessage) { + std::stringstream message; + message << "NoEffectCoordinate: (empty)"; + stopCriteria[7].second = message.str(); } } - if (countevals >= params.stopMaxFunEvals) { - message << "MaxFunEvals: conducted function evaluations " << countevals - << " >= " << params.stopMaxFunEvals << std::endl; + // 8. Maximum function evaluations counts is reached + { + stopCriteria[8].first = countevals >= params.stopMaxFunEvals; + if (stopCriteria[8].first || makeMessage) { + std::stringstream message; + message << "MaxFunEvals: conducted function evaluations " << countevals + << " >= " << params.stopMaxFunEvals; + stopCriteria[8].second = message.str(); + } } - if (gen >= params.stopMaxIter) { - message << "MaxIter: number of iterations " << gen - << " >= " << params.stopMaxIter << std::endl; + + // 9. Maximum generation counts is reached + { + stopCriteria[9].first = gen >= params.stopMaxIter; + if (stopCriteria[9].first || makeMessage) { + std::stringstream message; + message << "MaxIter: number of iterations " << gen + << " >= " << params.stopMaxIter; + stopCriteria[9].second = message.str(); + } } - stopMessage = message.str(); - return stopMessage != ""; + return std::accumulate( + stopCriteria.begin(), stopCriteria.end(), false, + [](bool accu, const std::pair crit) { + return accu || crit.first; + }); } /** * A message that contains a detailed description of the matched stop * criteria. */ - std::string getStopMessage() { return stopMessage; } + std::string getStopMessage() { + testForTermination(true); /* generate all messages */ + std::stringstream message; + for (const std::pair &p : stopCriteria) + message << "[" << (p.first ? "x" : " ") << "] " << p.second << std::endl; + return message.str(); + } /** * @param filename Output file name. diff --git a/include/color.h b/include/color.h index 4528a1a..10fd456 100644 --- a/include/color.h +++ b/include/color.h @@ -60,8 +60,6 @@ RGB XYZtoRGB(const XYZ &xyz); XYZ LABtoXYZ(const LAB &lab); RGB LABtoRGB(const LAB &lab); -bool offRGB(const LAB &lab); - } // namespace color std::ostream &operator<<(std::ostream &s, const color::LAB &labColor); diff --git a/include/fitness.h b/include/fitness.h new file mode 100644 index 0000000..66b24a1 --- /dev/null +++ b/include/fitness.h @@ -0,0 +1,17 @@ +#ifndef FITNESS_H_ +#define FITNESS_H_ + +#include "color.h" +#include "lexi.h" +#include + +inline double offRange(double x, double a, double b); + +double offRGB(const color::LAB &lab); + +LexiProduct fitnessFunc(double L, const std::vector &x); + +std::pair, LexiProduct> +perceptionL(double L, size_t M, bool quiet = false); + +#endif \ No newline at end of file diff --git a/include/lexi.h b/include/lexi.h index 695e7c8..ea7e2e8 100644 --- a/include/lexi.h +++ b/include/lexi.h @@ -1,3 +1,6 @@ +#ifndef LEXI_PRODUCT_H_ +#define LEXI_PRODUCT_H_ + #include #include #include @@ -57,11 +60,15 @@ template class LexiProduct { return LexiProduct(std::move(d)); } + Scalar const &operator[](size_t i) const { return this->prd[i]; } + + Scalar &operator[](size_t i) { return this->prd[i]; } + private: Product prd; - Scalar lexicmp(const LexiProduct &lhs, - const LexiProduct &rhs) const { + inline 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) { @@ -82,4 +89,6 @@ template class LexiProduct { os << lexi.prd[n - 1] << "]"; return os; } -}; \ No newline at end of file +}; + +#endif \ No newline at end of file diff --git a/src/conversion.cpp b/src/conversion.cpp index 03d94be..13afb09 100644 --- a/src/conversion.cpp +++ b/src/conversion.cpp @@ -42,11 +42,4 @@ XYZ LABtoXYZ(const LAB &lab) { RGB LABtoRGB(const LAB &lab) { return XYZtoRGB(LABtoXYZ(lab)); } -bool offRGB(const LAB &lab) { - const RGB rgb{LABtoRGB(lab)}; - - return rgb.r >= 0 && rgb.r <= 1 && rgb.g >= 0 && rgb.g <= 1 && rgb.b >= 0 && - rgb.b <= 1; -} - } // namespace color diff --git a/src/fitness.cpp b/src/fitness.cpp new file mode 100644 index 0000000..b29843a --- /dev/null +++ b/src/fitness.cpp @@ -0,0 +1,131 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +inline double offRange(double x, double a, double b) { + if (x < a) { + return a - x; + } + if (x > b) { + return x - b; + } + return 0; +} + +double offRGB(const color::LAB &lab) { + const color::RGB rgb{color::LABtoRGB(lab)}; + const auto dr = offRange(rgb.r, 0, 1); + const auto dg = offRange(rgb.g, 0, 1); + const auto db = offRange(rgb.b, 0, 1); + + return std::sqrt(dr * dr + dg * dg + db * db); +} + +LexiProduct fitnessFunc(double L, const std::vector &x) { + // number of colors + size_t M = x.size() / 2; + std::vector lab(M); + // off-RGB penalty + std::vector penalty(M); + for (size_t i = 0; i < M; ++i) { + lab[i] = color::LAB{L, x[2 * i], x[2 * i + 1]}; + penalty[i] = + offRGB(lab[i]) * 300 * M; /* 300 is the approx. diameter of ab-plane */ + } + + size_t K = M * (M - 1) / 2; /* number of combinations between colors */ + static size_t cached_K = 0; + static std::vector cached_ref; /* for almost-sorted O(n) sorting */ + static std::vector> ref_combi; + + if (cached_K != K) { + cached_K = K; + cached_ref.resize(K); + for (size_t i = 0; i < K; ++i) { + cached_ref[i] = i; + } + ref_combi.resize(K); + size_t combi_i = 0; + for (size_t i = 0; i < M - 1; ++i) { + for (size_t j = i + 1; j < M; ++j) { + ref_combi[combi_i] = std::pair(i, j); + ++combi_i; + } + } + } + + std::vector> fitness_ref(K); + for (size_t i = 0; i < K; ++i) { + const auto ref = cached_ref[i]; + const auto &refs = ref_combi[ref]; + fitness_ref[i] = std::pair( + penalty[refs.first] + penalty[refs.second] - + color::CIEDE2000(lab[refs.first], lab[refs.second]), + ref); + } + std::sort( + fitness_ref.begin(), fitness_ref.end(), + [](const std::pair &a, + const std::pair &b) { return a.first > b.first; }); + std::vector fitness(K); + for (size_t i = 0; i < K; ++i) { + const auto &pair = fitness_ref[i]; + fitness[i] = pair.first; + cached_ref[i] = pair.second; + } + return LexiProduct(std::move(fitness)); +} + +/* + * @param L luminosity + * @param M numbers of colors + * @param quiet write info to stdout + * @return pair(RGBs, bestFitness) + */ +std::pair, LexiProduct> +perceptionL(double L, size_t M, bool quiet) { + CMAES> evo; + Individual> *pop; + + const size_t N = M * 2; //!< number of variables + const size_t K = M * (M - 1) / 2; //!< number of combinations + Parameters> parameters; + parameters.lambda = (int)(100. * log(N)); + // parameters.diagonalCov = 1; + // parameters.stopTolFun = 1e-12; + // parameters.stopTolFunHist = 1e-13; + const std::vector xstart(N, 0.), stddev(N, std::min(100. - L, L)); + parameters.init((int)N, xstart.data(), stddev.data()); + evo.init(parameters); + + if (!quiet) + std::cout << evo.sayHello() << std::endl; + + std::vector> populationFitness((size_t)parameters.lambda); + while (!evo.testForTermination()) { + pop = evo.samplePopulation(); // Do not change content of pop + for (size_t i = 0; i < evo.get(CMAES>::Lambda); + ++i) + populationFitness[i] = fitnessFunc(50., pop[i].x); + evo.updateDistribution(populationFitness); + } + + if (!quiet) + std::cout << "Stop:" << std::endl << evo.getStopMessage(); + + auto x = evo.getVector(CMAES>::XMean); + auto fitness = evo.getValue(CMAES>::FBestEver); + + std::vector rgb(M); + for (size_t i = 0; i < M; ++i) { + rgb[i] = color::LABtoRGB(color::LAB{50., x[2 * i], x[2 * i + 1]}); + } + + return std::pair, LexiProduct>( + std::move(rgb), std::move(fitness)); +} \ No newline at end of file diff --git a/test/example1.cpp b/test/example1.cpp index 9c91442..bb23ab5 100644 --- a/test/example1.cpp +++ b/test/example1.cpp @@ -28,9 +28,7 @@ Y fitfun(const std::vector &x, int N) { */ int main(int, char **) { CMAES evo; - Y *arFunvals; Individual *pop; - double *xfinal; // Initialize everything const int dim = 22; @@ -45,10 +43,11 @@ int main(int, char **) { parameters.stopTolFunHist = std::vector{0, 1e-13}; // TODO Adjust parameters here parameters.init(dim, xstart, stddev); - arFunvals = evo.init(parameters); + evo.init(parameters); std::cout << evo.sayHello() << std::endl; + std::vector arFunvals((size_t)parameters.lambda); // Iterate until stop criterion holds while (!evo.testForTermination()) { // Generate lambda new search points, sample population @@ -69,7 +68,7 @@ int main(int, char **) { */ // evaluate the new search points using fitfun from above - for (int i = 0; i < evo.get(CMAES::Lambda); ++i) + for (size_t i = 0; i < evo.get(CMAES::Lambda); ++i) arFunvals[i] = fitfun(pop[i].x, (int)evo.get(CMAES::Dimension)); @@ -85,13 +84,12 @@ int main(int, char **) { "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 + const auto xfinal = evo.getVector(CMAES::XMean); // "XBestEver" + // might be + // used as + // well // do something with final solution and finally release memory - delete[] xfinal; return 0; } diff --git a/test/example2.cpp b/test/example2.cpp index a9a7308..34679b6 100644 --- a/test/example2.cpp +++ b/test/example2.cpp @@ -32,9 +32,9 @@ 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); +std::vector optimize(double (*pFun)(const std::vector &), + int number_of_restarts, + double increment_factor_for_population_size); int main(int, char **) { typedef double (*pfun_t)(const std::vector &); @@ -42,7 +42,6 @@ int main(int, char **) { std::vector rgpFun( maxTest); // vector (range) of pointer to objective function double incpopsize = 2; - double *x; // Put together objective functions rgpFun[0] = f_sphere; @@ -71,12 +70,9 @@ int main(int, char **) { if (f == nullptr) { continue; } - x = optimize(f, nbrestarts, incpopsize); + const auto x = optimize(f, nbrestarts, incpopsize); } - // here we could utilize the solution x, and finally free memory - delete[] x; - return 0; } @@ -84,12 +80,12 @@ int main(int, char **) { * Somewhat extended interface for optimizing pFun with CMAES implementing a * restart procedure with increasing population size. */ -double *optimize(double (*pFun)(const std::vector &), int nrestarts, - double incpopsize) { +std::vector optimize(double (*pFun)(const std::vector &), + int nrestarts, double incpopsize) { 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 fbestever = 0; + std::vector xbestever; // store best solution double fmean; int irun, lambda = 0, // offspring population size, 0 invokes default @@ -117,10 +113,11 @@ double *optimize(double (*pFun)(const std::vector &), int nrestarts, parameters.lambda = lambda; parameters.init(dim, xstart, stddev); - fitvals = evo.init(parameters); // allocs fitvals + evo.init(parameters); // allocs fitvals std::cout << evo.sayHello() << std::endl; evo.countevals = countevals; // a hack, effects the output and termination + std::vector fitvals((size_t)parameters.lambda); while (!evo.testForTermination()) { // Generate population of new candidate solutions pop = evo.samplePopulation(); // do not change content of pop @@ -135,7 +132,7 @@ double *optimize(double (*pFun)(const std::vector &), int nrestarts, */ // Compute fitness value for each candidate solution - for (int i = 0; i < evo.get(CMAES::PopSize); ++i) { + for (size_t 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 @@ -186,14 +183,15 @@ double *optimize(double (*pFun)(const std::vector &), int nrestarts, if (irun == 0 || evo.getValue(CMAES::FBestEver) < fbestever) { fbestever = evo.getValue(CMAES::FBestEver); - xbestever = evo.getInto(CMAES::XBestEver, - xbestever); // allocates memory if needed + xbestever = evo.getVector( + CMAES::XBestEver); // allocates memory if needed } // best estimator for the optimum is xmean, therefore check - const std::vector xmean = evo.getPtr(CMAES::XMean); + const std::vector xmean = + evo.getVector(CMAES::XMean); if ((fmean = (*pFun)(xmean)) < fbestever) { fbestever = fmean; - xbestever = evo.getInto(CMAES::XMean, xbestever); + xbestever = evo.getVector(CMAES::XMean); } // abandon restarts if target fitness value was achieved or MaxFunEvals diff --git a/test/testCIEDE2000.cpp b/test/testCIEDE2000.cpp index ac36c91..8f80bb7 100644 --- a/test/testCIEDE2000.cpp +++ b/test/testCIEDE2000.cpp @@ -24,35 +24,8 @@ int main(int argc, char *argv[]) { return (testcolor()); } int testcolor() { color::LAB lab1, lab2; - color::RGB rgb1, rgb_expected; unsigned int line = 1; std::vector passFail; - - lab1 = {100, 0, 0}; - rgb1 = color::LABtoRGB(lab1); - rgb_expected = {1, 1, 1}; - - passFail.push_back(std::abs(rgb1.r - rgb_expected.r) < 1e-3); - std::cout << line << ": " << rgb1.r << " vs " << rgb_expected.r << '\t' - << (passFail[line - 1] ? "PASS" : "FAIL") << std::endl; - ++line; - - passFail.push_back(std::abs(rgb1.g - rgb_expected.g) < 1e-3); - std::cout << line << ": " << rgb1.g << " vs " << rgb_expected.g << '\t' - << (passFail[line - 1] ? "PASS" : "FAIL") << std::endl; - ++line; - - passFail.push_back(std::abs(rgb1.b - rgb_expected.b) < 1e-3); - std::cout << line << ": " << rgb1.b << " vs " << rgb_expected.b << '\t' - << (passFail[line - 1] ? "PASS" : "FAIL") << std::endl; - ++line; - - std::cout << color::LABtoRGB(color::LAB{50, 0, 0}) << std::endl; - std::cout << color::LABtoRGB(color::LAB{50, 1000, 0}) << std::endl; - std::cout << color::LABtoRGB(color::LAB{50, -1000, 0}) << std::endl; - std::cout << color::LABtoRGB(color::LAB{50, 0, 1000}) << std::endl; - std::cout << color::LABtoRGB(color::LAB{50, 0, -1000}) << std::endl; - double expectedResult, myResult; char myBuffer[64], expectedBuffer[64]; std::string myStr, expectedStr; diff --git a/test/testconversion.cpp b/test/testconversion.cpp new file mode 100644 index 0000000..711b008 --- /dev/null +++ b/test/testconversion.cpp @@ -0,0 +1,75 @@ +/* + * testcolor.cpp + * Part of http://github.com/gfiumara/color by Gregory Fiumara. + * See LICENSE for details. + */ + +#include +#include +#include +#include + +#include + +static size_t line = 0; +static std::vector passFail; + +void check(bool res, std::string msg) { + passFail.push_back(res); + std::cout << line << ": " << msg << '\t' << (passFail[line] ? "PASS" : "FAIL") + << std::endl; + ++line; +} + +/** + * @brief + * Run the test dataset run from the color paper. + * + * @return + * EXIT_SUCCESS if all tests pass, EXIT_FAILURE otherwise. + */ +int testcolor(); + +int main(int argc, char *argv[]) { return (testcolor()); } + +int testcolor() { + unsigned int line = 1; + std::vector passFail; + + color::LAB lab1 = {100, 0, 0}; + color::RGB rgb1 = color::LABtoRGB(lab1); + color::RGB rgb_expected = {1, 1, 1}; + + check(std::abs(rgb1.r - rgb_expected.r) < 1e-3, "LAB->RGB"); + check(std::abs(rgb1.g - rgb_expected.g) < 1e-3, "LAB->RGB"); + check(std::abs(rgb1.b - rgb_expected.b) < 1e-3, "LAB->RGB"); + + // std::cout << color::LABtoRGB(color::LAB{50, 0, 0}) << std::endl; + + color::RGB rgb2 = color::LABtoRGB(color::LAB{50, 1000, 0}); + check(rgb2.r > 1, "RGB limit"); + check(rgb2.g < 0, "RGB limit"); + check(rgb2.b < 0, "RGB limit"); + + color::RGB rgb3 = color::LABtoRGB(color::LAB{50, -1000, 0}); + check(rgb3.r < 0, "RGB limit"); + + color::RGB rgb4 = color::LABtoRGB(color::LAB{50, 0, 1000}); + check(rgb4.b < 0, "RGB limit"); + + color::RGB rgb5 = color::LABtoRGB(color::LAB{50, 0, -1000}); + check(rgb5.r < 0, "RGB limit"); + check(rgb5.g < 0, "RGB limit"); + check(rgb5.b > 1, "RGB limit"); + + std::cout << std::endl; + int ret = EXIT_SUCCESS; + for (unsigned int i = 0; i < passFail.size(); i++) { + if (!passFail[i]) { + std::cout << "Test failed on line " << std::to_string(i + 1) << std::endl; + ret = EXIT_FAILURE; + } + } + + return (ret); +} diff --git a/test/testfitness.cpp b/test/testfitness.cpp new file mode 100644 index 0000000..28cc519 --- /dev/null +++ b/test/testfitness.cpp @@ -0,0 +1,62 @@ +/* + * testcolor.cpp + * Part of http://github.com/gfiumara/color by Gregory Fiumara. + * See LICENSE for details. + */ + +#include +#include +#include +#include + +#include +#include + +static size_t line = 0; +static std::vector passFail; + +void check(bool res, const std::string &msg) { + passFail.push_back(res); + std::cout << line << ": " << msg << '\t' << (passFail[line] ? "PASS" : "FAIL") + << std::endl; + ++line; +} + +void check(double a, double b, const std::string &msg, double delta = 1e-3) { + check(std::abs(a - b) < delta, msg); +} + +/** + * @brief + * Run the test dataset run from the color paper. + * + * @return + * EXIT_SUCCESS if all tests pass, EXIT_FAILURE otherwise. + */ +int testcolor(); + +int main(int argc, char *argv[]) { return (testcolor()); } + +int testcolor() { + unsigned int line = 1; + std::vector passFail; + + check(offRGB(color::LAB{200, 0, 0}), 2.147648383, "offRGB"); + + check(fitnessFunc(50, std::vector{0, 0, 10, 10, -10, 10, -10, -10, 10, -10})[0], + -12.8001, "fitnessFunc"); + + check(fitnessFunc(50, std::vector{-200, -200, +200, +200})[0], + 1626.445599, "fitnessFunc (off boundary)"); + + std::cout << std::endl; + int ret = EXIT_SUCCESS; + for (unsigned int i = 0; i < passFail.size(); i++) { + if (!passFail[i]) { + std::cout << "Test failed on line " << std::to_string(i + 1) << std::endl; + ret = EXIT_FAILURE; + } + } + + return (ret); +} diff --git a/test/testfitness2.cpp b/test/testfitness2.cpp new file mode 100644 index 0000000..67a776b --- /dev/null +++ b/test/testfitness2.cpp @@ -0,0 +1,19 @@ +#include +#include +#include +#include + +#include +#include +#include +#include + +int main(int argc, char *argv[]) { + const auto pair = perceptionL(50., 7); + + for (size_t i = 0; i < pair.first.size(); ++i) { + std::cout << pair.first[i] << " "; + } + std::cout << std::endl << pair.second << std::endl; + return 0; +} \ No newline at end of file