diff --git a/CMakeLists.txt b/CMakeLists.txt index 48e300e..98a4009 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,14 +2,24 @@ cmake_minimum_required (VERSION 3.1) project(color) -# Compiler-specific C++11 activation. set(CMAKE_CXX_STANDARD 11) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") + +find_package(PythonLibs REQUIRED 3) +find_package(PythonInterp REQUIRED 3) +find_package(NumPy REQUIRED) -find_library(PythonLibs 3) - -add_library(color SHARED src/debug.cpp src/conversion.cpp src/CIEDE2000.cpp) -target_include_directories(color PUBLIC include PRIVATE ${PYTHON_INCLUDE_DIR}) +add_library(color SHARED src/debug.cpp src/conversion.cpp src/CIEDE2000.cpp src/fitness.cpp) +target_include_directories(color PUBLIC include PRIVATE ${PYTHON_INCLUDE_DIR} ${NUMPY_INCLUDE_DIR}) add_executable(testCIEDE2000 test/testCIEDE2000.cpp) target_link_libraries(testCIEDE2000 color) + +add_executable(evo1 test/example1.cpp) +target_link_libraries(evo1 color) +add_executable(evo2 test/example2.cpp) +target_link_libraries(evo2 color) + +add_executable(testlexi test/testlexi.cpp) +target_link_libraries(testlexi color) \ No newline at end of file diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..94a9ed0 --- /dev/null +++ b/COPYING @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/include/cma-es/cmaes.h b/include/cma-es/cmaes.h new file mode 100644 index 0000000..98d89ec --- /dev/null +++ b/include/cma-es/cmaes.h @@ -0,0 +1,1658 @@ +/** + * @file cmaes.h + * @author Nikolaus Hansen, ported to C++ by Alexander Fabisch + * + * \mainpage + * CMA-ES for non-linear function minimization. + * + * Copyright of C implementation by Nikolaus Hansen (e-mail: + * hansen .AT. bionik.tu-berlin.de, hansen .AT. lri.fr), ported to C++ by + * Alexander Fabisch. + * + * \section lgpl License + * + * Copyright 1996, 2003, 2007, 2011 Nikolaus Hansen, Alexander Fabisch + * + * This file is part of CMA-ESpp. + * + * CMA-ESpp is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CMA-ESpp is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CMA-ESpp. If not, see . + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + * \section purpose General Purpose + * + * The CMA-ES (Evolution Strategy with Covariance Matrix Adaptation) is a + * robust search/optimization method. The goal is to minimize a given + * objective function, \f$f: R^N \rightarrow R\f$. The CMA-ES should be + * applied, if e.g. BFGS and/or conjugate gradient methods fail due to a + * rugged search landscape (e.g. discontinuities, outliers, noise, local + * optima, etc.). Learning the covariance matrix in the CMA-ES is similar + * to learning the inverse Hessian matrix in a quasi-Newton method. On + * smooth landscapes the CMA-ES is roughly ten times slower than BFGS, + * assuming derivatives are not directly available. For up to \f$N=10\f$ + * parameters the simplex direct search method (Nelder & Mead) is + * sometimes faster, but less robust than CMA-ES. On considerably hard + * problems the search (a single run) is expected to take between + * \f$100\cdot N\f$ and \f$300\cdot N^2\f$ function evaluations. But you + * might be lucky... + * + * \section application Application Remark + * + * The adaptation of the covariance matrix (e.g. by the CMA) is + * equivalent to a general linear transformation of the problem + * variables. Nevertheless, every problem specific knowledge about the + * best problem transformation should be exploited before starting the + * search procedure and an appropriate a priori transformation should be + * applied to the problem. In particular a decision should be taken + * whether variables, which are positive by nature, should be taken in + * the log scale. A hard lower variable bound can also be realized by + * taking the square. All variables should be re-scaled such that they + * "live" in a similar search range width (for example, but not + * necessarily between zero and one), such that the initial standard + * deviation can be chosen the same for all variables. + * + * + * \section links Links + * - http://www.lri.fr/~hansen/cmaesintro.html + * - http://www.lri.fr/~hansen/publications.html + * + * \section tut Tutorial + * - http://www.lri.fr/~hansen/cmatutorial.pdf + * + * \section references References + * + * - Hansen, N, and S. Kern (2004). Evaluating the CMA Evolution + * Strategy on Multimodal Test Functions. In: Eighth International + * Conference on Parallel Problem Solving from Nature PPSN VIII, + * Proceedings, pp. 282-291, Berlin: Springer + * + * - Hansen, N., S.D. Müller and P. Koumoutsakos (2003): Reducing the + * Time Complexity of the Derandomized Evolution Strategy with + * Covariance Matrix Adaptation (CMA-ES). Evolutionary Computation, + * 11(1). + * + * - Hansen, N. and A. Ostermeier (2001). Completely Derandomized + * Self-Adaptation in Evolution Strategies. Evolutionary Computation, + * 9(2), pp. 159-195. + * + * - Hansen, N. and A. Ostermeier (1996). Adapting arbitrary normal + * mutation distributions in evolution strategies: The covariance + * matrix adaptation. In Proceedings of the 1996 IEEE International + * Conference on Evolutionary Computation, pp. 312-317. + */ + +#pragma once + +#include "parameters.h" +#include "random.h" +#include "timings.h" +#include "utils.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @class CMAES + * Evolution Strategies with Covariance Matrix Adaptation. The public interface + * of the optimization algorithm. + */ +template class CMAES { +public: + /** + * Keys for get(). + */ + enum GetScalar { + NoScalar = 0, + AxisRatio = 1, + Eval = 2, + Evaluations = 2, + FctValue = 3, + FuncValue = 3, + FunValue = 3, + Fitness = 3, + FBestEver = 4, + Generation = 5, + Iteration = 5, + MaxEval = 6, + MaxFunEvals = 6, + StopMaxFunEvals = 6, + MaxGen = 7, + MaxIter = 7, + StopMaxIter = 7, + MaxAxisLength = 8, + MinAxisLength = 9, + MaxStdDev = 10, + MinStdDev = 11, + Dim = 12, + Dimension = 12, + Lambda = 13, + SampleSize = 13, + PopSize = 13, + Sigma = 14 + }; + + /** + * Keys for getPtr(). + */ + enum GetVector { + NoVector = 0, + DiagC = 1, + DiagD = 2, + StdDev = 3, + XBestEver = 4, + XBest = 5, + XMean = 6 + }; + + /** + * Keys for writeToFile(). + */ + enum WriteKey { + WCNone = 0, + WKResume = 1, + WKXMean = 2, + WKC = 4, + WKAll = 8, + WKFewInfo = 16, + WKFew = 32, + WKEval = 64, + WKFitness = 128, + WKFBestEver = 256, + WKCGeneration = 512, + WKSigma = 1024, + WKLambda = 2048, + WKB = 4096, + WKXBest = 8192, + WKClock = 16384, + WKDim = 32768 + }; + +private: + //! Implementation version. + std::string version; + //!< Random number generator. + Random rand; + //!< CMA-ES parameters. + Parameters params; + + //! Step size. + T sigma; + //! Mean x vector, "parent". + T *xmean; + //! Best sample ever. + T *xBestEver; + //! x-vectors, lambda offspring. + T **population; + //! Sorting index of sample population. + int *index; + //! History of function values. + T *funcValueHistory; + + T chiN; + //! Lower triangular matrix: i>=j for C[i][j]. + T **C; + //! Matrix with normalize eigenvectors in columns. + T **B; + //! Axis lengths. + T *rgD; + + //! Anisotropic evolution path (for covariance). + T *pc; + //! Isotropic evolution path (for step length). + T *ps; + //! Last mean. + T *xold; + //! Output vector. + T *output; + //! B*D*z. + T *BDz; + //! Temporary (random) vector used in different places. + T *tempRandom; + //! Objective function values of the population. + T *functionValues; + //!< Public objective function value array returned by init(). + T *publicFitness; + + //! Generation number. + T gen; + //! Algorithm state. + enum { INITIALIZED, SAMPLED, UPDATED } state; + + // repeatedly used for output + T maxdiagC; + T mindiagC; + T maxEW; + T minEW; + + bool eigensysIsUptodate; + bool doCheckEigen; //!< control via signals.par + T genOfEigensysUpdate; + + T dMaxSignifKond; + T dLastMinEWgroesserNull; + + bool isResumeDone; + + time_t printtime; + time_t writetime; //!< ideally should keep track for each output file + time_t firstwritetime; + time_t firstprinttime; + + std::string + stopMessage; //!< A message that contains all matched stop criteria. + + std::string getTimeStr(void) { + time_t tm = time(0); + std::string timeStr(ctime(&tm)); + return timeStr.substr(0, 24); + } + + /** + * Calculating eigenvalues and vectors. + * @param rgtmp (input) N+1-dimensional vector for temporal use. + * @param diag (output) N eigenvalues. + * @param Q (output) Columns are normalized eigenvectors. + */ + void eigen(T *diag, T **Q, T *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) + Q[i][j] = Q[j][i] = C[i][j]; + } + + householder(Q, diag, rgtmp); + ql(diag, rgtmp, Q); + } + + /** + * Exhaustive test of the output of the eigendecomposition, needs O(n^3) + * 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 res = 0; + for (int i = 0; i < params.N; ++i) + for (int j = 0; j < params.N; ++j) { + T 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]; + } + // 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); + const bool cond2 = + fabs(cc - C[i > j ? i : j][i > j ? j : i]) > T(3e-14); + if (cond1 && cond2) { + std::stringstream s; + s << i << " " << j << ": " << cc << " " + << C[i > j ? i : j][i > j ? j : i] << ", " + << cc - C[i > j ? i : j][i > j ? j : i]; + if (params.logWarnings) + params.logStream << "eigen(): imprecise result detected " << s.str() + << std::endl; + ++res; + } + if (std::fabs(dd - (i == j)) > T(1e-10)) { + std::stringstream s; + s << i << " " << j << " " << dd; + if (params.logWarnings) + params.logStream + << "eigen(): imprecise result detected (Q not orthog.)" + << s.str() << std::endl; + ++res; + } + } + return res; + } + + /** + * Symmetric tridiagonal QL algorithm, iterative. + * Computes the eigensystem from a tridiagonal matrix in roughtly 3N^3 + * operations code adapted from Java JAMA package, function tql2. + * @param d input: Diagonale of tridiagonal matrix. output: eigenvalues. + * @param e input: [1..n-1], off-diagonal, output from Householder + * @param V input: matrix output of Householder. output: basis of + * eigenvectors, according to d + */ + void ql(T *d, T *e, T **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 + + // shift input e + T *ep1 = e; + for (T *ep2 = e + 1, *const end = e + n; ep2 != end; ep1++, ep2++) + *ep1 = *ep2; + *ep1 = T(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); + if (tst1 < smallSDElement) + tst1 = smallSDElement; + const T epsTst1 = eps * tst1; + int m = l; + while (m < n) { + if (std::fabs(e[m]) <= epsTst1) + break; + m++; + } + + // 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)); + + // compute implicit shift + if (p < 0) + r = -r; + const T pr = p + r; + dl = el / pr; + h = g - dl; + const T dl1 = el * pr; + dl1r = dl1; + for (int i = l + 2; i < n; i++) + d[i] -= h; + f += h; + + // 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); + for (int i = m - 1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + const T &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]; + 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]; + h = Vki1; + T &Vki = V[k][i]; + Vki1 = s * Vki + c * h; + Vki *= c; + Vki -= s * h; + } + } + p = -s * s2 * c3 * el1 * el / dl1; + el = s * p; + dl = c * p; + } while (std::fabs(el) > epsTst1); + } + dl += f; + el = 0.0; + } + } + + /** + * Householder transformation of a symmetric matrix V into tridiagonal form. + * Code slightly adapted from the Java JAMA package, function private tred2(). + * @param V input: symmetric nxn-matrix. output: orthogonal transformation + * matrix: tridiag matrix == V* V_in* V^t. + * @param d output: diagonal + * @param e output: [0..n-1], off diagonal (elements 1..n-1) + */ + void householder(T **V, T *d, T *e) { + const int n = params.N; + + for (int j = 0; j < n; j++) { + d[j] = V[n - 1][j]; + } + + // Householder reduction to tridiagonal form + + 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++) { + scale += std::fabs(*pd); + } + if (scale == 0.0) { + e[i] = d[i - 1]; + for (int j = 0; j < i; j++) { + d[j] = V[i - 1][j]; + V[i][j] = 0.0; + V[j][i] = 0.0; + } + } else { + // generate Householder vector + for (T *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); + e[i] = scale * g; + h = h - f * g; + dim1 = f - g; + memset((void *)e, 0, (size_t)i * sizeof(T)); + + // apply similarity transformation to remaining columns + for (int j = 0; j < i; j++) { + f = d[j]; + V[j][i] = f; + T &ej = e[j]; + g = ej + V[j][j] * f; + for (int k = j + 1; k <= i - 1; k++) { + T &Vkj = V[k][j]; + g += Vkj * d[k]; + e[k] += Vkj * f; + } + ej = g; + } + f = 0.0; + for (int j = 0; j < i; j++) { + T &ej = e[j]; + ej /= h; + f += ej * d[j]; + } + T 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]; + f = dj; + g = e[j]; + for (int k = j; k <= i - 1; k++) { + V[k][j] -= f * e[k] + g * d[k]; + } + dj = V[i - 1][j]; + V[i][j] = 0.0; + } + } + d[i] = h; + } + + // accumulate transformations + const int nm1 = n - 1; + for (int i = 0; i < nm1; i++) { + T h; + T &Vii = V[i][i]; + V[n - 1][i] = Vii; + Vii = 1.0; + h = d[i + 1]; + if (h != 0.0) { + for (int k = 0; k <= i; k++) { + d[k] = V[k][i + 1] / h; + } + for (int j = 0; j <= i; j++) { + T g = 0.0; + for (int k = 0; k <= i; k++) { + T *Vk = V[k]; + g += Vk[i + 1] * Vk[j]; + } + for (int k = 0; k <= i; k++) { + V[k][j] -= g * d[k]; + } + } + } + for (int k = 0; k <= i; k++) { + V[k][i + 1] = 0.0; + } + } + for (int j = 0; j < n; j++) { + T &Vnm1j = V[n - 1][j]; + d[j] = Vnm1j; + Vnm1j = 0.0; + } + V[n - 1][n - 1] = 1.0; + e[0] = 0.0; + } + + /** + * Dirty index sort. + */ + void sortIndex(const T *rgFunVal, int *iindex, int n) { + int i, j; + for (i = 1, iindex[0] = 0; i < n; ++i) { + for (j = i; j > 0; --j) { + if (rgFunVal[iindex[j - 1]] < rgFunVal[i]) + break; + iindex[j] = iindex[j - 1]; // shift up + } + iindex[j] = i; + } + } + + void adaptC2(const int hsig) { + const int N = params.N; + bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; + + if (params.ccov != T(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); + + 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]; + 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]]; + Cij += ccovmu * params.weights[k] * (rgrgxindexk[i] - xold[i]) * + (rgrgxindexk[j] - xold[j]) / sigmasquare; + } + } + // update maximal and minimal diagonal value + maxdiagC = mindiagC = C[0][0]; + for (int i = 1; i < N; ++i) { + const T &Cii = C[i][i]; + if (maxdiagC < Cii) + maxdiagC = Cii; + else if (mindiagC > Cii) + mindiagC = Cii; + } + } + } + + /** + * Treats minimal standard deviations and numeric problems. Increases sigma. + */ + void testMinStdDevs(void) { + if (!this->params.rgDiffMinChange) + return; + + 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); + } + + /** + * Adds the mutation sigma*B*(D*z). + * @param x Search space vector. + * @param eps Mutation factor. + */ + void addMutation(T *x, T 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; + for (int j = 0; j < params.N; ++j) + sum += B[i][j] * tempRandom[j]; + x[i] = xmean[i] + eps * sigma * sum; + } + } + + /** + * This hack reads key words from input key for data to be written to + * a file, see file signals.par as input file. The length of the keys + * is mostly fixed. If the key phrase does not match the expectation the + * output might be strange. + */ + void writeToStream(int key, std::ostream &file) { + if (key & WKResume) { + file << std::endl << "# resume " << params.N << std::endl; + file << "xmean" << std::endl; + writeToStream(WKXMean, file); + file << "path for sigma" << std::endl; + for (int i = 0; i < params.N; ++i) + file << ps[i] << (i == params.N - 1 ? "\n" : "\t"); + file << "path for C" << std::endl; + for (int i = 0; i < params.N; ++i) + file << pc[i] << (i == params.N - 1 ? "\n" : "\t"); + file << "sigma " << sigma << std::endl; + // note than B and D might not be up-to-date + file << "covariance matrix" << std::endl; + writeToStream(WKC, file); + } + if (key & WKXMean) { + for (int i = 0; i < params.N; ++i) + file << (i == 0 ? "" : "\t") << xmean[i]; + file << std::endl; + } + if (key & WKC) { + for (int i = 0; i < params.N; ++i) + for (int j = 0; j <= i; ++j) { + file << C[i][j]; + if (j == i) + file << std::endl; + else + file << '\t'; + } + file << std::endl; + } + if (key & WKAll) { + time_t ti = time(0); + file << std::endl + << "# --------- " << asctime(localtime(&ti)) << std::endl; + file << " N " << params.N << std::endl; + 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 << "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)) + << 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] + << (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] + << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); + file << "xmean" << std::endl; + for (int 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) + 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) + 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] + << (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) + 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) + file << " " << std::setw(12) << B[i][k] + << (i % 5 == 4 || i == params.N - 1 ? '\n' : ' '); + file << std::endl; + } + if (key & WKFewInfo) { + file << " Iter\tFevals\tFunction Value\tSigma\tMaxCoorDev\tMinCoorDev\t" + << "AxisRatio\tMinDii\tTime in eig" << std::endl; + file << std::endl; + } + if (key & WKFew) { + file << (int)gen << "\t" << (int)countevals << "\t" + << functionValues[index[0]] << "\t\t" << sigma << " " + << sigma * std::sqrt(maxdiagC) << "\t" << sigma * std::sqrt(mindiagC) + << "\t" << std::scientific << std::setprecision(2) + << std::sqrt(maxEW / minEW) << "\t" << std::sqrt(minEW) << " " + << eigenTimings.totaltotaltime; + file << std::endl; + } + if (key & WKEval) { + file << countevals; + file << std::endl; + } + if (key & WKFitness) { + for (int i = 0; i < params.N; ++i) + file << (i == 0 ? "" : "\t") << functionValues[index[i]]; + file << std::endl; + } + if (key & WKFBestEver) { + file << xBestEver[params.N] << std::endl; + } + if (key & WKCGeneration) { + file << gen << std::endl; + } + if (key & WKSigma) { + file << sigma << std::endl; + } + if (key & WKLambda) { + file << params.lambda << std::endl; + } + if (key & WKB) { + int *iindex = new int[params.N]; + sortIndex(rgD, iindex, 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]]; + if (j != params.N - 1) + file << '\t'; + else + file << std::endl; + } + delete[] iindex; + iindex = 0; + file << std::endl; + } + if (key & WKXBest) { + for (int i = 0; i < params.N; ++i) + file << (i == 0 ? "" : "\t") << population[index[0]][i]; + file << std::endl; + } + if (key & WKClock) { + eigenTimings.update(); + file << eigenTimings.totaltotaltime << " " << eigenTimings.tictoctime + << std::endl; + } + if (key & WKDim) { + file << params.N; + file << std::endl; + } + } + +public: + T countevals; //!< objective function evaluations + Timing eigenTimings; + + CMAES() : version("1.0alpha") {} + + /** + * Releases the dynamically allocated memory, including that of the return + * value of init(). + */ + ~CMAES() { + delete[] pc; + delete[] ps; + delete[] tempRandom; + delete[] BDz; + 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]; + delete[] population; + delete[] C; + delete[] B; + delete[] index; + delete[] publicFitness; + delete[]-- functionValues; + delete[]-- funcValueHistory; + } + + /** + * 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. + */ + T *init(const Parameters ¶meters) { + params = parameters; + + stopMessage = ""; + + T 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)); + 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)) + break; + dMaxSignifKond = dtest / T(1000); // not sure whether this is really save, + // 100 does not work well enough + + gen = 0; + countevals = 0; + state = INITIALIZED; + dLastMinEWgroesserNull = T(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 + 1]; + 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++; + + for (int i = 0; i < params.N; ++i) { + C[i] = new T[i + 1]; + B[i] = new T[params.N]; + } + index = new int[params.lambda]; + for (int i = 0; i < params.lambda; ++i) + index[i] = i; + population = new T *[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; + } + + // initialize newed space + for (int i = 0; i < params.lambda; i++) { + functionValues[i] = std::numeric_limits::max(); + } + for (int i = 0; i < historySize; i++) { + 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); + C[i][i] *= C[i][i]; + pc[i] = ps[i] = T(0); + } + minEW = minElement(rgD, params.N); + minEW = minEW * minEW; + maxEW = maxElement(rgD, params.N); + maxEW = maxEW * maxEW; + + maxdiagC = C[0][0]; + for (int i = 1; i < params.N; ++i) + if (maxdiagC < C[i][i]) + maxdiagC = C[i][i]; + mindiagC = C[0][0]; + for (int i = 1; i < params.N; ++i) + if (mindiagC > C[i][i]) + mindiagC = C[i][i]; + + for (int i = 0; i < params.N; ++i) + xmean[i] = xold[i] = params.xstart[i]; + // use in case xstart as typicalX + if (params.typicalXcase) + for (int i = 0; i < params.N; ++i) + xmean[i] += sigma * rgD[i] * rand.gauss(); + + if (params.resumefile != "") + resumeDistribution(params.resumefile); + + return publicFitness; + } + + /** + * Well, says hello. + * @return eg. "(5,10)-CMA-ES(mu_eff=3.4), Ver="1.0alpha", dimension=9" + */ + 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() << ")"; + return stream.str(); + } + + /** + * Allows to restart with saved internal state (distribution) variables (use + * writeToFile() for saving). Keyword "resume" followed by a filename in + * initials.par invokes this function during initialization. Searches in + * filename for the last occurrence of word "resume", followed by a dimension + * number, and reads the subsequent values for xmean, evolution paths ps and + * pc, sigma and covariance matrix. Note that init() needs to be called + * before calling resume_distribution() explicitely. In the former all the + * remaining (strategy-)parameters are set. It can be useful to edit the + * written parameters, in particular to increase sigma, before resume. + * + * Not all internal state parameters are recovered. In particular generation + * number and xbestever are not restored. For covariance matrices with large + * condition numbers the writing precision of 6 digits is not sufficient and + * resume will lead to poor result. + * @param filename A file, that was written presumably by writeToFile(). + */ + void resumeDistribution(const std::string &filename) { + std::ifstream file(filename.c_str()); + if (!file.is_open()) + throw std::runtime_error("resumeDistribution(): could not open '" + + filename + "'"); + + std::streampos lastResume = 0; + std::string entry = ""; + while (!file.eof()) { + file >> entry; + if (entry == "resume") { + lastResume = file.tellg(); + break; + } + } + file.clear(); + file.seekg(lastResume); + + int n = 0; + file >> n; + if (n != params.N) + throw std::runtime_error( + "resumeDistribution(): Dimension numbers do not match"); + + // find next "xmean" entry + while (!file.eof()) { + file >> entry; + if (entry == "xmean") + break; + } + // read xmean + if (file.eof()) + throw std::runtime_error("resumeDistribution(): 'xmean' not found"); + for (int i = 0; i < n; i++) + file >> xmean[i]; + file.clear(); + file.seekg(lastResume); + + // find next "path for sigma" entry + while (!file.eof()) { + file >> entry; + if (entry == "path") { + std::string temp = ""; + file >> temp; + entry += " " + temp; + file >> temp; + entry += " " + temp; + if (entry == "path for sigma") + break; + } + } + // read ps + if (file.eof()) + throw std::runtime_error( + "resumeDistribution(): 'path for sigma' not found"); + for (int i = 0; i < n; i++) + file >> ps[i]; + file.clear(); + file.seekg(lastResume); + + // find next "path for C" entry + while (!file.eof()) { + file >> entry; + if (entry == "path") { + std::string temp = ""; + file >> temp; + entry += " " + temp; + file >> temp; + entry += " " + temp; + if (entry == "path for C") + break; + } + } + // read pc + if (file.eof()) + throw std::runtime_error("resumeDistribution(): 'path for C' not found"); + for (int i = 0; i < n; i++) + file >> pc[i]; + file.clear(); + file.seekg(lastResume); + + // find next "sigma" entry + while (!file.eof()) { + file >> entry; + if (entry == "sigma") + break; + } + // read pc + if (file.eof()) + throw std::runtime_error("resumeDistribution(): 'sigma' not found"); + file >> sigma; + file.clear(); + file.seekg(lastResume); + + // find next "covariance matrix" entry + while (!file.eof()) { + file >> entry; + if (entry == "covariance") { + std::string temp = ""; + file >> temp; + entry += " " + temp; + if (entry == "covariance matrix") + break; + } + } + // read C + if (file.eof()) + throw std::runtime_error( + "resumeDistribution(): 'covariance matrix' not found"); + for (int i = 0; i < params.N; ++i) + for (int j = 0; j <= i; ++j) + file >> C[i][j]; + + eigensysIsUptodate = false; + isResumeDone = true; + updateEigensystem(true); + } + + /** + * 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. + */ + T *const *samplePopulation() { + bool diag = params.diagonalCov == 1 || params.diagonalCov >= gen; + + // calculate eigensystem + if (!eigensysIsUptodate) { + if (!diag) + updateEigensystem(false); + 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)); + eigensysIsUptodate = true; + eigenTimings.start(); + } + } + + testMinStdDevs(); + + 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) + 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) + { + T sum = 0.0; + for (int j = 0; j < params.N; ++j) + sum += B[i][j] * tempRandom[j]; + rgrgxink[i] = xmean[i] + sigma * sum; + } + } + + if (state == UPDATED || gen == 0) + ++gen; + state = SAMPLED; + + return population; + } + + /** + * Can be called after samplePopulation() to resample single solutions of the + * population as often as desired. Useful to implement a box constraints + * (boundary) handling. + * @param i Index to an element of the returned value of samplePopulation(). + * population[index] will be resampled where \f$0\leq i<\lambda\f$ + * must hold. + * @return A pointer to the resampled "population". + */ + T *const *reSampleSingle(int i) { + T *x; + assert(i >= 0 && i < params.lambda && + "reSampleSingle(): index must be between 0 and sp.lambda"); + x = population[i]; + addMutation(x); + return population; + } + + /** + * Can be called after samplePopulation() to resample single solutions. In + * general, the function can be used to sample as many independent + * mean+sigma*Normal(0,C) distributed vectors as desired. + * + * Input x can be a pointer to an element of the vector returned by + * samplePopulation() but this is inconsistent with the const qualifier of the + * returned value and therefore rather reSampleSingle() should be used. + * @param x Solution vector that gets sampled a new value. If x == NULL new + * memory is allocated and must be released by the user using + * delete[]. + * @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]; + addMutation(x); + return x; + } + + /** + * Can be called after samplePopulation() to resample single solutions. In + * general, the function can be used to sample as many independent + * mean+sigma*Normal(0,C) distributed vectors as desired. + * @param x Element of the return value of samplePopulation(), that is + * pop[0..\f$\lambda\f$]. This solution vector of the population gets + * sampled a new value. + * @return A pointer to the resampled "population" member. + */ + T const *reSampleSingleOld(T *x) { + assert(x && "reSampleSingleOld(): Missing input x"); + addMutation(x); + return x; + } + + /** + * Used to reevaluate a slightly disturbed solution for an uncertaintly + * measurement. In case if x == NULL on input, the memory of the returned x + * must be released. + * @param x Solution vector that gets sampled a new value. If x == NULL new + * memory is allocated and must be released by the user using + * delete[] x. + * @param pxmean Mean vector \f$\mu\f$ for perturbation. + * @param eps Scale factor \f$\epsilon\f$ for perturbation: + * \f$x \sim \mu + \epsilon \sigma N(0,C)\f$. + * @return A pointer to the perturbed solution vector, equals input x for + * x != NULL. + */ + T *perturbSolutionInto(T *x, T const *pxmean, T eps) { + if (!x) + x = new T[params.N]; + assert(pxmean && "perturbSolutionInto(): pxmean was not given"); + addMutation(x, eps); + return x; + } + + /** + * 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. + * @return Mean value of the new distribution. + */ + T *updateDistribution(const T *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."); + + if (state == SAMPLED) // function values are delivered here + countevals += params.lambda; + else if (params.logWarnings) + 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]; + + // Generate index + sortIndex(fitnessValues, index, params.lambda); + + // 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); + if (params.logWarnings) { + params.logStream + << "Warning: sigma increased due to equal function values" + << std::endl + << " Reconsider the formulation of the objective function"; + } + } + + // update function value history + for (int i = (int)*(funcValueHistory - 1) - 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; + } + + const T sqrtmueffdivsigma = std::sqrt(params.mueff) / sigma; + // calculate xmean and rgBDz~N(0,C) + for (int 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]; + 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; + if (diag) + sum = BDz[i]; + else { + sum = 0.; + for (int j = 0; j < N; ++j) + sum += B[j][i] * BDz[j]; + } + tempRandom[i] = sum / rgD[i]; + } + + // 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; + for (int i = 0; i < N; ++i) { + T sum; + if (diag) + sum = tempRandom[i]; + else { + sum = T(0); + T *Bi = B[i]; + for (int j = 0; j < N; ++j) + sum += Bi[j] * tempRandom[j]; + } + ps[i] = invps * ps[i] + sqrtFactor * sum; + } + + // calculate norm(ps)^2 + T psxps(0); + for (int i = 0; i < N; ++i) { + const T &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)) / + 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)); + for (int i = 0; i < N; ++i) + pc[i] = ccumcovinv * pc[i] + hsigFactor * BDz[i]; + + // update of C + adaptC2(hsig); + + // update of sigma + sigma *= + std::exp(((std::sqrt(psxps) / chiN) - T(1)) * params.cs / params.damps); + + state = UPDATED; + return xmean; + } + + /** + * Request a scalar parameter from CMA-ES. + * @param key Key of the requested scalar. + * @return The desired value. + */ + T get(GetScalar key) { + switch (key) { + case AxisRatio: + return maxElement(rgD, params.N) / minElement(rgD, params.N); + case Eval: + return countevals; + case Fitness: + return functionValues[index[0]]; + case FBestEver: + return xBestEver[params.N]; + case Generation: + return gen; + case MaxEval: + return params.stopMaxFunEvals; + case MaxIter: + return std::ceil(params.stopMaxIter); + case MaxAxisLength: + return sigma * std::sqrt(maxEW); + case MinAxisLength: + return sigma * std::sqrt(minEW); + case MaxStdDev: + return sigma * std::sqrt(maxdiagC); + case MinStdDev: + return sigma * std::sqrt(mindiagC); + case Dimension: + return params.N; + case SampleSize: + return params.lambda; + case Sigma: + return sigma; + default: + throw std::runtime_error("get(): No match found for key"); + } + } + + /** + * Request a vector parameter from CMA-ES. + * @param key Key of the requested vector. + * @return Pointer to the desired value array. Its content might be + * overwritten during the next call to any member functions other + * than get(). + */ + const T *getPtr(GetVector key) { + switch (key) { + case DiagC: { + for (int i = 0; i < params.N; ++i) + output[i] = C[i][i]; + return output; + } + case DiagD: + return rgD; + case StdDev: { + for (int i = 0; i < params.N; ++i) + output[i] = sigma * std::sqrt(C[i][i]); + return output; + } + case XBestEver: + return xBestEver; + case XBest: + return population[index[0]]; + case XMean: + return xmean; + default: + throw std::runtime_error("getPtr(): No match found for key"); + } + } + + /** + * Request a vector parameter from CMA-ES. + * @param key Key of the requested vector. + * @return Pointer to the desired value array with unlimited reading and + * writing access to its elements. The memory must be explicitly + * released using delete[]. + */ + T *getNew(GetVector key) { return getInto(key, 0); } + + /** + * Request a vector parameter from CMA-ES. + * @param key Key of the requested vector. + * @param res Memory of size N == dimension, where the desired values are + * 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); + if (!res) + res = new T[params.N]; + for (int i = 0; i < params.N; ++i) + res[i] = res0[i]; + return res; + } + + /** + * Some stopping criteria can be set in initials.par, with names starting + * with stop... Internal stopping criteria include a maximal condition number + * of about 10^15 for the covariance matrix and situations where the numerical + * discretisation error in x-space becomes noticeably. You can get a message + * that contains the matched stop criteria via getStopMessage(). + * @return Does any stop criterion match? + */ + bool testForTermination() { + T range, fac; + int iAchse, iKoo; + int diag = params.diagonalCov == 1 || params.diagonalCov >= gen; + int N = params.N; + std::stringstream message; + + 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; + } + + // TolFun + range = std::max(maxElement(funcValueHistory, + (int)std::min(gen, *(funcValueHistory - 1))), + maxElement(functionValues, params.lambda)) - + std::min(minElement(funcValueHistory, + (int)std::min(gen, *(funcValueHistory - 1))), + minElement(functionValues, params.lambda)); + + if (gen > 0 && range <= params.stopTolFun) { + message << "TolFun: function value differences " << range + << " < stopTolFun=" << params.stopTolFun << std::endl; + } + + // TolFunHist + if (gen > *(funcValueHistory - 1)) { + range = maxElement(funcValueHistory, (int)*(funcValueHistory - 1)) - + minElement(funcValueHistory, (int)*(funcValueHistory - 1)); + if (range <= params.stopTolFunHist) + message << "TolFunHist: history of function value changes " << range + << " stopTolFunHist=" << params.stopTolFunHist << std::endl; + } + + // 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; + } + + // TolUpX + for (int i = 0; i < N; ++i) { + if (sigma * std::sqrt(C[i][i]) > + params.stopTolUpXFactor * params.rgInitialStds[i]) { + message << "TolUpX: standard deviation increased by more than " + << params.stopTolUpXFactor + << ", larger initial standard deviation recommended." + << std::endl; + break; + } + } + + // 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; + } + + // 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]) + break; + } + if (iKoo == N) { + message << "NoEffectAxis: standard deviation 0.1*" << (fac / 0.1) + << " in principal axis " << iAchse << " without effect" + << std::endl; + break; + } + } + } + // 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)) { + message << "NoEffectCoordinate: standard deviation 0.2*" + << (sigma * std::sqrt(C[iKoo][iKoo])) << " in coordinate " + << iKoo << " without effect" << std::endl; + break; + } + } + + if (countevals >= params.stopMaxFunEvals) { + message << "MaxFunEvals: conducted function evaluations " << countevals + << " >= " << params.stopMaxFunEvals << std::endl; + } + if (gen >= params.stopMaxIter) { + message << "MaxIter: number of iterations " << gen + << " >= " << params.stopMaxIter << std::endl; + } + + stopMessage = message.str(); + return stopMessage != ""; + } + + /** + * A message that contains a detailed description of the matched stop + * criteria. + */ + std::string getStopMessage() { return stopMessage; } + + /** + * @param filename Output file name. + * @param key Key of type WriteKey that indicates the content that should be + * written. You can combine multiple keys with |. + */ + void writeToFile(int key, const std::string &filename) { + std::ofstream file; + file.open(filename.c_str(), std::ios_base::app); + + if (file.is_open()) { + if (gen > 0 || filename.substr(0, 11) != "outcmaesfit") + writeToStream(key, file); /* do not write fitness for gen==0 */ + file.close(); + } else { + throw std::runtime_error("writeToFile(): could not open '" + filename + + "'"); + } + } + + /** + * 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$ + * 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. + */ + void updateEigensystem(bool force) { + eigenTimings.update(); + + if (!force) { + if (eigensysIsUptodate) + return; + // return on modulo generation number + if (gen < genOfEigensysUpdate + params.updateCmode.modulo) + return; + // return on time percentage + if (params.updateCmode.maxtime < 1.00 && + eigenTimings.tictoctime > + params.updateCmode.maxtime * eigenTimings.totaltime && + eigenTimings.tictoctime > 0.0002) + return; + } + + eigenTimings.tic(); + eigen(rgD, B, tempRandom); + eigenTimings.toc(); + + // find largest and smallest eigenvalue, they are supposed to be sorted + // anyway + minEW = minElement(rgD, params.N); + maxEW = maxElement(rgD, params.N); + + if (doCheckEigen) // needs O(n^3)! writes, in case, error message in error + // file + checkEigen(rgD, B); + + for (int i = 0; i < params.N; ++i) + rgD[i] = std::sqrt(rgD[i]); + + eigensysIsUptodate = true; + genOfEigensysUpdate = gen; + } + + /** + * Distribution mean could be changed before samplePopulation(). This might + * lead to unexpected behaviour if done repeatedly. + * @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) { + assert(state != SAMPLED && "setMean: mean cannot be set inbetween the calls" + "of samplePopulation and updateDistribution"); + + if (newxmean && newxmean != xmean) + for (int i = 0; i < params.N; ++i) + xmean[i] = newxmean[i]; + else + newxmean = xmean; + + return newxmean; + } +}; diff --git a/include/cma-es/parameters.h b/include/cma-es/parameters.h new file mode 100644 index 0000000..501aa57 --- /dev/null +++ b/include/cma-es/parameters.h @@ -0,0 +1,391 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +template class CMAES; + +/** + * @class Parameters + * Holds all parameters that can be adjusted by the user. + */ +template class Parameters { + friend class CMAES; + +public: + /* Input parameters. */ + //! Problem dimension, must stay constant. + int N; + //! Initial search space vector. + T *xstart; + //! A typical value for a search space vector. + T *typicalX; + //! Indicates that the typical x is the initial point. + bool typicalXcase; + //! Initial standard deviations. + T *rgInitialStds; + T *rgDiffMinChange; + + /* Termination parameters. */ + //! Maximal number of objective function evaluations. + T stopMaxFunEvals; + T facmaxeval; + //! Maximal number of iterations. + T stopMaxIter; + //! Minimal fitness value. Only activated if flg is true. + struct { + bool flg; + T val; + } stStopFitness; + //! Minimal value difference. + T stopTolFun; + //! Minimal history value difference. + T stopTolFunHist; + //! Minimal search space step size. + T stopTolX; + //! Defines the maximal condition number. + T stopTolUpXFactor; + + /* internal evolution strategy parameters */ + /** + * Population size. Number of samples per iteration, at least two, + * generally > 4. + */ + int lambda; + /** + * Number of individuals used to recompute the mean. + */ + int mu; + T mucov; + /** + * Variance effective selection mass, should be lambda/4. + */ + T mueff; + /** + * Weights used to recombinate the mean sum up to one. + */ + T *weights; + /** + * Damping parameter for step-size adaption, d = inifinity or 0 means adaption + * is turned off, usually close to one. + */ + T damps; + /** + * cs^-1 (approx. n/3) is the backward time horizon for the evolution path + * ps and larger than one. + */ + T cs; + T ccumcov; + /** + * ccov^-1 (approx. n/4) is the backward time horizon for the evolution path + * pc and larger than one. + */ + T ccov; + T diagonalCov; + struct { + T modulo; + T maxtime; + } updateCmode; + T facupdateCmode; + + /** + * Determines the method used to initialize the weights. + */ + enum Weights { + UNINITIALIZED_WEIGHTS, + LINEAR_WEIGHTS, + EQUAL_WEIGHTS, + LOG_WEIGHTS + } weightMode; + + //! File that contains an optimization state that should be resumed. + std::string resumefile; + + //! Set to true to activate logging warnings. + bool logWarnings; + //! Output stream that is used to log warnings, usually std::cerr. + std::ostream &logStream; + + Parameters() + : N(-1), xstart(0), typicalX(0), typicalXcase(false), rgInitialStds(0), + rgDiffMinChange(0), stopMaxFunEvals(-1), facmaxeval(1.0), + stopMaxIter(-1.0), stopTolFun(1e-12), stopTolFunHist(1e-13), + stopTolX(0), // 1e-11*insigma would also be reasonable + stopTolUpXFactor(1e3), lambda(-1), mu(-1), mucov(-1), mueff(-1), + weights(0), damps(-1), cs(-1), ccumcov(-1), ccov(-1), facupdateCmode(1), + weightMode(UNINITIALIZED_WEIGHTS), resumefile(""), logWarnings(false), + logStream(std::cerr) { + stStopFitness.flg = false; + stStopFitness.val = -std::numeric_limits::max(); + updateCmode.modulo = -1; + updateCmode.maxtime = -1; + } + + Parameters(const Parameters ¶meters) { assign(parameters); } + + ~Parameters() { + if (xstart) + delete[] xstart; + if (typicalX) + delete[] typicalX; + if (rgInitialStds) + delete[] rgInitialStds; + if (rgDiffMinChange) + delete[] rgDiffMinChange; + if (weights) + delete[] weights; + } + + Parameters &operator=(const Parameters ¶meters) { + assign(parameters); + return *this; + } + + /** + * @param dimension Dimension of the search space \f$N\f$. No default + * 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$. + * @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$. + */ + void init(int dimension = 0, const T *inxstart = 0, const T *inrgsigma = 0) { + if (logWarnings) { + if (!(xstart || inxstart || typicalX)) + logStream << "Warning: initialX undefined. typicalX = 0.5...0.5." + << std::endl; + if (!(rgInitialStds || inrgsigma)) + logStream << "Warning: initialStandardDeviations undefined. 0.3...0.3." + << std::endl; + } + + if (dimension <= 0 && N <= 0) + throw std::runtime_error("Problem dimension N undefined."); + else if (dimension > 0) + N = dimension; + + if (weightMode == UNINITIALIZED_WEIGHTS) + weightMode = LOG_WEIGHTS; + + diagonalCov = 0; // default is 0, but this might change in future + + if (!xstart) { + xstart = new T[N]; + if (inxstart) { + for (int i = 0; i < N; ++i) + xstart[i] = inxstart[i]; + } else if (typicalX) { + typicalXcase = true; + for (int i = 0; i < N; ++i) + xstart[i] = typicalX[i]; + } else { + typicalXcase = true; + for (int i = 0; i < N; i++) + xstart[i] = 0.5; + } + } + + if (!rgInitialStds) { + rgInitialStds = new T[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); + } + } + + supplementDefaults(); + } + +private: + void assign(const Parameters &p) { + N = p.N; + + if (xstart) + delete[] xstart; + if (p.xstart) { + xstart = new T[N]; + for (int i = 0; i < N; i++) + xstart[i] = p.xstart[i]; + } + + if (typicalX) + delete[] typicalX; + if (p.typicalX) { + typicalX = new T[N]; + for (int i = 0; i < N; i++) + typicalX[i] = p.typicalX[i]; + } + + typicalXcase = p.typicalXcase; + + if (rgInitialStds) + delete[] rgInitialStds; + if (p.rgInitialStds) { + rgInitialStds = new T[N]; + for (int i = 0; i < N; i++) + rgInitialStds[i] = p.rgInitialStds[i]; + } + + if (rgDiffMinChange) + delete[] rgDiffMinChange; + if (p.rgDiffMinChange) { + rgDiffMinChange = new T[N]; + for (int i = 0; i < N; i++) + rgDiffMinChange[i] = p.rgDiffMinChange[i]; + } + + stopMaxFunEvals = p.stopMaxFunEvals; + facmaxeval = p.facmaxeval; + stopMaxIter = p.stopMaxIter; + + stStopFitness.flg = p.stStopFitness.flg; + stStopFitness.val = p.stStopFitness.val; + + stopTolFun = p.stopTolFun; + stopTolFunHist = p.stopTolFunHist; + stopTolX = p.stopTolX; + stopTolUpXFactor = p.stopTolUpXFactor; + + lambda = p.lambda; + mu = p.mu; + mucov = p.mucov; + mueff = p.mueff; + + if (weights) + delete[] weights; + if (p.weights) { + weights = new T[mu]; + for (int i = 0; i < mu; i++) + weights[i] = p.weights[i]; + } + + damps = p.damps; + cs = p.cs; + ccumcov = p.ccumcov; + ccov = p.ccov; + diagonalCov = p.diagonalCov; + + updateCmode.modulo = p.updateCmode.modulo; + updateCmode.maxtime = p.updateCmode.maxtime; + + facupdateCmode = p.facupdateCmode; + + weightMode = p.weightMode; + + resumefile = p.resumefile; + } + + /** + * Supplements default parameter values. + */ + void supplementDefaults() { + if (lambda < 2) + lambda = 4 + (int)(3.0 * log((double)N)); + if (mu <= 0) + mu = lambda / 2; + if (!weights) + setWeights(weightMode); + + if (cs > 0) + cs *= (mueff + 2.) / (N + mueff + 3.); + if (cs <= 0 || cs >= 1) + cs = (mueff + 2.) / (N + mueff + 3.); + + if (ccumcov <= 0 || ccumcov > 1) + 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); + t2 = (t2 > 1) ? 1 : t2; + t2 = (1. / mucov) * t1 + (1. - 1. / mucov) * t2; + if (ccov >= 0) + ccov *= t2; + if (ccov < 0 || ccov > 1) + ccov = t2; + + if (diagonalCov < 0) + diagonalCov = 2 + 100. * N / sqrt((double)lambda); + + if (stopMaxFunEvals <= 0) + 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); + 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))) + + cs; + + if (updateCmode.modulo < 0) + updateCmode.modulo = 1. / ccov / (double)N / 10.; + updateCmode.modulo *= facupdateCmode; + if (updateCmode.maxtime < 0) + updateCmode.maxtime = 0.20; // maximal 20% of CPU-time + } + + /** + * Initializes the offspring weights. + */ + void setWeights(Weights mode) { + if (weights) + delete[] weights; + weights = new T[mu]; + switch (mode) { + case LINEAR_WEIGHTS: + for (int i = 0; i < mu; ++i) + weights[i] = mu - i; + break; + case EQUAL_WEIGHTS: + for (int i = 0; i < mu; ++i) + weights[i] = 1; + break; + case LOG_WEIGHTS: + default: + for (int i = 0; i < mu; ++i) + weights[i] = log(mu + 1.) - log(i + 1.); + break; + } + + // normalize weights vector and set mueff + T s1 = 0, s2 = 0; + for (int i = 0; i < mu; ++i) { + s1 += weights[i]; + s2 += weights[i] * weights[i]; + } + mueff = s1 * s1 / s2; + for (int i = 0; i < mu; ++i) + weights[i] /= s1; + + if (mu < 1 || mu > lambda || + (mu == lambda && weights[0] == weights[mu - 1])) + throw std::runtime_error("setWeights(): invalid setting of mu or lambda"); + } +}; diff --git a/include/cma-es/random.h b/include/cma-es/random.h new file mode 100644 index 0000000..50dcd34 --- /dev/null +++ b/include/cma-es/random.h @@ -0,0 +1,83 @@ +#pragma once + +#include +#include + +/** + * @class Random + * A pseudo random number generator. + */ +template class Random { + // variables for uniform() + long int startseed; + long int aktseed; + long int aktrand; + long int rgrand[32]; + // variables for gauss() + bool stored; + T hold; + +public: + /** + * @param seed use clock if 0 + */ + Random(long unsigned seed = 0) : hold(0.0) { + stored = false; + if (seed < 1) { + long int t = 100 * time(0) + clock(); + seed = (long unsigned)(t < 0 ? -t : t); + } + start(seed); + } + /** + * @param seed 0 == 1 + */ + void start(long unsigned seed) { + stored = false; + startseed = seed; + if (seed < 1) + seed = 1; + aktseed = seed; + for (int i = 39; i >= 0; --i) { + long tmp = aktseed / 127773; + aktseed = 16807 * (aktseed - tmp * 127773) - 2836 * tmp; + if (aktseed < 0) + aktseed += 2147483647; + if (i < 32) + rgrand[i] = aktseed; + } + aktrand = rgrand[0]; + } + /** + * @return (0,1)-normally distributed random number + */ + T gauss(void) { + if (stored) { + stored = false; + return hold; + } + stored = true; + T x1, x2, rquad; + do { + x1 = 2.0 * uniform() - 1.0; + x2 = 2.0 * uniform() - 1.0; + rquad = x1 * x1 + x2 * x2; + } while (rquad >= 1 || rquad <= 0); + const T fac = std::sqrt(T(-2) * std::log(rquad) / rquad); + hold = fac * x1; + return fac * x2; + } + /** + * @return (0,1)-uniform distributed random number + */ + T uniform(void) { + long tmp = aktseed / 127773; + aktseed = 16807 * (aktseed - tmp * 127773) - 2836 * tmp; + if (aktseed < 0) + aktseed += 2147483647; + tmp = aktrand / 67108865; + aktrand = rgrand[tmp]; + rgrand[tmp] = aktseed; + return (T)aktrand / T(2.147483647e9); + } +}; diff --git a/include/cma-es/timings.h b/include/cma-es/timings.h new file mode 100644 index 0000000..a255e31 --- /dev/null +++ b/include/cma-es/timings.h @@ -0,0 +1,95 @@ +#pragma once + +#include "utils.h" +#include +#include + +/** + * @class Timing + * A class for time measurements of the eigen decomposition. + * Timing measures overall time and times between calls of tic and toc. For + * small time spans (up to 1000 seconds) CPU time via clock() is used. For large + * time spans the fall-back to elapsed time from time() is used. + * timings_update() must be called often enough to prevent the fallback. + */ +class Timing { + clock_t lastclock; + time_t lasttime; + clock_t ticclock; + time_t tictime; + short istic; + short isstarted; + + double lastdiff; + double tictoczwischensumme; + +public: + double totaltime; //! zeroed by re-calling timings_start + double totaltotaltime; + double tictoctime; + double lasttictoctime; + + Timing() { + totaltotaltime = 0; + start(); + } + + void start() { + totaltime = 0; + tictoctime = 0; + lasttictoctime = 0; + istic = 0; + lastclock = clock(); + lasttime = time(NULL); + lastdiff = 0; + tictoczwischensumme = 0; + isstarted = 1; + } + + /** + * @return time between last call of timings_*() and now + */ + double update() { + double diffc, difft; + clock_t lc = lastclock; // measure CPU in 1e-6s + time_t lt = lasttime; // measure time in s + + assert( + isstarted == 1 && + "timings_started() must be called before using timings... functions"); + + lastclock = + clock(); // measures at most 2147 seconds, where 1s = 1e6 CLOCKS_PER_SEC + lasttime = time(NULL); + diffc = (double)(lastclock - lc) / + CLOCKS_PER_SEC; // is presumably in [-21??, 21??] + difft = difftime(lasttime, lt); // is presumably an integer + lastdiff = difft; // on the "save" side + // use diffc clock measurement if appropriate + if (diffc > 0 && difft < 1000) + lastdiff = diffc; + assert(lastdiff >= 0 && "BUG in time measurement"); + totaltime += lastdiff; + totaltotaltime += lastdiff; + if (istic) { + tictoczwischensumme += lastdiff; + tictoctime += lastdiff; + } + return lastdiff; + } + + void tic() { + assert(!istic && "Timingic called twice without toc"); + update(); + istic = 1; + } + + double toc() { + assert(istic && "Timingoc called without tic"); + update(); + lasttictoctime = tictoczwischensumme; + tictoczwischensumme = 0; + istic = 0; + return lasttictoctime; + } +}; diff --git a/include/cma-es/utils.h b/include/cma-es/utils.h new file mode 100644 index 0000000..94d349b --- /dev/null +++ b/include/cma-es/utils.h @@ -0,0 +1,41 @@ +/** + * @file utils.h + * Contains some utility functions. + */ + +#pragma once + +#include +#include +#include + +template T square(T d) { return d * d; } + +template T maxElement(const T *rgd, int len) { + return *std::max_element(rgd, rgd + len); +} + +template T minElement(const T *rgd, int 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 int minIndex(const T *rgd, int len) { + return 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); + 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); + } else + return T(0); +} diff --git a/include/lexi.h b/include/lexi.h new file mode 100644 index 0000000..2536141 --- /dev/null +++ b/include/lexi.h @@ -0,0 +1,51 @@ +#include + +template class LexiProduct { + using Product = std::vector; + +public: + LexiProduct() : LexiProduct(0) {} + LexiProduct(size_t size) : prd(size) {} + LexiProduct &operator=(const Product &prd2) { + prd = prd2; + return *this; + } + LexiProduct &operator=(Product &&prd2) { + prd = std::move(prd2); + return *this; + } + + bool operator==(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) == 0; + } + bool operator!=(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) != 0; + } + bool operator<(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) < 0; + } + bool operator<=(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) <= 0; + } + bool operator>(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) > 0; + } + bool operator>=(const LexiProduct &lexi2) { + return lexicmp(*this, lexi2) >= 0; + } + +private: + Product prd; + + int lexicmp(const LexiProduct &lhs, const LexiProduct &rhs) { + const int lencmp = 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]; + if (d != 0) { + return d; + } + } + return lencmp; + } +}; diff --git a/src/conversion.cpp b/src/conversion.cpp index a493f2c..bbb54e8 100644 --- a/src/conversion.cpp +++ b/src/conversion.cpp @@ -12,9 +12,12 @@ RGB XYZtoRGB(const XYZ &xyz) { double g = xyz.x * A21 + xyz.y * A22 + xyz.z * A23; double b = xyz.x * A31 + xyz.y * A32 + xyz.z * A33; - r = ((abs(r) > 0.0031308) ? sgn(r) * (1.055 * pow(abs(r), 1 / 2.4) - 0.055) : (12.92 * r)); - g = ((abs(g) > 0.0031308) ? sgn(g) * (1.055 * pow(abs(g), 1 / 2.4) - 0.055) : (12.92 * g)); - b = ((abs(b) > 0.0031308) ? sgn(b) * (1.055 * pow(abs(b), 1 / 2.4) - 0.055) : (12.92 * b)); + r = ((abs(r) > 0.0031308) ? sgn(r) * (1.055 * pow(abs(r), 1 / 2.4) - 0.055) + : (12.92 * r)); + g = ((abs(g) > 0.0031308) ? sgn(g) * (1.055 * pow(abs(g), 1 / 2.4) - 0.055) + : (12.92 * g)); + b = ((abs(b) > 0.0031308) ? sgn(b) * (1.055 * pow(abs(b), 1 / 2.4) - 0.055) + : (12.92 * b)); return RGB{r, g, b}; } diff --git a/test/example1.cpp b/test/example1.cpp new file mode 100644 index 0000000..a1253c8 --- /dev/null +++ b/test/example1.cpp @@ -0,0 +1,81 @@ +/** + * @file example1.cpp + * Very short example source code. The purpose of the example codes is to be + * edited/extended. + */ + +#include +#include +#include + +/** + * The objective (fitness) function to be minized. "cigtab" function. + */ +double fitfun(double const *x, int N) { + double sum = 1e4 * x[0] * x[0] + 1e-4 * x[1] * x[1]; + for (int i = 2; i < N; ++i) + sum += x[i] * x[i]; + return sum; +} + +/** + * The optimization loop. + */ +int main(int, char **) { + CMAES evo; + double *arFunvals, *const *pop, *xfinal; + + // Initialize everything + const int dim = 22; + double xstart[dim]; + for (int i = 0; i < dim; i++) + xstart[i] = 0.5; + double stddev[dim]; + for (int i = 0; i < dim; i++) + stddev[i] = 0.3; + Parameters parameters; + // TODO Adjust parameters here + parameters.init(dim, xstart, stddev); + arFunvals = evo.init(parameters); + + std::cout << evo.sayHello() << std::endl; + + // Iterate until stop criterion holds + while (!evo.testForTermination()) { + // Generate lambda new search points, sample population + pop = evo.samplePopulation(); // Do not change content of pop + + /* Here you may resample each solution point pop[i] until it + becomes feasible, e.g. for box constraints (variable + boundaries). function is_feasible(...) needs to be + user-defined. + Assumptions: the feasible domain is convex, the optimum is + not on (or very close to) the domain boundary, initialX is + feasible and initialStandardDeviations are sufficiently small + to prevent quasi-infinite looping. + */ + /* for (i = 0; i < evo.get(CMAES::PopSize); ++i) + while (!is_feasible(pop[i])) + evo.reSampleSingle(i); + */ + + // 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)); + + // update the search distribution used for sampleDistribution() + evo.updateDistribution(arFunvals); + } + std::cout << "Stop:" << std::endl << evo.getStopMessage(); + 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 + + // do something with final solution and finally release memory + delete[] xfinal; + + return 0; +} diff --git a/test/example2.cpp b/test/example2.cpp new file mode 100644 index 0000000..6ad1915 --- /dev/null +++ b/test/example2.cpp @@ -0,0 +1,483 @@ +/** + * @file example2.cpp + * Implements additional restarts with increasing population + * (Auger & Hansen 2005). + */ + +#include +#include +#include +#include +#include +#include +#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 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 + double incpopsize = 2; + double *x; + + // Put together objective functions + rgpFun[0] = f_sphere; + rgpFun[1] = f_elli; + rgpFun[2] = f_cigar; + rgpFun[3] = f_cigtab; + rgpFun[4] = f_tablet; + rgpFun[5] = f_rosenbrock; + rgpFun[6] = f_parabR; + rgpFun[7] = f_sharpR; + rgpFun[8] = f_diffpow; + rgpFun[9] = f_kugelmin1; + rgpFun[10] = f_ellinumtest; + rgpFun[11] = f_elli100; + rgpFun[18] = f_gleichsys5; + rgpFun[19] = f_rand; + rgpFun[20] = f_constant; + rgpFun[21] = f_stepsphere; + rgpFun[22] = f_ellirot; + rgpFun[23] = f_diffpowrot; + + int nb = 5; + int nbrestarts = 0; + + // Optimize function + x = optimize(rgpFun[nb], nbrestarts, incpopsize); + + // here we could utilize the solution x, and finally free memory + delete[] x; + + return 0; +} + +/** + * 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 incpopsize) { + CMAES evo; // the optimizer + double *const *pop; // sampled population + double *fitvals; // objective function values of sampled population + double fbestever = 0, *xbestever = NULL; // store best solution + double fmean; + int irun, + lambda = 0, // offspring population size, 0 invokes default + countevals = 0; // used to set for restarts + + for (irun = 0; irun < nrestarts + 1; ++irun) { + /* Parameters can be set in two ways. Here as input parameter to evo.init() + and as value read from signals.par by calling evo.readSignals + explicitely. + */ + const int dim = 22; + double xstart[dim]; + for (int i = 0; i < dim; i++) + xstart[i] = 0.5; + double stddev[dim]; + for (int i = 0; i < dim; i++) + stddev[i] = 0.3; + Parameters parameters; + // You can resume a previous run by specifying a file that contains the + // resume data: + // parameters.resumefile = "resumeevo2.dat"; + parameters.logWarnings = true; // warnings will be printed on std::cerr + parameters.stopTolX = 1e-11; + parameters.updateCmode.maxtime = 1.0; + parameters.lambda = lambda; + parameters.init(dim, xstart, stddev); + + fitvals = evo.init(parameters); // allocs fitvals + std::cout << evo.sayHello() << std::endl; + evo.countevals = countevals; // a hack, effects the output and termination + + while (!evo.testForTermination()) { + // Generate population of new candidate solutions + pop = evo.samplePopulation(); // do not change content of pop + + /* Here optionally handle constraints etc. on pop. You may + * call evo.reSampleSingle(i) to resample the i-th + * vector pop[i], see below. Do not change pop in any other + * way. You may also copy and modify (repair) pop[i] only + * for the evaluation of the fitness function and consider + * adding a penalty depending on the size of the + * modification. + */ + + // Compute fitness value for each candidate solution + 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 + is_feasible() needs to be user-defined. + Assumptions: the feasible domain is convex, the optimum + is not on (or very close to) the domain boundary, + initialX is feasible (or in case typicalX +- + 2*initialStandardDeviations is feasible) and + initialStandardDeviations is (are) sufficiently small to prevent + quasi-infinite looping. + */ + /* while (!is_feasible(pop[i])) + evo.reSampleSingle(i); + */ + fitvals[i] = (*pFun)(pop[i]); + } + + // update search distribution + evo.updateDistribution(fitvals); + + fflush(stdout); + } + + 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 + << " (axis-ratio=" + << 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"); + + // keep best ever solution + if (irun == 0 || evo.get(CMAES::FBestEver) < fbestever) { + fbestever = evo.get(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) { + fbestever = fmean; + xbestever = evo.getInto(CMAES::XMean, xbestever); + } + + // abandon restarts if target fitness value was achieved or MaxFunEvals + // reached + std::string stop = evo.getStopMessage(); + if (strncmp(stop.c_str(), "Fitness", 7) == 0 || + strncmp(stop.c_str(), "MaxFunEvals", 11) == 0) + break; + if (strncmp(stop.c_str(), "Manual", 6) == 0) { + printf("Press RETURN to start next run\n"); + fflush(stdout); + getchar(); + } + } + + return xbestever; // was dynamically allocated, should be deleted in the end +} + +double f_rand(double const *) { + double d = (double)rand() / RAND_MAX; + while (d == 0.) + d = (double)rand() / RAND_MAX; + return d; +} + +double f_constant(double const *) { return 1; } + +static double SQR(double d) { return d * d; } + +double f_stepsphere(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + for (i = 0; i < DIM; ++i) + sum += floor(x[i] * x[i]); + return sum; +} + +double f_sphere(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + for (i = 0; i < DIM; ++i) + sum += x[i] * x[i]; + return sum; +} + +double f_cigar(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + for (i = 1; i < DIM; ++i) + sum += x[i] * x[i]; + sum *= 1e6; + sum += x[0] * x[0]; + return sum; +} + +double f_cigtab(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + sum = x[0] * x[0] + 1e8 * x[DIM - 1] * x[DIM - 1]; + for (i = 1; i < DIM - 1; ++i) + sum += 1e4 * x[i] * x[i]; + return sum; +} + +double f_tablet(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + sum = 1e6 * x[0] * x[0]; + for (i = 1; i < DIM; ++i) + sum += x[i] * x[i]; + return sum; +} + +/* a hack, memory is never released */ +double **OrthogonalBasis(int DIM) { + static int b_dim; + static double **b; + double sp; + int i, j, k; + Random R(2); + + if (b_dim != 0) { // Initialization was done + if (b_dim != DIM) { + printf("function OrthogonalBasis cannot change dimensionality in file " + "example2.cpp"); + exit(0); + } + return b; + } + + // Otherwise initialize basis b + + // allocate b + b = (double **)calloc((unsigned)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)); + if (!b[i]) { + printf("calloc failed in function Orthogonalbasis in file example2.cpp"); + exit(0); + } + } + b_dim = DIM; + + // generate orthogonal basis + for (i = 0; i < DIM; ++i) { + // sample components gaussian + for (j = 0; j < DIM; ++j) + b[i][j] = R.gauss(); + // substract projection of previous vectors + for (j = i - 1; j >= 0; --j) { + for (sp = 0., k = 0; k < DIM; ++k) + sp += b[i][k] * b[j][k]; // scalar product + for (k = 0; k < DIM; ++k) + b[i][k] -= sp * b[j][k]; // substract + } + // normalize + for (sp = 0., k = 0; k < DIM; ++k) + sp += b[i][k] * b[i][k]; // squared norm + for (k = 0; k < DIM; ++k) + b[i][k] /= sqrt(sp); + } + + return b; +} + +double f_ellirot(double const *x) { + int i, k; + double sum = 0., y; + int DIM = (int)(x[-1]); + double **B = OrthogonalBasis(DIM); + + if (DIM == 1) + return x[0] * x[0]; + for (i = 0; i < DIM; ++i) { + for (y = 0., k = 0; k < DIM; ++k) + y += B[i][k] * x[k]; + sum += exp(log(1e6) * 2. * (double)(i) / (DIM - 1)) * y * y; + } + return sum; +} + +double f_elli(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + if (DIM == 1) + return x[0] * x[0]; + for (i = 0; i < DIM; ++i) + sum += exp(log(1000.) * 2. * (double)(i) / (DIM - 1)) * x[i] * x[i]; + return sum; +} + +double f_elli100(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + if (DIM == 1) + return x[0] * x[0]; + for (i = 0; i < DIM; ++i) + sum += exp(log(100.) * 2. * (double)(i) / (DIM - 1)) * x[i] * x[i]; + return sum; +} + +double f_diffpow(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + if (DIM == 1) + return x[0] * x[0]; + for (i = 0; i < DIM; ++i) + sum += pow(fabs(x[i]), 2. + 10 * (double)(i) / (DIM - 1)); + return sum; +} + +double f_diffpowrot(double const *x) { + int i, k; + double sum = 0., y; + int DIM = (int)(x[-1]); + double **B = OrthogonalBasis(DIM); + + if (DIM == 1) + return x[0] * x[0]; + for (i = 0; i < DIM; ++i) { + for (y = 0., k = 0; k < DIM; ++k) + y += B[i][k] * x[k]; + sum += pow(fabs(y), 2. + 10 * (double)(i) / (DIM - 1)); + } + return sum; +} + +double f_kugelmin1(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + + for (i = 1; i < DIM; ++i) + sum += x[i] * x[i]; + return sum; +} + +/** + * Rosenbrock's Function, generalized. + */ +double f_rosenbrock(double const *x) { + double qualitaet; + int i; + int DIM = (int)(x[-1]); + qualitaet = 0.0; + + for (i = DIM - 2; i >= 0; --i) + qualitaet += 100. * SQR(SQR(x[i]) - x[i + 1]) + SQR(1. - x[i]); + return (qualitaet); +} + +double f_parabR(double const *x) { + int i; + double sum = 0.; + int DIM = (int)(x[-1]); + 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 sum = 0.; + int DIM = (int)(x[-1]); + 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 sum = 0.; + int DIM = (int)(x[-1]); + static double maxVerhaeltnis = 0.; + if (maxVerhaeltnis == 0.) { + for (maxVerhaeltnis = 1.; + maxVerhaeltnis < 1e99 && maxVerhaeltnis < 2. * maxVerhaeltnis; + maxVerhaeltnis *= 2.) + if (maxVerhaeltnis == maxVerhaeltnis + 1.) + break; + maxVerhaeltnis *= 10.; + maxVerhaeltnis = sqrt(maxVerhaeltnis); + } + if (DIM < 3) + return x[0] * x[0]; + for (i = 1; i < DIM; ++i) + sum += exp(log(maxVerhaeltnis) * 2. * (double)(i - 1) / (DIM - 2)) * x[i] * + x[i]; + return sum; +} + +/** + * 5-dimensional system of equations (Juergen Bremer). + * For each row the following equation holds: + * + * c_1*x[1] + c_2*x[2] + c_3*x[3] + c_4*x[4] + c_5*x[5] + c_0 = 0 + * + * This is the reason why the square of the left side is in the quality + * function. + */ +double f_gleichsys5(double const *x) { + double qualitaet = 0.0; + + static double koeff[5][6] = {/* c_1, c_2, c_3, c_4, c_5, c_0 */ + {4, 191, 27, 199, 21, 172}, + {191, 10883, 1413, 5402, 684, -8622}, + {27, 1413, 191, 1032, 118, -94}, + {199, 5402, 1032, 29203, 2331, 78172}, + {21, 684, 118, 2331, 199, 5648}}; + int i, j; + double sum; + + for (i = 0; i < 5; ++i) { + sum = koeff[i][5]; + for (j = 0; j < 5; ++j) { + sum += koeff[i][j] * x[j]; + } + qualitaet += sum * sum; + } + return qualitaet; +} diff --git a/test/testlexi.cpp b/test/testlexi.cpp new file mode 100644 index 0000000..a5577c6 --- /dev/null +++ b/test/testlexi.cpp @@ -0,0 +1,26 @@ +#include +#include +#include + +int main() { + LexiProduct p1, p2, p2a, p3, p4; + p1 = std::vector{11, 3, 4}; + p2 = 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); + + std::cout << "testlexi: All tests passed" << std::endl; +}; \ No newline at end of file