diff --git a/HaskellNN.cabal b/HaskellNN.cabal new file mode 100644 --- /dev/null +++ b/HaskellNN.cabal @@ -0,0 +1,62 @@ +Name: HaskellNN +Version: 0.1 +License: GPL +License-file: LICENSE +Author: Kiet Lam +Maintainer: Kiet Lam +Synopsis: High Performance Neural Network in Haskell +Description: High Performance Neural Network in Haskell + . + Provides fast training algorithms using + hmatrix's bindings to GSL and custom bindings + to the liblbfgs C-library + . + Supported training algorithms: Gradient Descent, Conjugate Gradient, BFGS, LBFGS + . + - Users should focus on "AI.Model" for most usages (classification / regression) + . + - Other modules are provided for user expansion if needed + . + Go to for examples and tests for usage + +Category: AI + +Build-type: Simple +Cabal-version: >= 1.6 + +Library + + Build-depends: base >= 4 && < 5, + hmatrix >= 0.13.0.0, + random + + Extensions: ForeignFunctionInterface + + hs-source-dirs: src + Exposed-modules: AI.Calculation, + AI.Calculation.Activation, + AI.Calculation.Cost, + AI.Calculation.Gradients, + AI.Calculation.NetworkOutput, + AI.Signatures, + AI.Model, + AI.Model.Classification, + AI.Model.General, + AI.Model.GenericModel, + AI.Training, + AI.Network + + Other-modules: AI.Training.Internal, + AI.Training.Internal.LBFGSAux + + + ghc-prof-options: -prof -auto-all + Include-Dirs: cbits + C-sources: src/AI/Training/Internal/lbfgs_aux.c, + cbits/lbfgs.c + Includes: lbfgs.h + + +source-repository head + type: git + location: https://github.com/ktklam9/HaskellNN diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cbits/lbfgs.c b/cbits/lbfgs.c new file mode 100644 --- /dev/null +++ b/cbits/lbfgs.c @@ -0,0 +1,1371 @@ +/* + * Limited memory BFGS (L-BFGS). + * + * Copyright (c) 1990, Jorge Nocedal + * Copyright (c) 2007-2010 Naoaki Okazaki + * All rights reserved. + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +/* $Id$ */ + +/* +This library is a C port of the FORTRAN implementation of Limited-memory +Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method written by Jorge Nocedal. +The original FORTRAN source code is available at: +http://www.ece.northwestern.edu/~nocedal/lbfgs.html + +The L-BFGS algorithm is described in: + - Jorge Nocedal. + Updating Quasi-Newton Matrices with Limited Storage. + Mathematics of Computation, Vol. 35, No. 151, pp. 773--782, 1980. + - Dong C. Liu and Jorge Nocedal. + On the limited memory BFGS method for large scale optimization. + Mathematical Programming B, Vol. 45, No. 3, pp. 503-528, 1989. + +The line search algorithms used in this implementation are described in: + - John E. Dennis and Robert B. Schnabel. + Numerical Methods for Unconstrained Optimization and Nonlinear + Equations, Englewood Cliffs, 1983. + - Jorge J. More and David J. Thuente. + Line search algorithm with guaranteed sufficient decrease. + ACM Transactions on Mathematical Software (TOMS), Vol. 20, No. 3, + pp. 286-307, 1994. + +This library also implements Orthant-Wise Limited-memory Quasi-Newton (OWL-QN) +method presented in: + - Galen Andrew and Jianfeng Gao. + Scalable training of L1-regularized log-linear models. + In Proceedings of the 24th International Conference on Machine + Learning (ICML 2007), pp. 33-40, 2007. + +I would like to thank the original author, Jorge Nocedal, who has been +distributing the effieicnt and explanatory implementation in an open source +licence. +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif/*HAVE_CONFIG_H*/ + +#include +#include +#include +#include + +#include + +#ifdef _MSC_VER +#define inline __inline +#endif/*_MSC_VER*/ + +#if defined(USE_SSE) && defined(__SSE2__) && LBFGS_FLOAT == 64 +/* Use SSE2 optimization for 64bit double precision. */ +#include "arithmetic_sse_double.h" + +#elif defined(USE_SSE) && defined(__SSE__) && LBFGS_FLOAT == 32 +/* Use SSE optimization for 32bit float precision. */ +#include "arithmetic_sse_float.h" + +#else +/* No CPU specific optimization. */ +#include "arithmetic_ansi.h" + +#endif + +#define min2(a, b) ((a) <= (b) ? (a) : (b)) +#define max2(a, b) ((a) >= (b) ? (a) : (b)) +#define max3(a, b, c) max2(max2((a), (b)), (c)); + +struct tag_callback_data { + int n; + void *instance; + lbfgs_evaluate_t proc_evaluate; + lbfgs_progress_t proc_progress; +}; +typedef struct tag_callback_data callback_data_t; + +struct tag_iteration_data { + lbfgsfloatval_t alpha; + lbfgsfloatval_t *s; /* [n] */ + lbfgsfloatval_t *y; /* [n] */ + lbfgsfloatval_t ys; /* vecdot(y, s) */ +}; +typedef struct tag_iteration_data iteration_data_t; + +static const lbfgs_parameter_t _defparam = { + 6, 1e-5, 0, 1e-5, + 0, LBFGS_LINESEARCH_DEFAULT, 40, + 1e-20, 1e20, 1e-4, 0.9, 0.9, 1.0e-16, + 0.0, 0, -1, +}; + +/* Forward function declarations. */ + +typedef int (*line_search_proc)( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wa, + callback_data_t *cd, + const lbfgs_parameter_t *param + ); + +static int line_search_backtracking( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wa, + callback_data_t *cd, + const lbfgs_parameter_t *param + ); + +static int line_search_backtracking_owlqn( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wp, + callback_data_t *cd, + const lbfgs_parameter_t *param + ); + +static int line_search_morethuente( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wa, + callback_data_t *cd, + const lbfgs_parameter_t *param + ); + +static int update_trial_interval( + lbfgsfloatval_t *x, + lbfgsfloatval_t *fx, + lbfgsfloatval_t *dx, + lbfgsfloatval_t *y, + lbfgsfloatval_t *fy, + lbfgsfloatval_t *dy, + lbfgsfloatval_t *t, + lbfgsfloatval_t *ft, + lbfgsfloatval_t *dt, + const lbfgsfloatval_t tmin, + const lbfgsfloatval_t tmax, + int *brackt + ); + +static lbfgsfloatval_t owlqn_x1norm( + const lbfgsfloatval_t* x, + const int start, + const int n + ); + +static void owlqn_pseudo_gradient( + lbfgsfloatval_t* pg, + const lbfgsfloatval_t* x, + const lbfgsfloatval_t* g, + const int n, + const lbfgsfloatval_t c, + const int start, + const int end + ); + +static void owlqn_project( + lbfgsfloatval_t* d, + const lbfgsfloatval_t* sign, + const int start, + const int end + ); + + +#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__)) +static int round_out_variables(int n) +{ + n += 7; + n /= 8; + n *= 8; + return n; +} +#endif/*defined(USE_SSE)*/ + +lbfgsfloatval_t* lbfgs_malloc(int n) +{ +#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__)) + n = round_out_variables(n); +#endif/*defined(USE_SSE)*/ + return (lbfgsfloatval_t*)vecalloc(sizeof(lbfgsfloatval_t) * n); +} + +void lbfgs_free(lbfgsfloatval_t *x) +{ + vecfree(x); +} + +void lbfgs_parameter_init(lbfgs_parameter_t *param) +{ + memcpy(param, &_defparam, sizeof(*param)); +} + +int lbfgs( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *ptr_fx, + lbfgs_evaluate_t proc_evaluate, + lbfgs_progress_t proc_progress, + void *instance, + lbfgs_parameter_t *_param + ) +{ + int ret; + int i, j, k, ls, end, bound; + lbfgsfloatval_t step; + + /* Constant parameters and their default values. */ + lbfgs_parameter_t param = (_param != NULL) ? (*_param) : _defparam; + const int m = param.m; + + lbfgsfloatval_t *xp = NULL; + lbfgsfloatval_t *g = NULL, *gp = NULL, *pg = NULL; + lbfgsfloatval_t *d = NULL, *w = NULL, *pf = NULL; + iteration_data_t *lm = NULL, *it = NULL; + lbfgsfloatval_t ys, yy; + lbfgsfloatval_t xnorm, gnorm, beta; + lbfgsfloatval_t fx = 0.; + lbfgsfloatval_t rate = 0.; + line_search_proc linesearch = line_search_morethuente; + + /* Construct a callback data. */ + callback_data_t cd; + cd.n = n; + cd.instance = instance; + cd.proc_evaluate = proc_evaluate; + cd.proc_progress = proc_progress; + +#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__)) + /* Round out the number of variables. */ + n = round_out_variables(n); +#endif/*defined(USE_SSE)*/ + + /* Check the input parameters for errors. */ + if (n <= 0) { + return LBFGSERR_INVALID_N; + } +#if defined(USE_SSE) && (defined(__SSE__) || defined(__SSE2__)) + if (n % 8 != 0) { + return LBFGSERR_INVALID_N_SSE; + } + if ((uintptr_t)(const void*)x % 16 != 0) { + return LBFGSERR_INVALID_X_SSE; + } +#endif/*defined(USE_SSE)*/ + if (param.epsilon < 0.) { + return LBFGSERR_INVALID_EPSILON; + } + if (param.past < 0) { + return LBFGSERR_INVALID_TESTPERIOD; + } + if (param.delta < 0.) { + return LBFGSERR_INVALID_DELTA; + } + if (param.min_step < 0.) { + return LBFGSERR_INVALID_MINSTEP; + } + if (param.max_step < param.min_step) { + return LBFGSERR_INVALID_MAXSTEP; + } + if (param.ftol < 0.) { + return LBFGSERR_INVALID_FTOL; + } + if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE || + param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE) { + if (param.wolfe <= param.ftol || 1. <= param.wolfe) { + return LBFGSERR_INVALID_WOLFE; + } + } + if (param.gtol < 0.) { + return LBFGSERR_INVALID_GTOL; + } + if (param.xtol < 0.) { + return LBFGSERR_INVALID_XTOL; + } + if (param.max_linesearch <= 0) { + return LBFGSERR_INVALID_MAXLINESEARCH; + } + if (param.orthantwise_c < 0.) { + return LBFGSERR_INVALID_ORTHANTWISE; + } + if (param.orthantwise_start < 0 || n < param.orthantwise_start) { + return LBFGSERR_INVALID_ORTHANTWISE_START; + } + if (param.orthantwise_end < 0) { + param.orthantwise_end = n; + } + if (n < param.orthantwise_end) { + return LBFGSERR_INVALID_ORTHANTWISE_END; + } + if (param.orthantwise_c != 0.) { + switch (param.linesearch) { + case LBFGS_LINESEARCH_BACKTRACKING: + linesearch = line_search_backtracking_owlqn; + break; + default: + /* Only the backtracking method is available. */ + return LBFGSERR_INVALID_LINESEARCH; + } + } else { + switch (param.linesearch) { + case LBFGS_LINESEARCH_MORETHUENTE: + linesearch = line_search_morethuente; + break; + case LBFGS_LINESEARCH_BACKTRACKING_ARMIJO: + case LBFGS_LINESEARCH_BACKTRACKING_WOLFE: + case LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE: + linesearch = line_search_backtracking; + break; + default: + return LBFGSERR_INVALID_LINESEARCH; + } + } + + /* Allocate working space. */ + xp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + g = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + gp = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + d = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + w = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + if (xp == NULL || g == NULL || gp == NULL || d == NULL || w == NULL) { + ret = LBFGSERR_OUTOFMEMORY; + goto lbfgs_exit; + } + + if (param.orthantwise_c != 0.) { + /* Allocate working space for OW-LQN. */ + pg = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + if (pg == NULL) { + ret = LBFGSERR_OUTOFMEMORY; + goto lbfgs_exit; + } + } + + /* Allocate limited memory storage. */ + lm = (iteration_data_t*)vecalloc(m * sizeof(iteration_data_t)); + if (lm == NULL) { + ret = LBFGSERR_OUTOFMEMORY; + goto lbfgs_exit; + } + + /* Initialize the limited memory. */ + for (i = 0;i < m;++i) { + it = &lm[i]; + it->alpha = 0; + it->ys = 0; + it->s = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + it->y = (lbfgsfloatval_t*)vecalloc(n * sizeof(lbfgsfloatval_t)); + if (it->s == NULL || it->y == NULL) { + ret = LBFGSERR_OUTOFMEMORY; + goto lbfgs_exit; + } + } + + /* Allocate an array for storing previous values of the objective function. */ + if (0 < param.past) { + pf = (lbfgsfloatval_t*)vecalloc(param.past * sizeof(lbfgsfloatval_t)); + } + + /* Evaluate the function value and its gradient. */ + fx = cd.proc_evaluate(cd.instance, x, g, cd.n, 0); + if (0. != param.orthantwise_c) { + /* Compute the L1 norm of the variable and add it to the object value. */ + xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end); + fx += xnorm * param.orthantwise_c; + owlqn_pseudo_gradient( + pg, x, g, n, + param.orthantwise_c, param.orthantwise_start, param.orthantwise_end + ); + } + + /* Store the initial value of the objective function. */ + if (pf != NULL) { + pf[0] = fx; + } + + /* + Compute the direction; + we assume the initial hessian matrix H_0 as the identity matrix. + */ + if (param.orthantwise_c == 0.) { + vecncpy(d, g, n); + } else { + vecncpy(d, pg, n); + } + + /* + Make sure that the initial variables are not a minimizer. + */ + vec2norm(&xnorm, x, n); + if (param.orthantwise_c == 0.) { + vec2norm(&gnorm, g, n); + } else { + vec2norm(&gnorm, pg, n); + } + if (xnorm < 1.0) xnorm = 1.0; + if (gnorm / xnorm <= param.epsilon) { + ret = LBFGS_ALREADY_MINIMIZED; + goto lbfgs_exit; + } + + /* Compute the initial step: + step = 1.0 / sqrt(vecdot(d, d, n)) + */ + vec2norminv(&step, d, n); + + k = 1; + end = 0; + for (;;) { + /* Store the current position and gradient vectors. */ + veccpy(xp, x, n); + veccpy(gp, g, n); + + /* Search for an optimal step. */ + if (param.orthantwise_c == 0.) { + ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, ¶m); + } else { + ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, ¶m); + owlqn_pseudo_gradient( + pg, x, g, n, + param.orthantwise_c, param.orthantwise_start, param.orthantwise_end + ); + } + if (ls < 0) { + /* Revert to the previous point. */ + veccpy(x, xp, n); + veccpy(g, gp, n); + ret = ls; + goto lbfgs_exit; + } + + /* Compute x and g norms. */ + vec2norm(&xnorm, x, n); + if (param.orthantwise_c == 0.) { + vec2norm(&gnorm, g, n); + } else { + vec2norm(&gnorm, pg, n); + } + + /* Report the progress. */ + if (cd.proc_progress) { + if ((ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls))) { + goto lbfgs_exit; + } + } + + /* + Convergence test. + The criterion is given by the following formula: + |g(x)| / \max(1, |x|) < \epsilon + */ + if (xnorm < 1.0) xnorm = 1.0; + if (gnorm / xnorm <= param.epsilon) { + /* Convergence. */ + ret = LBFGS_SUCCESS; + break; + } + + /* + Test for stopping criterion. + The criterion is given by the following formula: + (f(past_x) - f(x)) / f(x) < \delta + */ + if (pf != NULL) { + /* We don't test the stopping criterion while k < past. */ + if (param.past <= k) { + /* Compute the relative improvement from the past. */ + rate = (pf[k % param.past] - fx) / fx; + + /* The stopping criterion. */ + if (rate < param.delta) { + ret = LBFGS_STOP; + break; + } + } + + /* Store the current value of the objective function. */ + pf[k % param.past] = fx; + } + + if (param.max_iterations != 0 && param.max_iterations < k+1) { + /* Maximum number of iterations. */ + ret = LBFGSERR_MAXIMUMITERATION; + break; + } + + /* + Update vectors s and y: + s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}. + y_{k+1} = g_{k+1} - g_{k}. + */ + it = &lm[end]; + vecdiff(it->s, x, xp, n); + vecdiff(it->y, g, gp, n); + + /* + Compute scalars ys and yy: + ys = y^t \cdot s = 1 / \rho. + yy = y^t \cdot y. + Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor). + */ + vecdot(&ys, it->y, it->s, n); + vecdot(&yy, it->y, it->y, n); + it->ys = ys; + + /* + Recursive formula to compute dir = -(H \cdot g). + This is described in page 779 of: + Jorge Nocedal. + Updating Quasi-Newton Matrices with Limited Storage. + Mathematics of Computation, Vol. 35, No. 151, + pp. 773--782, 1980. + */ + bound = (m <= k) ? m : k; + ++k; + end = (end + 1) % m; + + /* Compute the steepest direction. */ + if (param.orthantwise_c == 0.) { + /* Compute the negative of gradients. */ + vecncpy(d, g, n); + } else { + vecncpy(d, pg, n); + } + + j = end; + for (i = 0;i < bound;++i) { + j = (j + m - 1) % m; /* if (--j == -1) j = m-1; */ + it = &lm[j]; + /* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */ + vecdot(&it->alpha, it->s, d, n); + it->alpha /= it->ys; + /* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */ + vecadd(d, it->y, -it->alpha, n); + } + + vecscale(d, ys / yy, n); + + for (i = 0;i < bound;++i) { + it = &lm[j]; + /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */ + vecdot(&beta, it->y, d, n); + beta /= it->ys; + /* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */ + vecadd(d, it->s, it->alpha - beta, n); + j = (j + 1) % m; /* if (++j == m) j = 0; */ + } + + /* + Constrain the search direction for orthant-wise updates. + */ + if (param.orthantwise_c != 0.) { + for (i = param.orthantwise_start;i < param.orthantwise_end;++i) { + if (d[i] * pg[i] >= 0) { + d[i] = 0; + } + } + } + + /* + Now the search direction d is ready. We try step = 1 first. + */ + step = 1.0; + } + +lbfgs_exit: + /* Return the final value of the objective function. */ + if (ptr_fx != NULL) { + *ptr_fx = fx; + } + + vecfree(pf); + + /* Free memory blocks used by this function. */ + if (lm != NULL) { + for (i = 0;i < m;++i) { + vecfree(lm[i].s); + vecfree(lm[i].y); + } + vecfree(lm); + } + vecfree(pg); + vecfree(w); + vecfree(d); + vecfree(gp); + vecfree(g); + vecfree(xp); + + return ret; +} + + + +static int line_search_backtracking( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wp, + callback_data_t *cd, + const lbfgs_parameter_t *param + ) +{ + int count = 0; + lbfgsfloatval_t width, dg; + lbfgsfloatval_t finit, dginit = 0., dgtest; + const lbfgsfloatval_t dec = 0.5, inc = 2.1; + + /* Check the input parameters for errors. */ + if (*stp <= 0.) { + return LBFGSERR_INVALIDPARAMETERS; + } + + /* Compute the initial gradient in the search direction. */ + vecdot(&dginit, g, s, n); + + /* Make sure that s points to a descent direction. */ + if (0 < dginit) { + return LBFGSERR_INCREASEGRADIENT; + } + + /* The initial value of the objective function. */ + finit = *f; + dgtest = param->ftol * dginit; + + for (;;) { + veccpy(x, xp, n); + vecadd(x, s, *stp, n); + + /* Evaluate the function and gradient values. */ + *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp); + + ++count; + + if (*f > finit + *stp * dgtest) { + width = dec; + } else { + /* The sufficient decrease condition (Armijo condition). */ + if (param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_ARMIJO) { + /* Exit with the Armijo condition. */ + return count; + } + + /* Check the Wolfe condition. */ + vecdot(&dg, g, s, n); + if (dg < param->wolfe * dginit) { + width = inc; + } else { + if(param->linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE) { + /* Exit with the regular Wolfe condition. */ + return count; + } + + /* Check the strong Wolfe condition. */ + if(dg > -param->wolfe * dginit) { + width = dec; + } else { + /* Exit with the strong Wolfe condition. */ + return count; + } + } + } + + if (*stp < param->min_step) { + /* The step is the minimum value. */ + return LBFGSERR_MINIMUMSTEP; + } + if (*stp > param->max_step) { + /* The step is the maximum value. */ + return LBFGSERR_MAXIMUMSTEP; + } + if (param->max_linesearch <= count) { + /* Maximum number of iteration. */ + return LBFGSERR_MAXIMUMLINESEARCH; + } + + (*stp) *= width; + } +} + + + +static int line_search_backtracking_owlqn( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wp, + callback_data_t *cd, + const lbfgs_parameter_t *param + ) +{ + int i, count = 0; + lbfgsfloatval_t width = 0.5, norm = 0.; + lbfgsfloatval_t finit = *f, dgtest; + + /* Check the input parameters for errors. */ + if (*stp <= 0.) { + return LBFGSERR_INVALIDPARAMETERS; + } + + /* Choose the orthant for the new point. */ + for (i = 0;i < n;++i) { + wp[i] = (xp[i] == 0.) ? -gp[i] : xp[i]; + } + + for (;;) { + /* Update the current point. */ + veccpy(x, xp, n); + vecadd(x, s, *stp, n); + + /* The current point is projected onto the orthant. */ + owlqn_project(x, wp, param->orthantwise_start, param->orthantwise_end); + + /* Evaluate the function and gradient values. */ + *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp); + + /* Compute the L1 norm of the variables and add it to the object value. */ + norm = owlqn_x1norm(x, param->orthantwise_start, param->orthantwise_end); + *f += norm * param->orthantwise_c; + + ++count; + + dgtest = 0.; + for (i = 0;i < n;++i) { + dgtest += (x[i] - xp[i]) * gp[i]; + } + + if (*f <= finit + param->ftol * dgtest) { + /* The sufficient decrease condition. */ + return count; + } + + if (*stp < param->min_step) { + /* The step is the minimum value. */ + return LBFGSERR_MINIMUMSTEP; + } + if (*stp > param->max_step) { + /* The step is the maximum value. */ + return LBFGSERR_MAXIMUMSTEP; + } + if (param->max_linesearch <= count) { + /* Maximum number of iteration. */ + return LBFGSERR_MAXIMUMLINESEARCH; + } + + (*stp) *= width; + } +} + + + +static int line_search_morethuente( + int n, + lbfgsfloatval_t *x, + lbfgsfloatval_t *f, + lbfgsfloatval_t *g, + lbfgsfloatval_t *s, + lbfgsfloatval_t *stp, + const lbfgsfloatval_t* xp, + const lbfgsfloatval_t* gp, + lbfgsfloatval_t *wa, + callback_data_t *cd, + const lbfgs_parameter_t *param + ) +{ + int count = 0; + int brackt, stage1, uinfo = 0; + lbfgsfloatval_t dg; + lbfgsfloatval_t stx, fx, dgx; + lbfgsfloatval_t sty, fy, dgy; + lbfgsfloatval_t fxm, dgxm, fym, dgym, fm, dgm; + lbfgsfloatval_t finit, ftest1, dginit, dgtest; + lbfgsfloatval_t width, prev_width; + lbfgsfloatval_t stmin, stmax; + + /* Check the input parameters for errors. */ + if (*stp <= 0.) { + return LBFGSERR_INVALIDPARAMETERS; + } + + /* Compute the initial gradient in the search direction. */ + vecdot(&dginit, g, s, n); + + /* Make sure that s points to a descent direction. */ + if (0 < dginit) { + return LBFGSERR_INCREASEGRADIENT; + } + + /* Initialize local variables. */ + brackt = 0; + stage1 = 1; + finit = *f; + dgtest = param->ftol * dginit; + width = param->max_step - param->min_step; + prev_width = 2.0 * width; + + /* + The variables stx, fx, dgx contain the values of the step, + function, and directional derivative at the best step. + The variables sty, fy, dgy contain the value of the step, + function, and derivative at the other endpoint of + the interval of uncertainty. + The variables stp, f, dg contain the values of the step, + function, and derivative at the current step. + */ + stx = sty = 0.; + fx = fy = finit; + dgx = dgy = dginit; + + for (;;) { + /* + Set the minimum and maximum steps to correspond to the + present interval of uncertainty. + */ + if (brackt) { + stmin = min2(stx, sty); + stmax = max2(stx, sty); + } else { + stmin = stx; + stmax = *stp + 4.0 * (*stp - stx); + } + + /* Clip the step in the range of [stpmin, stpmax]. */ + if (*stp < param->min_step) *stp = param->min_step; + if (param->max_step < *stp) *stp = param->max_step; + + /* + If an unusual termination is to occur then let + stp be the lowest point obtained so far. + */ + if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->xtol * stmax))) { + *stp = stx; + } + + /* + Compute the current value of x: + x <- x + (*stp) * s. + */ + veccpy(x, xp, n); + vecadd(x, s, *stp, n); + + /* Evaluate the function and gradient values. */ + *f = cd->proc_evaluate(cd->instance, x, g, cd->n, *stp); + vecdot(&dg, g, s, n); + + ftest1 = finit + *stp * dgtest; + ++count; + + /* Test for errors and convergence. */ + if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0)) { + /* Rounding errors prevent further progress. */ + return LBFGSERR_ROUNDING_ERROR; + } + if (*stp == param->max_step && *f <= ftest1 && dg <= dgtest) { + /* The step is the maximum value. */ + return LBFGSERR_MAXIMUMSTEP; + } + if (*stp == param->min_step && (ftest1 < *f || dgtest <= dg)) { + /* The step is the minimum value. */ + return LBFGSERR_MINIMUMSTEP; + } + if (brackt && (stmax - stmin) <= param->xtol * stmax) { + /* Relative width of the interval of uncertainty is at most xtol. */ + return LBFGSERR_WIDTHTOOSMALL; + } + if (param->max_linesearch <= count) { + /* Maximum number of iteration. */ + return LBFGSERR_MAXIMUMLINESEARCH; + } + if (*f <= ftest1 && fabs(dg) <= param->gtol * (-dginit)) { + /* The sufficient decrease condition and the directional derivative condition hold. */ + return count; + } + + /* + In the first stage we seek a step for which the modified + function has a nonpositive value and nonnegative derivative. + */ + if (stage1 && *f <= ftest1 && min2(param->ftol, param->gtol) * dginit <= dg) { + stage1 = 0; + } + + /* + A modified function is used to predict the step only if + we have not obtained a step for which the modified + function has a nonpositive function value and nonnegative + derivative, and if a lower function value has been + obtained but the decrease is not sufficient. + */ + if (stage1 && ftest1 < *f && *f <= fx) { + /* Define the modified function and derivative values. */ + fm = *f - *stp * dgtest; + fxm = fx - stx * dgtest; + fym = fy - sty * dgtest; + dgm = dg - dgtest; + dgxm = dgx - dgtest; + dgym = dgy - dgtest; + + /* + Call update_trial_interval() to update the interval of + uncertainty and to compute the new step. + */ + uinfo = update_trial_interval( + &stx, &fxm, &dgxm, + &sty, &fym, &dgym, + stp, &fm, &dgm, + stmin, stmax, &brackt + ); + + /* Reset the function and gradient values for f. */ + fx = fxm + stx * dgtest; + fy = fym + sty * dgtest; + dgx = dgxm + dgtest; + dgy = dgym + dgtest; + } else { + /* + Call update_trial_interval() to update the interval of + uncertainty and to compute the new step. + */ + uinfo = update_trial_interval( + &stx, &fx, &dgx, + &sty, &fy, &dgy, + stp, f, &dg, + stmin, stmax, &brackt + ); + } + + /* + Force a sufficient decrease in the interval of uncertainty. + */ + if (brackt) { + if (0.66 * prev_width <= fabs(sty - stx)) { + *stp = stx + 0.5 * (sty - stx); + } + prev_width = width; + width = fabs(sty - stx); + } + } + + return LBFGSERR_LOGICERROR; +} + + + +/** + * Define the local variables for computing minimizers. + */ +#define USES_MINIMIZER \ + lbfgsfloatval_t a, d, gamma, theta, p, q, r, s; + +/** + * Find a minimizer of an interpolated cubic function. + * @param cm The minimizer of the interpolated cubic. + * @param u The value of one point, u. + * @param fu The value of f(u). + * @param du The value of f'(u). + * @param v The value of another point, v. + * @param fv The value of f(v). + * @param du The value of f'(v). + */ +#define CUBIC_MINIMIZER(cm, u, fu, du, v, fv, dv) \ + d = (v) - (u); \ + theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \ + p = fabs(theta); \ + q = fabs(du); \ + r = fabs(dv); \ + s = max3(p, q, r); \ + /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \ + a = theta / s; \ + gamma = s * sqrt(a * a - ((du) / s) * ((dv) / s)); \ + if ((v) < (u)) gamma = -gamma; \ + p = gamma - (du) + theta; \ + q = gamma - (du) + gamma + (dv); \ + r = p / q; \ + (cm) = (u) + r * d; + +/** + * Find a minimizer of an interpolated cubic function. + * @param cm The minimizer of the interpolated cubic. + * @param u The value of one point, u. + * @param fu The value of f(u). + * @param du The value of f'(u). + * @param v The value of another point, v. + * @param fv The value of f(v). + * @param du The value of f'(v). + * @param xmin The maximum value. + * @param xmin The minimum value. + */ +#define CUBIC_MINIMIZER2(cm, u, fu, du, v, fv, dv, xmin, xmax) \ + d = (v) - (u); \ + theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \ + p = fabs(theta); \ + q = fabs(du); \ + r = fabs(dv); \ + s = max3(p, q, r); \ + /* gamma = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \ + a = theta / s; \ + gamma = s * sqrt(max2(0, a * a - ((du) / s) * ((dv) / s))); \ + if ((u) < (v)) gamma = -gamma; \ + p = gamma - (dv) + theta; \ + q = gamma - (dv) + gamma + (du); \ + r = p / q; \ + if (r < 0. && gamma != 0.) { \ + (cm) = (v) - r * d; \ + } else if (a < 0) { \ + (cm) = (xmax); \ + } else { \ + (cm) = (xmin); \ + } + +/** + * Find a minimizer of an interpolated quadratic function. + * @param qm The minimizer of the interpolated quadratic. + * @param u The value of one point, u. + * @param fu The value of f(u). + * @param du The value of f'(u). + * @param v The value of another point, v. + * @param fv The value of f(v). + */ +#define QUARD_MINIMIZER(qm, u, fu, du, v, fv) \ + a = (v) - (u); \ + (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a; + +/** + * Find a minimizer of an interpolated quadratic function. + * @param qm The minimizer of the interpolated quadratic. + * @param u The value of one point, u. + * @param du The value of f'(u). + * @param v The value of another point, v. + * @param dv The value of f'(v). + */ +#define QUARD_MINIMIZER2(qm, u, du, v, dv) \ + a = (u) - (v); \ + (qm) = (v) + (dv) / ((dv) - (du)) * a; + +/** + * Update a safeguarded trial value and interval for line search. + * + * The parameter x represents the step with the least function value. + * The parameter t represents the current step. This function assumes + * that the derivative at the point of x in the direction of the step. + * If the bracket is set to true, the minimizer has been bracketed in + * an interval of uncertainty with endpoints between x and y. + * + * @param x The pointer to the value of one endpoint. + * @param fx The pointer to the value of f(x). + * @param dx The pointer to the value of f'(x). + * @param y The pointer to the value of another endpoint. + * @param fy The pointer to the value of f(y). + * @param dy The pointer to the value of f'(y). + * @param t The pointer to the value of the trial value, t. + * @param ft The pointer to the value of f(t). + * @param dt The pointer to the value of f'(t). + * @param tmin The minimum value for the trial value, t. + * @param tmax The maximum value for the trial value, t. + * @param brackt The pointer to the predicate if the trial value is + * bracketed. + * @retval int Status value. Zero indicates a normal termination. + * + * @see + * Jorge J. More and David J. Thuente. Line search algorithm with + * guaranteed sufficient decrease. ACM Transactions on Mathematical + * Software (TOMS), Vol 20, No 3, pp. 286-307, 1994. + */ +static int update_trial_interval( + lbfgsfloatval_t *x, + lbfgsfloatval_t *fx, + lbfgsfloatval_t *dx, + lbfgsfloatval_t *y, + lbfgsfloatval_t *fy, + lbfgsfloatval_t *dy, + lbfgsfloatval_t *t, + lbfgsfloatval_t *ft, + lbfgsfloatval_t *dt, + const lbfgsfloatval_t tmin, + const lbfgsfloatval_t tmax, + int *brackt + ) +{ + int bound; + int dsign = fsigndiff(dt, dx); + lbfgsfloatval_t mc; /* minimizer of an interpolated cubic. */ + lbfgsfloatval_t mq; /* minimizer of an interpolated quadratic. */ + lbfgsfloatval_t newt; /* new trial value. */ + USES_MINIMIZER; /* for CUBIC_MINIMIZER and QUARD_MINIMIZER. */ + + /* Check the input parameters for errors. */ + if (*brackt) { + if (*t <= min2(*x, *y) || max2(*x, *y) <= *t) { + /* The trival value t is out of the interval. */ + return LBFGSERR_OUTOFINTERVAL; + } + if (0. <= *dx * (*t - *x)) { + /* The function must decrease from x. */ + return LBFGSERR_INCREASEGRADIENT; + } + if (tmax < tmin) { + /* Incorrect tmin and tmax specified. */ + return LBFGSERR_INCORRECT_TMINMAX; + } + } + + /* + Trial value selection. + */ + if (*fx < *ft) { + /* + Case 1: a higher function value. + The minimum is brackt. If the cubic minimizer is closer + to x than the quadratic one, the cubic one is taken, else + the average of the minimizers is taken. + */ + *brackt = 1; + bound = 1; + CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt); + QUARD_MINIMIZER(mq, *x, *fx, *dx, *t, *ft); + if (fabs(mc - *x) < fabs(mq - *x)) { + newt = mc; + } else { + newt = mc + 0.5 * (mq - mc); + } + } else if (dsign) { + /* + Case 2: a lower function value and derivatives of + opposite sign. The minimum is brackt. If the cubic + minimizer is closer to x than the quadratic (secant) one, + the cubic one is taken, else the quadratic one is taken. + */ + *brackt = 1; + bound = 0; + CUBIC_MINIMIZER(mc, *x, *fx, *dx, *t, *ft, *dt); + QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt); + if (fabs(mc - *t) > fabs(mq - *t)) { + newt = mc; + } else { + newt = mq; + } + } else if (fabs(*dt) < fabs(*dx)) { + /* + Case 3: a lower function value, derivatives of the + same sign, and the magnitude of the derivative decreases. + The cubic minimizer is only used if the cubic tends to + infinity in the direction of the minimizer or if the minimum + of the cubic is beyond t. Otherwise the cubic minimizer is + defined to be either tmin or tmax. The quadratic (secant) + minimizer is also computed and if the minimum is brackt + then the the minimizer closest to x is taken, else the one + farthest away is taken. + */ + bound = 1; + CUBIC_MINIMIZER2(mc, *x, *fx, *dx, *t, *ft, *dt, tmin, tmax); + QUARD_MINIMIZER2(mq, *x, *dx, *t, *dt); + if (*brackt) { + if (fabs(*t - mc) < fabs(*t - mq)) { + newt = mc; + } else { + newt = mq; + } + } else { + if (fabs(*t - mc) > fabs(*t - mq)) { + newt = mc; + } else { + newt = mq; + } + } + } else { + /* + Case 4: a lower function value, derivatives of the + same sign, and the magnitude of the derivative does + not decrease. If the minimum is not brackt, the step + is either tmin or tmax, else the cubic minimizer is taken. + */ + bound = 0; + if (*brackt) { + CUBIC_MINIMIZER(newt, *t, *ft, *dt, *y, *fy, *dy); + } else if (*x < *t) { + newt = tmax; + } else { + newt = tmin; + } + } + + /* + Update the interval of uncertainty. This update does not + depend on the new step or the case analysis above. + + - Case a: if f(x) < f(t), + x <- x, y <- t. + - Case b: if f(t) <= f(x) && f'(t)*f'(x) > 0, + x <- t, y <- y. + - Case c: if f(t) <= f(x) && f'(t)*f'(x) < 0, + x <- t, y <- x. + */ + if (*fx < *ft) { + /* Case a */ + *y = *t; + *fy = *ft; + *dy = *dt; + } else { + /* Case c */ + if (dsign) { + *y = *x; + *fy = *fx; + *dy = *dx; + } + /* Cases b and c */ + *x = *t; + *fx = *ft; + *dx = *dt; + } + + /* Clip the new trial value in [tmin, tmax]. */ + if (tmax < newt) newt = tmax; + if (newt < tmin) newt = tmin; + + /* + Redefine the new trial value if it is close to the upper bound + of the interval. + */ + if (*brackt && bound) { + mq = *x + 0.66 * (*y - *x); + if (*x < *y) { + if (mq < newt) newt = mq; + } else { + if (newt < mq) newt = mq; + } + } + + /* Return the new trial value. */ + *t = newt; + return 0; +} + + + + + +static lbfgsfloatval_t owlqn_x1norm( + const lbfgsfloatval_t* x, + const int start, + const int n + ) +{ + int i; + lbfgsfloatval_t norm = 0.; + + for (i = start;i < n;++i) { + norm += fabs(x[i]); + } + + return norm; +} + +static void owlqn_pseudo_gradient( + lbfgsfloatval_t* pg, + const lbfgsfloatval_t* x, + const lbfgsfloatval_t* g, + const int n, + const lbfgsfloatval_t c, + const int start, + const int end + ) +{ + int i; + + /* Compute the negative of gradients. */ + for (i = 0;i < start;++i) { + pg[i] = g[i]; + } + + /* Compute the psuedo-gradients. */ + for (i = start;i < end;++i) { + if (x[i] < 0.) { + /* Differentiable. */ + pg[i] = g[i] - c; + } else if (0. < x[i]) { + /* Differentiable. */ + pg[i] = g[i] + c; + } else { + if (g[i] < -c) { + /* Take the right partial derivative. */ + pg[i] = g[i] + c; + } else if (c < g[i]) { + /* Take the left partial derivative. */ + pg[i] = g[i] - c; + } else { + pg[i] = 0.; + } + } + } + + for (i = end;i < n;++i) { + pg[i] = g[i]; + } +} + +static void owlqn_project( + lbfgsfloatval_t* d, + const lbfgsfloatval_t* sign, + const int start, + const int end + ) +{ + int i; + + for (i = start;i < end;++i) { + if (d[i] * sign[i] <= 0) { + d[i] = 0; + } + } +} diff --git a/src/AI/Calculation.hs b/src/AI/Calculation.hs new file mode 100644 --- /dev/null +++ b/src/AI/Calculation.hs @@ -0,0 +1,25 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides common calculation functions +-- +-- +---------------------------------------------------- + + +module AI.Calculation ( + module AI.Calculation.Activation, + module AI.Calculation.Cost, + module AI.Calculation.Gradients, + module AI.Calculation.NetworkOutput + ) where + +import AI.Calculation.Activation +import AI.Calculation.Cost +import AI.Calculation.Gradients +import AI.Calculation.NetworkOutput diff --git a/src/AI/Calculation/Activation.hs b/src/AI/Calculation/Activation.hs new file mode 100644 --- /dev/null +++ b/src/AI/Calculation/Activation.hs @@ -0,0 +1,66 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides common activation functions +-- and their derivative +-- +-- +---------------------------------------------------- + + +module AI.Calculation.Activation ( + Activation(..), + getActivation, + getDerivative + ) where + +import AI.Signatures + + +-- | Represents the activation of +-- each neuron in the neural network +data Activation = Sigmoid -- ^ The sigmoid activation function + | HyperbolicTangent -- ^ The hyperbolic tangent activation function + + +-- | Get the activation function associated with an activation +getActivation :: Activation -> ActivationFunction +getActivation Sigmoid = sigmoid +getActivation HyperbolicTangent = hTangent + + +-- | Get the derivative function associated with an activation +getDerivative :: Activation -> DerivativeFunction +getDerivative Sigmoid = sigmoidDeriv +getDerivative HyperbolicTangent = hTangentDeriv + + +-- The sigmoid function +sigmoid :: ActivationFunction +sigmoid x = (1 / (1 + exp(-x))) + + +-- The derivative of the sigmoid function +-- +-- NOTE: The derivative is (sigmoid x) * (1 - sigmoid x) +-- NOT (x * (1 - x)) +sigmoidDeriv :: DerivativeFunction +sigmoidDeriv x = (sigmoid x) * (1 - (sigmoid x)) + + +-- The hyperbolic tangent function +hTangent :: ActivationFunction +hTangent x = tanh x + + +-- The derivative of the hyperbolic tangent +-- +-- NOTE: The derivative is 1 - (tanh x)^2 +-- NOT 1 - x^2 +hTangentDeriv :: DerivativeFunction +hTangentDeriv x = 1 - ((tanh x) ** 2) diff --git a/src/AI/Calculation/Cost.hs b/src/AI/Calculation/Cost.hs new file mode 100644 --- /dev/null +++ b/src/AI/Calculation/Cost.hs @@ -0,0 +1,155 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides common cost functions +-- and their derivatives +-- +-- +---------------------------------------------------- + + +module AI.Calculation.Cost ( + Cost(..), + getCostFunction, + getCostDerivative, +) where + +import Data.Packed.Matrix +import Data.Packed.Vector +import Numeric.Container + +import AI.Signatures +import AI.Network + + +-- | Represents the cost model +-- of the Neural Network +data Cost = MeanSquared -- ^ The mean-squared cost + | Logistic -- ^ The logistic cost + + +-- | Gets the cost function associated +-- with the cost model +getCostFunction :: Cost -> CostFunction +getCostFunction = generalCost . getErrorFunction + + +-- | Gets the cost derivative associated +-- with the cost model +getCostDerivative :: Cost -> CostDerivative +getCostDerivative MeanSquared = meanSquaredDerivative +getCostDerivative Logistic = logisticDerivative + + +-- The general cost function that can be extended by +-- partial function application +generalCost :: ErrorFunction -- The error function to be used + -> Network -- The neural network of interest + -> Matrix Double -- The matrix of the inputs, where the ith row + -- is the input vector of a training set + -> Matrix Double -- The matrix of the expected output, where the ith + -- row is the expected output vector of a + -- training set + -> Double -- Returns the cost by comparing the network's + -- output neurons and the expected output matrix +generalCost errorF nn inMatrix exMatrix = + let activF = toActivation nn -- The activation function + ws = toWeightMatrices nn -- The list of weight matrices + la = toLambda nn -- The regularization constant + n = rows exMatrix -- Get us the number of training sets + + -- Set up the bias neurons for each forward propagation + fBias = \m -> (fromColumns . ((fromList (replicate n 1)):) . toColumns) m + + -- This is forward propagation right here + -- We propagate forward by using the functional + -- foldl accumulating over the weight matrices + f = \a w -> fBias (mapMatrix activF (a `multiply` w)) + + -- Prepare the inMatrix to be used + inMatrix' = fBias inMatrix + + -- We fold f over the weight matrices accumulating + -- the activation matrix. The resulting activation matrix + -- is our output. + -- We take out the bias neurons in our output layer + outMatrix = (fromColumns . tail . toColumns) $ foldl f inMatrix' ws + + -- We first convert the outMatrix and the exMatrix into + -- list of rows. Then we zip them into oVec and eVec + -- Then get the error between the oVec and the eVec using our + -- errorF (error function) and the result + -- is a list of errors. Then we create a vector + -- from the list, thus resulting in an errorVec + errorVec = fromList [errorF oVec eVec | + (oVec, eVec) <- zip (toRows outMatrix) (toRows exMatrix)] + + -- Now we get the un-regularized cost by + -- summing the elements and taking the average + -- by dividing by the number of the training sets + j = (1 / fromIntegral n) * sumElements errorVec + + -- Now we get the vector representation of the + -- weights to prepare for regularization + wsFlattened = toWeights nn + in + -- Finally, we regularize our cost by using the lambda constant + j + (la / (2.0 * fromIntegral n)) * (sumElements $ mapVector (**2) wsFlattened) + + +-- This is the general error function +-- It requires a function that will calculate an +-- error when given a calculated value and an +-- expected value +generalErrorFunction :: (Double -> Double -> Double) -- The function to calculate an "error" + -- Given a calculated value and an + -- expected value + -> ErrorFunction -- Returns the error function +generalErrorFunction errF calVec exVec = + let n = dim calVec -- Get us the size of the vectors + + -- Get us the error vector between the calculated + -- vector and the expected vector by zipping + -- the error function + errorVec = zipVectorWith errF calVec exVec + in + -- Now we take the average of the sum of the errors + (1 / fromIntegral n) * (sumElements errorVec) + + +-- Returns the error function given a cost detail +getErrorFunction :: Cost -> ErrorFunction +getErrorFunction MeanSquared = generalErrorFunction meanSquaredError +getErrorFunction Logistic = generalErrorFunction logisticError + + +-- The mean squared error function +meanSquaredError :: Double -> Double -> Double +meanSquaredError cal ex = + (cal - ex) ** 2 + + +-- The derivative of the mean squared error function +-- with respect to each parameter +meanSquaredDerivative :: CostDerivative +meanSquaredDerivative (Network {derivative = df}) inMat actMat exMat = + mapMatrix (*2) ((mapMatrix df inMat) `mul` (actMat `sub` exMat)) + + +-- The logistic error function +-- Also known as the cross-entropy +-- error function +logisticError :: Double -> Double -> Double +logisticError cal ex= + (-1) * ((ex) * (log cal) + (1 - ex) * (log (1 - cal))) + + +-- The derivative of the logistic error function +-- with respect to each parameter +logisticDerivative :: CostDerivative +logisticDerivative _ _ actMat exMat = actMat `sub` exMat diff --git a/src/AI/Calculation/Gradients.hs b/src/AI/Calculation/Gradients.hs new file mode 100644 --- /dev/null +++ b/src/AI/Calculation/Gradients.hs @@ -0,0 +1,205 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module represents ways to calculate the gradient +-- vector of the weights of the neural network +-- +-- Backpropagation should always be preferred over +-- the Numerical Gradient method +-- +-- +---------------------------------------------------- + + +module AI.Calculation.Gradients ( + backpropagation, + numericalGradients + ) where + +import Numeric.Container + +import AI.Signatures +import AI.Network + + +-- | Calculate the analytical gradient of the weights of the network +-- by using backpropagation +backpropagation :: GradientFunction +backpropagation _ outputDeltasF nn inMatrix exMatrix = + let n = rows exMatrix -- Get us the number of training sets + + -- Get important informations from the neural network + activF = toActivation nn + df = toDerivative nn + ws = toWeightMatrices nn + la = toLambda nn + + -- Function to set up the bias neurons for each forward propagation + fBias = fromColumns . ((fromList (replicate n 1)):) . toColumns + + -- This is forward propagation right here + -- We propagate forward by using the functional + -- foldl accumulating over the weight matrices + -- NOTE: Also, for the ability to use more complex activation + -- function with non-trivial derivatives, we also accumulate + -- the weighted inputs to each layer, so the accumulation + -- for each forward propagation is a tuple of the activation + -- of each neuron in the layer and the weighted inputs into the + -- layer. + f = \p w -> ((fBias . mapMatrix activF) (p `multiply` w), fBias $ p `multiply` w) + + -- Prepare the inMatrix to be used + inMatrix' = fBias inMatrix + + -- Forward propagate each row of the inputs matrix and matrix-multiply it with the + -- weight matrices calculated before. + -- NOTE: Forward propagation can be calculated efficiently by using + -- foldl where the initial value is the input matrix and we calculate + -- and simultaneously accumulate the activation values of each layer + activations = foldl (\(a@(p,_):as) w -> ((f p w):a:as)) [(inMatrix',inMatrix')] ws + + -- Helper function to remove bias neurons + fRemoveBias = fromColumns . tail . toColumns + + -- Because we cannot possibly calculate the outputs node deltas, + -- the user must supply a function that will do that + -- We pass in the weighted inputs to the output neurons, + -- the activation values of the output neurons, + -- the expected matrix of the training set, + -- the derivative of the activation function of the networks + -- And we expect it to return for us the output nodes + -- deltas for us to propagate backwards to each layer + initialDeltas = outputDeltasF nn ((fRemoveBias . snd . head) activations) + ((fRemoveBias . fst . head) activations) exMatrix + + -- This one line is basically the entire backpropagation + -- + -- NOTE: Backpropagation can be computed efficiently + -- using foldl. Because the activations we calculated above + -- are in reverse order, we can efficiently backpropagate + -- the initial output nodes deltas by simply folding + -- backwards on the reverse of the weights + -- The initial value for foldl is our initialDeltas + -- and we accumulate the node deltas of each + -- previous layer. + -- + -- NOTE: Because we also accumulated the weighted inputs + -- to each layer, we can use more exotic activation + -- function instead of the ones with trivial derivatives. + -- + -- Example: Instead of using the derivative of the + -- sigmoid function as x * (1 - x), where x is the + -- "sigmoided value", we can actually use the real + -- derivative, which is (sigmoid x) * (1 - sigmoid x) + -- where x is the weighted input + -- + -- This allows us to use more exotic activation + -- function whose derivatives is non-trivial + allDeltas = foldl (\(d:ds) (as,w) -> (fRemoveBias $ (d `multiply` (trans w)) `mul` (mapMatrix df as)):d:ds) + [initialDeltas] (zip ((tail . map snd) activations) (reverse ws)) + + -- Now this is where we finally calculate the gradients + -- by multiplying activations of each layer to the + -- node deltas of the next layer + grads = [[a `outer` d | (a, d) <- zip (toRows as) (toRows deltas)] + | (as, deltas) <- zip ((tail . map fst) activations) (reverse allDeltas)] + + -- zeroF gets us a zero matrix given a row and a column + -- I believe the HMatrix package must have a function to + -- create zero matrices, but I haven't fond it yet... >_< + zeroF = \m -> buildMatrix (rows m) (cols m) (\_ -> 0.0) + + -- Now we add all of the gradients together + -- and the average of the gradients by dividing by + -- the number of the training sets + gradsSums = map (mapMatrix (/(fromIntegral n))) [foldl add ((zeroF . head) g) g | g <- grads] + in + -- And after all that exhaustive work, we flatten the matrices into + -- one big vector and add regularization to it + zipVectorWith (\x y -> x + (la / fromIntegral n) * y) ((join . map flatten) (reverse gradsSums)) (toWeights nn) + + +-- | NOTE: This should only be used as a last resort +-- if for some reason (bugs?) the backpropagation +-- algorithm does not give you good gradients +-- +-- The numerical algorithm requires two forward +-- propagations, while the backpropagation algorithm +-- only requires one, so this is more costly +-- +-- Also, analytical gradients almost always perform +-- better than numerical gradients +-- +-- User must provide an epsilon value. +-- Make sure to use a very small value for the epsilon +-- for more accurate gradients +numericalGradients :: Double -- ^ The epsilon + -> GradientFunction -- ^ Returns a gradient function + -- that calculates the numerical + -- gradients of the weights +numericalGradients epsilon costF _ nn inMat exMat = + let plusE = \x -> x + epsilon -- Add epsilon to the argument + minusE = \x -> x - epsilon -- Subtract epsilon from the argument + + -- Get the vector representation of the weights + params = toWeights nn + + -- Calculate a matrix of parameters that have been + -- modified by adding and subtracting the epsilon value + -- The result is two lists whose element is a vector + -- of the parameters that have been modified + dx1s = (toRows . mapElementToVector (modifyElementAt plusE)) params + dx2s = (toRows . mapElementToVector (modifyElementAt minusE)) params + + f = \ws -> costF (setWeights nn ws) inMat exMat + -- Now we calculate the costs of each modified parameters + cost1 = (fromList . map f) dx1s + cost2 = (fromList . map f) dx2s + in + -- Use the (f(x+e) - (f(x-e)))/(2*e) to get the gradients + mapVector (/ (2 * epsilon)) $ cost1 `sub` cost2 + + +-- Map over every single element of a vector +-- and apply the function on each vector +-- This returns a matrix where the ith row is +-- a vector whose ith element is applied +-- to the function f +mapElementToVector :: (Vector Double -> Int -> Vector Double) + -> Vector Double + -> Matrix Double +mapElementToVector f vec = + let n = dim vec + -- Get the size of the vector + + -- Apply f over each element indexed by i + -- and we join the list of vectors into a giant + -- vector + flattened = join [f vec i | i <- [0..n - 1]] + in + -- Reshape the flattened vector into a matrix + -- by using the size of the vector calculated before + reshape n flattened + + +-- Modify one single element of the vector +-- by applying f to it +modifyElementAt :: (Double -> Double) -- The function to be applied + -> Vector Double -- The vector of interest + -> Int -- The index of the element ot be modified + -> Vector Double -- The resulting modified vector +modifyElementAt f vec index = + -- If the index passed into us is the index, we apply the + -- function to the element, otherwise we leave it alone + let g = \i v -> if' (i == index) (f v) (v) in + mapVectorWithIndex g vec + + +if' :: Bool -> a -> a -> a +if' True x _ = x +if' False _ y = y diff --git a/src/AI/Calculation/NetworkOutput.hs b/src/AI/Calculation/NetworkOutput.hs new file mode 100644 --- /dev/null +++ b/src/AI/Calculation/NetworkOutput.hs @@ -0,0 +1,35 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides forward propagation +-- to let the user get the output of the neural +-- network given an input vector +-- +-- +---------------------------------------------------- + + +module AI.Calculation.NetworkOutput ( + networkOutput + ) where + +import Numeric.Container + +import AI.Network + + +-- | Forward propagate to get the network's output +networkOutput :: Network -- ^ The neural network of interest + -> Vector Double -- ^ The input vector + -> Vector Double -- ^ The output vector of the output neurons +networkOutput nn inVec = + let activF = toActivation nn + ws = toWeightMatrices nn + fBias = \v -> fromList $ 1.0:(toList v :: [Double]) + in + foldl (\as w -> mapVector activF ((fBias as) `vXm` w)) inVec ws diff --git a/src/AI/Model.hs b/src/AI/Model.hs new file mode 100644 --- /dev/null +++ b/src/AI/Model.hs @@ -0,0 +1,24 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- Provides models interface for easy initialization and +-- training of neural networks +-- +-- +---------------------------------------------------- + + +module AI.Model ( + module AI.Model.GenericModel, + module AI.Model.General, + module AI.Model.Classification + ) where + +import AI.Model.GenericModel +import AI.Model.General +import AI.Model.Classification diff --git a/src/AI/Model/Classification.hs b/src/AI/Model/Classification.hs new file mode 100644 --- /dev/null +++ b/src/AI/Model/Classification.hs @@ -0,0 +1,41 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides an initialization for +-- a classification neural network model +-- +-- NOTE: This theoretically should be faster than +-- the General model if used for classification +-- +-- +---------------------------------------------------- + + +module AI.Model.Classification ( + initializeClassification + ) where + +import Data.Packed.Vector +import System.Random + +import AI.Calculation +import AI.Network +import AI.Model.GenericModel + + +-- | Make a neural network model +-- that should be used for classification +-- using the Sigmoid as the activation model +-- and Logistic as the cost model +initializeClassification :: [Int] -- ^ The architecture of the neural network + -> Double -- ^ The regularization constant + -> StdGen -- ^ The random generator + -> GenericModel -- ^ Returns the initialized classification + -- model +initializeClassification arch la gen = + initializeModel Sigmoid Logistic arch la gen diff --git a/src/AI/Model/General.hs b/src/AI/Model/General.hs new file mode 100644 --- /dev/null +++ b/src/AI/Model/General.hs @@ -0,0 +1,42 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides an initialization for +-- a general neural network model that can do +-- either regression or classification +-- +-- If for regression, the training +-- data must be normalized by user to have +-- range of [-1,1] +-- +---------------------------------------------------- + + +module AI.Model.General ( + initializeGeneral + ) where + +import Data.Packed.Vector +import System.Random + +import AI.Calculation +import AI.Network +import AI.Model.GenericModel + + +-- | This is a general neural network +-- model that can be used for classification +-- or regression using HyperbolicTangent +-- as the activation model and MeanSquared as +-- the cost model +initializeGeneral :: [Int] -- ^ The architecture of the neural network + -> Double -- ^ The regularization constant + -> StdGen -- ^ The random generator + -> GenericModel -- ^ Returns the initialized general model +initializeGeneral arch la gen = + initializeModel HyperbolicTangent MeanSquared arch la gen diff --git a/src/AI/Model/GenericModel.hs b/src/AI/Model/GenericModel.hs new file mode 100644 --- /dev/null +++ b/src/AI/Model/GenericModel.hs @@ -0,0 +1,88 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides a generic module for +-- initiialization and training of neural networks +-- +-- User must provide the needed functions +-- +-- +---------------------------------------------------- + + +module AI.Model.GenericModel ( + GenericModel(..), + initializeModel, + getOutput, + trainModel + ) where + +import System.Random +import Data.Packed.Vector +import Data.Packed.Matrix + +import AI.Signatures +import AI.Calculation +import AI.Network +import AI.Training + + +-- | Generic neural network model for expansion +data GenericModel = GenericModel + { + cost :: Cost, -- ^ The cost model of the model + net :: Network -- ^ The neural network to be used for modeling + } + + +-- | Initialize neural network model with the weights +-- randomized within [-1.0,1.0] +initializeModel :: Activation -- ^ The activation model of each neuron + -> Cost -- ^ The cost model of the output neurons + -- compared to the expected output + -> [Int] -- ^ The architecture of the network + -- e.g., a 2-3-1 architecture would be [2,3,1] + -> Double -- ^ The regularization constant + -- should be 0 if you do not want regularization + -> StdGen -- ^ The random generator + -> GenericModel -- ^ Returns the initialized model +initializeModel ac co arch la gen = + let n = foldl (+) 0 [((x + 1) * y) | (x,y) <- zip arch (tail arch)] + ws = (fromList . take n) (randomRs (-1.0, 1.0) gen :: [Double]) + in + GenericModel { cost = co, + net = Network { activation = getActivation ac, + derivative = getDerivative ac, + lambda = la, + weights = ws, + architecture = arch + } + } + + +-- | Get the output of the model +getOutput :: GenericModel -- ^ The model of interest + -> Vector Double -- ^ The input vector to the input layer + -> Vector Double -- ^ The output of the network model +getOutput (GenericModel {net = nn}) input = networkOutput nn input + + +-- | Train the model given the parameters and the training algorithm +trainModel :: GenericModel -- ^ The model to be trained + -> TrainingAlgorithm -- ^ The training algorithm to be used + -> Double -- ^ The precision to train with regards to + -- the cost function + -> Int -- ^ The maximum amount of epochs to train + -> Matrix Double -- ^ The input matrix + -> Matrix Double -- ^ The expected output matrix + -> GenericModel -- ^ Returns the trained model +trainModel (GenericModel {cost = co, net = nn}) algo prec epochs inMat exMat = + let trainedNet = trainNetwork algo co backpropagation nn prec epochs inMat exMat in + GenericModel { net = trainedNet, + cost = co + } diff --git a/src/AI/Network.hs b/src/AI/Network.hs new file mode 100644 --- /dev/null +++ b/src/AI/Network.hs @@ -0,0 +1,101 @@ +---------------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- Efficient representation of an Artificial Neural Network +-- using vector to represent the weights between each layer +-- +-- This module provides the neural network data representation +-- that will be used extensively +-- +-- +--------------------------------------------------------- + + +module AI.Network ( + Network(..), + toActivation, toDerivative, + toLambda, toWeights, + toWeightMatrices, toArchitecture, + setActivation, setDerivative, + setLambda, setWeights, + setArchitecture + ) where + +import Data.Packed.Vector +import Data.Packed.Matrix + + +-- | The representation of an Artificial Neural Network +data Network = Network + { + activation :: (Double -> Double), -- ^ The activation function for each + -- neuron + derivative :: (Double -> Double), -- ^ The derivative of the activation + -- function + lambda :: Double, -- ^ The regularization constant + weights :: Vector Double, -- ^ The vector of the weights between each + -- layer of the neural network + architecture :: [Int] -- ^ The architecture of the neural + -- networks. + -- + -- e.g., a network of an architecture + -- of 2-3-1 would have an architecture + -- representation of [2,3,1] + -- + -- NOTE: The library will automatically create + -- a bias neuron in each layer, so you do not need + -- to state them explicitly + } + + +-- Self-explanatory +toActivation :: Network -> (Double -> Double) +toActivation (Network {activation = f}) = f + +toDerivative :: Network -> (Double -> Double) +toDerivative (Network {derivative = df}) = df + +toLambda :: Network -> Double +toLambda (Network {lambda = la}) = la + +toWeights :: Network -> Vector Double +toWeights (Network {weights = w}) = w + +-- | Get the list of matrices of weights between +-- each layer. This can be more useful +-- than the barebone vector representation +-- of the weights +toWeightMatrices :: Network -> [Matrix Double] +toWeightMatrices (Network {weights = ws, architecture = arch}) = + let elems = 0:[((x + 1) * y) | (x,y) <- zip arch (tail arch)] in + [reshape i v | (i, v) <- zip (tail arch) (takesV (tail elems) ws)] + +toArchitecture :: Network -> [Int] +toArchitecture (Network {architecture = a}) = a + + +-- Self-explanatory +setActivation :: Network -> (Double -> Double) -> Network +setActivation (Network {derivative = df, lambda = la, weights = w, architecture = a}) f = + (Network {activation = f, derivative = df, lambda = la, weights = w, architecture = a}) + +setDerivative :: Network -> (Double -> Double) -> Network +setDerivative (Network {activation = f, lambda = la, weights = w, architecture = a}) df = + (Network {activation = f, derivative = df, lambda = la, weights = w, architecture = a}) + +setLambda :: Network -> Double -> Network +setLambda (Network {activation = f, derivative = df, weights = w, architecture = a}) la = + Network {activation = f, derivative = df, lambda = la, weights = w, architecture = a} + +setWeights :: Network -> Vector Double -> Network +setWeights (Network {activation = f, derivative = df, lambda = la, architecture = a}) w = + (Network {activation = f, derivative = df, lambda = la, weights = w, architecture = a}) + +setArchitecture :: Network -> [Int] -> Network +setArchitecture (Network {activation = f, derivative = df, lambda = la, weights = w}) a = + (Network {activation = f, derivative = df, lambda = la, weights = w, architecture = a}) diff --git a/src/AI/Signatures.hs b/src/AI/Signatures.hs new file mode 100644 --- /dev/null +++ b/src/AI/Signatures.hs @@ -0,0 +1,93 @@ +--------------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides the signatures for needed +-- functions in a neural network +-- +-- +--------------------------------------------------------- + + +module AI.Signatures ( + ActivationFunction, + DerivativeFunction, + ErrorFunction, + CostFunction, + CostDerivative, + GradientFunction + ) where + +import Data.Packed.Matrix +import Data.Packed.Vector + +import AI.Network + + +-- | Type that represents the activation function +type ActivationFunction = Double -> Double + + +-- | Type that represents the derivative of the activation function +-- +-- NOTE: The derivative can be non-trivial and must be continuous +type DerivativeFunction = Double -> Double + + +-- | Type that represents the error function +-- between the calculated output vector +-- and the expected output vector +type ErrorFunction = Vector Double -- ^ The calculated output vector + -> Vector Double -- ^ The expected output vector + -> Double -- ^ Returns the error of how far off + -- the calculated vector is from the + -- expected vector + + +-- | Type that represents the function +-- that can calculate the total cost of the neural networks +-- given the neural networks, the input matrix and an expected output matrix +type CostFunction = Network -- ^ The neural networks of interest + -> Matrix Double -- ^ The input matrix, where the ith row + -- is the input vector of a training set + -> Matrix Double -- ^ The expected output matrix, where the + -- ith row is the expected output vector + -- of a training set + -> Double -- ^ Returns the cost of the calculated output vector + -- from the neural network and the given + -- expected output vector + + +-- | Type that represents the cost function derivative. +-- on the output nodes +type CostDerivative = Network -- ^ The neural networks of interest + -> Matrix Double -- ^ The matrix of inputs where the ith row + -- is the ith training set + -> Matrix Double -- ^ The matrix of calculated outputs where the + -- ith row is the ith training set + -> Matrix Double -- ^ The matrix of expected outputs where the + -- ith row is the ith expected output of + -- of the training set + -> Matrix Double -- ^ Returns the matrix of the derivatives + -- of the cost of the output nodes + -- compared to the expected matrix + + +-- | The type to represent a function that +-- can calculate the gradient vector +-- of the weights of the neural network +-- +-- NOTE: Must be supplied a function to calculate the cost, the +-- cost derivative of the output neurons, the neural network +-- the input matrix, and the expected output matrix +type GradientFunction = CostFunction -- ^ The cost function + -> CostDerivative -- ^ The cost derivative + -> Network -- ^ The neural network + -> Matrix Double -- ^ The input matrix + -> Matrix Double -- ^ The expected output matrix + -> Vector Double -- ^ Returns the gradient vector + -- of the weights diff --git a/src/AI/Training.hs b/src/AI/Training.hs new file mode 100644 --- /dev/null +++ b/src/AI/Training.hs @@ -0,0 +1,122 @@ +---------------------------------------------------- +-- | +-- Module : AI.Network +-- License : GPL +-- +-- Maintainer : Kiet Lam +-- +-- +-- This module provides training algorithms to train +-- a neural network given training data. +-- +-- User should only use LBFGS though because +-- it uses custom bindings to the C-library liblbfgs +-- +-- GSL's multivariate minimization algorithms are known to be inefficient +-- +-- and LBFGS outperforms them on many (of my) tests +-- +-- +---------------------------------------------------- + + +module AI.Training ( + TrainingAlgorithm(..), + trainNetwork + ) where + +import Numeric.GSL.Minimization +import Data.Packed.Vector +import Data.Packed.Matrix + +import AI.Training.Internal +import AI.Signatures +import AI.Calculation +import AI.Network + + +-- | The types of training algorithm to use +-- +-- NOTE: These are all batch training algorithms +data TrainingAlgorithm = GradientDescent -- ^ hmatrix's binding to GSL + | ConjugateGradient -- ^ hmatrix's binding to GSL + | BFGS -- ^ hmatrix's binding to GSL + | LBFGS -- ^ home-made binding to liblbfgs + deriving (Show, Read, Enum) + + +-- This function is needed to work with HMatrix's +-- multivariate minimization algorithms +vectorWeightToCost :: CostFunction -- The cost function + -> Network -- The neural network + -> Matrix Double -- The input matrix + -> Matrix Double -- The output matrix + -> Vector Double -- The vector weights + -> Double -- Returns the calculated cost +vectorWeightToCost costF nn inMat exMat ws = costF (setWeights nn ws) inMat exMat + + +-- This function is needed to work with HMatrix's +-- multivariate minimization algorithms +vectorWeightToGradients :: GradientFunction -- The function can can calculate the + -- gradient vector given a cost model + -> Cost -- the cost model + -> Network -- The neural network + -> Matrix Double -- The input matrix + -> Matrix Double -- The output matrix + -> Vector Double -- The vector weights + -> Vector Double -- Returns the vector gradients +vectorWeightToGradients gradF cost nn inMat exMat ws = + gradF (getCostFunction cost) (getCostDerivative cost) (setWeights nn ws) inMat exMat + + +-- | Train the neural network given a training algorithm, +-- the training parameters and the training data +trainNetwork :: TrainingAlgorithm -- ^ The training algorithm to use + -> Cost -- ^ The cost model of the neural network + -> GradientFunction -- ^ The function that can calculate the + -- gradients vector + -> Network -- ^ The network to be trained + -> Double -- ^ The precision of the training with regards + -- to the cost function + -> Int -- ^ The maximum number of iterations + -> Matrix Double -- ^ The input matrix + -> Matrix Double -- ^ The expected output matrix + -> Network -- ^ Returns the trained network +trainNetwork algo cost gradF nn prec iterations inMat exMat = + let ws = toWeights nn -- Get the initial weights of the network + + -- f represents the cost function to minimize + f = vectorWeightToCost (getCostFunction cost) nn inMat exMat + + -- g represents the function that can calculate the gradient + -- vector of the parameters (the weights) + g = vectorWeightToGradients gradF cost nn inMat exMat + + -- Get the training algorithm + trainAlgo = getTrainAlgo algo + + -- Set the tol and initial step size + initStepSize = 0.1 + tol = 0.1 + + -- Use the training algorithm to train the weights + trainedWeights = trainAlgo prec iterations initStepSize tol f g ws + in + setWeights nn trainedWeights + + +-- Auxilary function for trainNetwork +getTrainAlgo :: TrainingAlgorithm + -> Double + -> Int + -> Double + -> Double + -> (Vector Double -> Double) + -> (Vector Double -> Vector Double) + -> Vector Double + -> Vector Double +getTrainAlgo GradientDescent prec iter step tol f df initVec = fst $ minimizeVD SteepestDescent prec iter step tol f df initVec +getTrainAlgo ConjugateGradient prec iter step tol f df initVec = fst $ minimizeVD ConjugatePR prec iter step tol f df initVec +getTrainAlgo BFGS prec iter step tol f df initVec = fst $ minimizeVD VectorBFGS2 prec iter step tol f df initVec +getTrainAlgo LBFGS prec iter step tol f df initVec = minimizeLBFGS prec iter step tol f df initVec diff --git a/src/AI/Training/Internal.hs b/src/AI/Training/Internal.hs new file mode 100644 --- /dev/null +++ b/src/AI/Training/Internal.hs @@ -0,0 +1,29 @@ +-- ------------------------------------------------- +-- +-- Author: Kiet Lam +-- File INTERNAL_HS +-- +-- ------------------------------------------------- +-- Last Updated: Time-stamp: <2012-01-18 18:37:35 (lam)> +-- +-- +-- +-- 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 . + + +module AI.Training.Internal ( + module AI.Training.Internal.LBFGSAux + ) where + +import AI.Training.Internal.LBFGSAux diff --git a/src/AI/Training/Internal/LBFGSAux.hs b/src/AI/Training/Internal/LBFGSAux.hs new file mode 100644 --- /dev/null +++ b/src/AI/Training/Internal/LBFGSAux.hs @@ -0,0 +1,123 @@ +-- ------------------------------------------------- +-- +-- Author: Kiet Lam +-- File LBFGSAUX_HS +-- +-- ------------------------------------------------- +-- Last Updated: Time-stamp: <2012-01-19 00:25:43 (lam)> +-- +-- +-- +-- 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 . + + +module AI.Training.Internal.LBFGSAux ( + minimizeLBFGS + ) where + + +import Data.Packed.Vector +import Foreign.C.Types +import Foreign.Ptr(Ptr, FunPtr) +import Foreign.Marshal.Array +import System.IO.Unsafe(unsafePerformIO) + + + +-- Don't make too much changes here + +type TV = CInt -> Ptr Double -> IO CInt +type TVV = CInt -> Ptr Double -> TV + + +aux_LToL :: ([Double] -> [Double]) -> TVV +aux_LToL f n1 p1 _ p2 = + do + v <- peekArray (fromIntegral n1) p1 + let vr = f v in + do + pokeArray p2 vr + return 0 + + +aux_LToD :: ([Double] -> Double) + -> CInt -> Ptr Double -> Double +aux_LToD f n p = + unsafePerformIO $ + do + v <- peekArray (fromIntegral n) p + return $ f v + + +foreign import ccall "wrapper" + mkListFun :: (CInt -> Ptr Double -> Double) + -> IO (FunPtr (CInt -> Ptr Double -> Double)) + + +foreign import ccall "wrapper" + mkListListFun :: (TVV) -> IO (FunPtr TVV) + + +foreign import ccall "lbfgs_aux.c minimizeLBFGS" + c_minimizeLBFGS :: Double + -> CInt + -> Double + -> Double + -> FunPtr (CInt -> Ptr Double -> Double) + -> FunPtr (CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt) + -> CInt -> Ptr Double + -> CInt -> Ptr Double + -> IO CInt + + +vecFuncToLFunc :: (Vector Double -> Vector Double) -> [Double] -> [Double] +vecFuncToLFunc f vec = (toList . f . fromList) vec + + +vecFuncToFunc :: (Vector Double -> Double) -> [Double] -> Double +vecFuncToFunc f vec = (f . fromList) vec + + +minimizeLBFGS_aux :: Double + -> Int + -> Double + -> Double + -> (Vector Double -> Double) + -> (Vector Double -> Vector Double) + -> Vector Double + -> [Double] +minimizeLBFGS_aux prec maxIter initStep tol f df initVec = + let f' = vecFuncToFunc f + df' = vecFuncToLFunc df + initVec' = toList initVec + n = length initVec' + in + unsafePerformIO $ withArray initVec' $ \ar -> allocaArray n $ \res -> + do + fp <- mkListFun (aux_LToD f') + dfp <- mkListListFun (aux_LToL df') + _ <- c_minimizeLBFGS prec (fromIntegral maxIter) initStep tol fp dfp (fromIntegral n) ar (fromIntegral n) res + peekArray n res + + +minimizeLBFGS :: Double + -> Int + -> Double + -> Double + -> (Vector Double -> Double) + -> (Vector Double -> Vector Double) + -> Vector Double + -> Vector Double +minimizeLBFGS prec maxIter initStep tol f df initVec = + fromList $ minimizeLBFGS_aux prec maxIter initStep tol f df initVec diff --git a/src/AI/Training/Internal/lbfgs_aux.c b/src/AI/Training/Internal/lbfgs_aux.c new file mode 100644 --- /dev/null +++ b/src/AI/Training/Internal/lbfgs_aux.c @@ -0,0 +1,85 @@ +//_______________________________________________________________________________ +// +// Author: Kiet Lam +// File lbfgs_aux.c +//_______________________________________________________________________________ +// Last Updated: Time-stamp: <2012-01-18 23:41:52 (lam)> +// +// +// 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 . + + +#include +#include + + +#define HASKELLARRAY(A) int A##n, double* A##p + + +typedef double (*fFunc) (int, const lbfgsfloatval_t*); +typedef int (*dfFunc) (int, const lbfgsfloatval_t*, int, lbfgsfloatval_t*); + + +typedef struct +{ + double (*f) (int, const lbfgsfloatval_t*); + int (*df) (int, const lbfgsfloatval_t*, int, double*); +} FdfData; + + +lbfgsfloatval_t lbfgs_evaluate_aux (void* data, + const lbfgsfloatval_t* x, + lbfgsfloatval_t* g, + const int n, + const lbfgsfloatval_t step) +{ + FdfData* fdf = (FdfData*) data; + fdf->df(n, x, n, g); + return fdf->f(n, x); +} + + + +int minimizeLBFGS (double precision, int max_iter, double init_step, double tol, + fFunc fun, dfFunc dfun, HASKELLARRAY(x), HASKELLARRAY(r)) +{ + FdfData fdfDat; + fdfDat.f = fun; + fdfDat.df = dfun; + + lbfgs_parameter_t param; + + lbfgs_parameter_init(¶m); + + param.epsilon = precision; + param.max_iterations = max_iter; + + lbfgsfloatval_t* v = lbfgs_malloc(xn); + + int i; + for (i = 0; i < xn; ++i) + { + v[i] = xp[i]; + } + + lbfgs(xn, v, NULL, lbfgs_evaluate_aux, NULL, &fdfDat, ¶m); + + int j; + for (j = 0; j < xn; ++j) + { + rp[j] = v[j]; + } + + lbfgs_free(v); +}