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. + + {one line to give the program's name and a brief idea of what it does.} + Copyright (C) {year} {name of author} + + 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: + + {project} Copyright (C) {year} {fullname} + 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/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,57 @@ +# sparse-linear-algebra + +Sparse linear algebra datastructures and algorithms in Haskell + +TravisCI : [![Build Status](https://travis-ci.org/ocramz/sparse-linear-algebra.png)](https://travis-ci.org/ocramz/sparse-linear-algebra) + +This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting. + +Algorithms : + +* Iterative linear solvers + + * Conjugate Gradient Squared (CGS) + + * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems) + +* Matrix decompositions + + * QR factorization + +* Eigenvalue algorithms + + * QR algorithm + + * Rayleigh quotient iteration + +* Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection + +* Predicates : Matrix orthogonality test (A^T A ~= I) + + + +---------- + +This is also an experiment in principled scientific programming : + +* set the stage by declaring typeclasses and some useful generic operations (normed linear vector spaces, i.e. finite-dimensional spaces equipped with an inner product that induces a distance function), + +* define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via `Data.IntMap`, are made instances of the VectorSpace and AdditiveGroup classes respectively). This allows to decouple the algorithms from the actual implementation of the backend, + +* implement the algorithms, following 1:1 the textbook [1] + + +## License + +GPL3, see LICENSE + +## Credits + +Inspired by + +* `linear` : https://hackage.haskell.org/package/linear +* `sparse-lin-alg` : https://github.com/laughedelic/sparse-lin-alg + +## References + +[1] : Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000 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/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,8 @@ +module Main where + +import Lib (ourAdd) + +import Text.Printf (printf) + +main :: IO () +main = printf "2 + 3 = %d\n" (ourAdd 2 3) diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal new file mode 100644 --- /dev/null +++ b/sparse-linear-algebra.cabal @@ -0,0 +1,66 @@ +name: sparse-linear-algebra +version: 0.1.0.0 +synopsis: Sparse linear algebra datastructures and algorithms +description: Please see README.md +homepage: https://github.com/ocramz/sparse-linear-algebra +license: BSD3 +license-file: LICENSE +author: Marco Zocca +maintainer: zocca.marco gmail +copyright: 2016 Marco Zocca +category: Math +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 +tested-with: GHC == 8.0.1 + +library + default-language: Haskell2010 + ghc-options: -Wall + hs-source-dirs: src + exposed-modules: Lib + Math.Linear.Sparse + Math.Linear.Sparse.IntMap + build-depends: QuickCheck + , base >= 4.7 && < 5 + , containers + , hspec + , primitive >= 0.6.1.0 + , transformers >= 0.5.2.0 + -- , lens + , mtl >= 2.2.1 + , mwc-random + , monad-loops + +executable sparse-linear-algebra + default-language: Haskell2010 + ghc-options: -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: app + main-is: Main.hs + build-depends: base + , mtl >= 2.2.1 + , mwc-random + , primitive >= 0.6.1.0 + , sparse-linear-algebra + , transformers >= 0.5.2.0 + +test-suite spec + default-language: Haskell2010 + ghc-options: -Wall + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Spec.hs + build-depends: base + , containers + , hspec + , mtl >= 2.2.1 + , mwc-random + , primitive >= 0.6.1.0 + , sparse-linear-algebra + , transformers >= 0.5.2.0 + , criterion + -- , QuickCheck + +source-repository head + type: git + location: https://github.com/ocramz/sparse-linear-algebra diff --git a/src/Lib.hs b/src/Lib.hs new file mode 100644 --- /dev/null +++ b/src/Lib.hs @@ -0,0 +1,11 @@ +-- | A library to do stuff. +module Lib + ( + ourAdd + ) where + +-- | Add two 'Int' values. +ourAdd :: Int -- ^ left + -> Int -- ^ right + -> Int -- ^ sum +ourAdd x y = x + y diff --git a/src/Math/Linear/Sparse.hs b/src/Math/Linear/Sparse.hs new file mode 100644 --- /dev/null +++ b/src/Math/Linear/Sparse.hs @@ -0,0 +1,1842 @@ +{-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-} +-- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-} + +module Math.Linear.Sparse where + +import Math.Linear.Sparse.IntMap + + +import Control.Monad.Primitive + +import Control.Monad (mapM_, forM_, replicateM) +import Control.Monad.Loops + +import Control.Monad.Cont +import Control.Monad.State.Strict +import Control.Monad.Writer +import Control.Monad.Trans + +import Control.Monad.Trans.State (runStateT) +import Control.Monad.Trans.Writer (runWriterT) + +import qualified Data.IntMap.Strict as IM +-- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers` + +import qualified System.Random.MWC as MWC +import qualified System.Random.MWC.Distributions as MWC + +import Data.Monoid +import qualified Data.Foldable as F +import qualified Data.Traversable as T + + +import qualified Data.List as L +import Data.Maybe + + + +-- | ========= CLASSES and common operations + +-- | Additive ring +class Functor f => Additive f where + -- | zero element + zero :: Num a => f a + + -- | componentwise operations + (^+^) :: Num a => f a -> f a -> f a + (^-^) :: Num a => f a -> f a -> f a + + + +-- | negate the values in a functor +negated :: (Num a, Functor f) => f a -> f a +negated = fmap negate + +x `minus` y = x ^+^ negated y + + + +-- | Vector space +class Additive f => VectorSpace f where + -- | multiplication by a scalar + (.*) :: Num a => a -> f a -> f a + + +-- |linear interpolation +lerp :: (VectorSpace f, Num a) => a -> f a -> f a -> f a +lerp a u v = a .* u ^+^ ((1-a) .* v) + + +-- | Hilbert space (inner product) +class VectorSpace f => Hilbert f where + -- | inner product + dot :: Num a => f a -> f a -> a + +class Hilbert f => Normed f where + norm :: (Floating a, Eq a) => a -> f a -> a + + +-- some norms and related results + +-- squared norm +normSq :: (Hilbert f, Num a) => f a -> a +normSq v = v `dot` v + + +-- L1 norm +norm1 :: (Foldable t, Num a, Functor t) => t a -> a +norm1 v = sum (fmap abs v) + +-- Euclidean norm +norm2 :: (Hilbert f, Floating a) => f a -> a +norm2 v = sqrt (normSq v) + +-- Lp norm (p > 0) +normP :: (Foldable t, Functor t, Floating a) => a -> t a -> a +normP p v = sum u**(1/p) where + u = fmap (**p) v + +-- infinity-norm +normInfty :: (Foldable t, Ord a) => t a -> a +normInfty = maximum + + + +-- normalize +normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a +normalize n v = (1 / norm n v) .* v + + + + + + + +-- -- Lp inner product (p > 0) +dotLp :: (Set t, Foldable t, Floating a) => a -> t a -> t a -> a +dotLp p v1 v2 = sum u**(1/p) where + f a b = (a*b)**p + u = liftI2 f v1 v2 + + +-- reciprocal +reciprocal :: (Functor f, Fractional b) => f b -> f b +reciprocal = fmap recip + + +-- scale +scale :: (Num b, Functor f) => b -> f b -> f b +scale n = fmap (* n) + + + + + + + +-- | FiniteDim : finite-dimensional objects + +class Additive f => FiniteDim f where + type FDSize f :: * + dim :: f a -> FDSize f + + +instance FiniteDim SpVector where + type FDSize SpVector = Int + dim = svDim + + +instance FiniteDim SpMatrix where + type FDSize SpMatrix = (Rows, Cols) + dim = smDim + + +-- unary dimension-checking bracket +withDim :: (FiniteDim f, Show e) => + f a + -> (FDSize f -> f a -> Bool) + -> (f a -> c) + -> String + -> (f a -> e) + -> c +withDim x p f e ef | p (dim x) x = f x + | otherwise = error e' where e' = e ++ show (ef x) + +-- binary dimension-checking bracket +withDim2 :: (FiniteDim f, FiniteDim g, Show e) => + f a + -> g b + -> (FDSize f -> FDSize g -> f a -> g b -> Bool) + -> (f a -> g b -> c) + -> String + -> (f a -> g b -> e) + -> c +withDim2 x y p f e ef | p (dim x) (dim y) x y = f x y + | otherwise = error e' where e' = e ++ show (ef x y) + + + + + + +-- | HasData : accessing inner data (do not export) + +class Additive f => HasData f a where + type HDData f a :: * + dat :: f a -> HDData f a + +instance HasData SpVector a where + type HDData SpVector a = IM.IntMap a + dat = svData + +instance HasData SpMatrix a where + type HDData SpMatrix a = IM.IntMap (IM.IntMap a) + dat = smData + + + + + + +-- | Sparse : sparse datastructures + +class (FiniteDim f, HasData f a) => Sparse f a where + spy :: Fractional b => f a -> b + + +instance Sparse SpVector a where + spy = spySV + +instance Sparse SpMatrix a where + spy = spySM + + + + +class Functor f => Set f where + -- |union binary lift + liftU2 :: (a -> a -> a) -> f a -> f a -> f a + + -- |intersection binary lift + liftI2 :: (a -> b -> c) -> f a -> f b -> f c + + + +-- class (Set f, Sparse f a) => SparseSet f a + +-- instance SparseSet SpVector a where + + + +-- | ======================================================= + +instance Set IM.IntMap where + liftU2 = IM.unionWith + {-# INLINE liftU2 #-} + liftI2 = IM.intersectionWith + {-# INLINE liftI2 #-} + +-- | IntMap implementation +instance Additive IM.IntMap where + zero = IM.empty + {-# INLINE zero #-} + (^+^) = liftU2 (+) + {-# INLINE (^+^) #-} + x ^-^ y = x ^+^ negated y + {-# INLINE (^-^) #-} + +instance VectorSpace IM.IntMap where + n .* im = IM.map (* n) im + +instance Hilbert IM.IntMap where + a `dot` b = sum $ liftI2 (*) a b + + +instance Normed IM.IntMap where + norm p v | p==1 = norm1 v + | p==2 = norm2 v + | otherwise = normP p v + + + +-- | ======================================================= + +-- | Sparse Vector +data SpVector a = SV { svDim :: Int , + svData :: IM.IntMap a} deriving Eq + +dimSV :: SpVector a -> Int +dimSV = svDim + +spySV :: Fractional b => SpVector a -> b +spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (svDim s) + + +-- internal : projection functions, do not export +imSV :: SpVector a -> IM.IntMap a +imSV = svData + + + + +-- | instances for SpVector +instance Functor SpVector where + fmap f (SV n x) = SV n (fmap f x) + +instance Set SpVector where + liftU2 f2 (SV n1 x1) (SV n2 x2) = SV (max n1 n2) (liftU2 f2 x1 x2) + liftI2 f2 (SV n1 x1) (SV n2 x2) = SV (max n1 n2) (liftI2 f2 x1 x2) + +instance Foldable SpVector where + foldr f d v = F.foldr f d (svData v) + +instance Additive SpVector where + zero = SV 0 IM.empty + (^+^) = liftU2 (+) + (^-^) = liftU2 (-) + + +instance VectorSpace SpVector where + n .* v = scale n v + +instance Hilbert SpVector where + a `dot` b | dim a == dim b = dot (dat a) (dat b) + | otherwise = + error $ "dot : sizes must coincide, instead we got " ++ + show (dim a, dim b) + + +instance Normed SpVector where + norm p (SV _ v) = norm p v + + + + + + +-- | empty sparse vector (size n, no entries) + +zeroSV :: Int -> SpVector a +zeroSV n = SV n IM.empty + + +singletonSV :: a -> SpVector a +singletonSV x = SV 1 (IM.singleton 0 x) + + +-- | create a sparse vector from an association list while discarding all zero entries +mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a +mkSpVector d im = SV d $ IM.filterWithKey (\k v -> v /= 0 && inBounds0 d k) im + +-- | ", from logically dense array (consecutive indices) +mkSpVectorD :: (Num a, Eq a) => Int -> [a] -> SpVector a +mkSpVectorD d ll = mkSpVector d (IM.fromList $ denseIxArray (take d ll)) + +-- ", don't filter zero elements +mkSpVector1 :: Int -> IM.IntMap a -> SpVector a +mkSpVector1 d ll = SV d $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll + +mkSpVector1D :: Int -> [a] -> SpVector a +mkSpVector1D d ll = mkSpVector1 d (IM.fromList $ denseIxArray (take d ll)) + + + +-- | DENSE vector of `1`s +onesSV :: Num a => Int -> SpVector a +onesSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 1 + +-- | DENSE vector of `0`s +zerosSV :: Num a => Int -> SpVector a +zerosSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 0 + + + + +-- insert +insertSpVector :: Int -> a -> SpVector a -> SpVector a +insertSpVector i x (SV d xim) + | inBounds0 d i = SV d (IM.insert i x xim) + | otherwise = error "insertSpVector : index out of bounds" + + +fromListSV :: Int -> [(Int, a)] -> SpVector a +fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix )) + +-- toList +toListSV :: SpVector a -> [(IM.Key, a)] +toListSV sv = IM.toList (imSV sv) + +-- to dense list (default = 0) +toDenseListSV :: Num b => SpVector b -> [b] +toDenseListSV (SV d im) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1] + + + + + + + + + + +instance Show a => Show (SpVector a) where + show (SV d x) = "SV (" ++ show d ++ ") "++ show (IM.toList x) + + +-- | lookup + +lookupDenseSV :: Num a => IM.Key -> SpVector a -> a +lookupDenseSV i (SV _ im) = IM.findWithDefault 0 i im + +findWithDefault0IM :: Num a => IM.Key -> IM.IntMap a -> a +findWithDefault0IM = IM.findWithDefault 0 + + + + +-- | SV manipulation + +tailSV :: SpVector a -> SpVector a +tailSV (SV n sv) = SV (n-1) ta where + ta = IM.mapKeys (\i -> i - 1) $ IM.delete 0 sv + + +headSV :: Num a => SpVector a -> a +headSV sv = fromMaybe 0 (IM.lookup 0 (imSV sv)) + + + +-- | concatenate SpVector + + +concatSV :: SpVector a -> SpVector a -> SpVector a +concatSV (SV n1 s1) (SV n2 s2) = SV (n1+n2) (IM.union s1 s2') where + s2' = IM.mapKeys (+ n1) s2 + + + + + + + + + + +-- | promote a SV to SM + +svToSM :: SpVector a -> SpMatrix a +svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d + + + + + +-- | outer vector product + +outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a +outerProdSV v1 v2 = fromListSM (m, n) ixy where + m = dim v1 + n = dim v2 + ixy = [(i,j, x * y) | (i,x) <- toListSV v1 , (j, y) <- toListSV v2] + +(><) = outerProdSV + + + + + + + + + +-- | ======================================================= + + + + +data SpMatrix a = SM {smDim :: (Rows, Cols), + smData :: IM.IntMap (IM.IntMap a)} deriving Eq + + + +-- | instances for SpMatrix +instance Show a => Show (SpMatrix a) where + show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x) + +instance Functor SpMatrix where + fmap f (SM d md) = SM d ((fmap . fmap) f md) + +instance Set SpMatrix where + liftU2 f2 (SM n1 x1) (SM n2 x2) = SM (maxTup n1 n2) ((liftU2.liftU2) f2 x1 x2) + liftI2 f2 (SM n1 x1) (SM n2 x2) = SM (minTup n1 n2) ((liftI2.liftI2) f2 x1 x2) + +instance Additive SpMatrix where + zero = SM (0,0) IM.empty + (^+^) = liftU2 (+) + (^-^) = liftU2 (-) + + +-- | TODO : use semilattice properties instead +maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t) +maxTup (x1,y1) (x2,y2) = (max x1 x2, max y1 y2) +minTup (x1,y1) (x2,y2) = (min x1 x2, min y1 y2) + +-- | empty matrix of size d +emptySpMatrix :: (Int, Int) -> SpMatrix a +emptySpMatrix d = SM d IM.empty + + + +-- multiply matrix by a scalar +matScale :: Num a => a -> SpMatrix a -> SpMatrix a +matScale a = fmap (*a) + +-- Frobenius norm (sqrt of trace of M^T M) +normFrobenius :: SpMatrix Double -> Double +normFrobenius m = sqrt $ foldlSM (+) 0 m' where + m' | nrows m > ncols m = transposeSM m ## m + | otherwise = m ## transposeSM m + + + + + + +-- | ========= MATRIX METADATA + +-- type synonyms +type Rows = Int +type Cols = Int + +type IxRow = Int +type IxCol = Int + +-- -- predicates +-- are the supplied indices within matrix bounds? +validIxSM :: SpMatrix a -> (Int, Int) -> Bool +validIxSM mm = inBounds02 (dim mm) + +-- is the matrix square? +isSquareSM :: SpMatrix a -> Bool +isSquareSM m = nrows m == ncols m + +-- is the matrix diagonal? +isDiagonalSM :: SpMatrix a -> Bool +isDiagonalSM m = IM.size d == nrows m where + d = IM.filterWithKey ff (immSM m) + ff irow row = IM.size row == 1 && + IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1 + + + + + + + + +-- -- internal projection functions, do not export: +immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t) +immSM (SM _ imm) = imm + +dimSM :: SpMatrix t -> (Rows, Cols) +dimSM (SM d _) = d + +nelSM :: SpMatrix t -> Int +nelSM (SM (nr,nc) _) = nr*nc + +-- | nrows, ncols : size accessors +nrows, ncols :: SpMatrix a -> Int +nrows = fst . dim +ncols = snd . dim + + + + +data SMInfo = SMInfo { smNz :: Int, + smSpy :: Double} deriving (Eq, Show) + +infoSM :: SpMatrix a -> SMInfo +infoSM s = SMInfo (nzSM s) (spySM s) + +nzSM :: SpMatrix a -> Int +nzSM s = sum $ fmap IM.size (immSM s) + +spySM :: Fractional b => SpMatrix a -> b +spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s) + + +-- # NZ in row i + +nzRowU :: SpMatrix a -> IM.Key -> Int +nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s) + +nzRow :: SpMatrix a -> IM.Key -> Int +nzRow s i | inBounds0 (nrows s) i = nzRowU s i + | otherwise = error "nzRow : index out of bounds" + + + + +-- | bandwidth bounds (min, max) + +bwMinSM :: SpMatrix a -> Int +bwMinSM = fst . bwBoundsSM + +bwMaxSM :: SpMatrix a -> Int +bwMaxSM = snd . bwBoundsSM + +bwBoundsSM :: SpMatrix a -> (Int, Int) +bwBoundsSM s = -- b + (snd $ IM.findMin b, + snd $ IM.findMax b) + where + ss = immSM s + fmi = fst . IM.findMin + fma = fst . IM.findMax + b = fmap (\x -> fma x - fmi x + 1:: Int) ss + + + + + + + + + +-- | ========= SPARSE MATRIX BUILDERS + +zeroSM :: Int -> Int -> SpMatrix a +zeroSM m n = SM (m,n) IM.empty + + +insertSpMatrix :: IxRow -> IxCol -> a -> SpMatrix a -> SpMatrix a +insertSpMatrix i j x s + | inBounds02 d (i,j) = SM d $ insertIM2 i j x smd + | otherwise = error "insertSpMatrix : index out of bounds" where + smd = immSM s + d = dim s + + +-- | from list (row, col, value) +fromListSM' :: Foldable t => t (IxRow, IxCol, a) -> SpMatrix a -> SpMatrix a +fromListSM' iix sm = foldl ins sm iix where + ins t (i,j,x) = insertSpMatrix i j x t + +fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a +fromListSM (m,n) iix = fromListSM' iix (zeroSM m n) + + +fromListDenseSM :: Int -> [a] -> SpMatrix a +fromListDenseSM m ll = fromListSM (m, n) $ denseIxArray2 m ll where + n = length ll `div` m + + + +-- | to List + +-- toDenseListSM : populate missing entries with 0 +toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)] +toDenseListSM m = + [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]] + + + + + +-- -- create diagonal and identity matrix +mkDiagonal :: Int -> [a] -> SpMatrix a +mkDiagonal n = mkSubDiagonal n 0 + + +eye :: Num a => Int -> SpMatrix a +eye n = mkDiagonal n (ones n) + +ones :: Num a => Int -> [a] +ones n = replicate n 1 + + + + + +-- super- and sub- diagonal + +mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a +mkSubDiagonal n o xx | abs o < n = if o >= 0 + then fz ii jj xx + else fz jj ii xx + | otherwise = error "mkSubDiagonal : offset > dimension" where + ii = [0 .. n-1] + jj = [abs o .. n - 1] + fz a b x = fromListSM (n,n) (zip3 a b x) + + +-- fromList :: [(Key,a)] -> IntMap a +-- fromList xs +-- = foldlStrict ins empty xs +-- where +-- ins t (k,x) = insert k x t + + + + +encode :: (Int, Int) -> (Rows, Cols) -> Int +encode (nr,_) (i,j) = i + (j * nr) + +decode :: (Int, Int) -> Int -> (Rows, Cols) +decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr + + + + + + + + + +-- | ========= SUB-MATRICES + + +extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a +extractSubmatrixSM (SM (r, c) im) (i1, i2) (j1, j2) + | q = SM (m', n') imm' + | otherwise = error $ "extractSubmatrixSM : invalid indexing " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where + imm' = mapKeysIM2 (\i -> i - i1) (\j -> j - j1) $ -- rebalance keys + IM.filter (not . IM.null) $ -- remove all-null rows + ifilterIM2 ff im -- keep `submatrix` + ff i j _ = i1 <= i && + i <= i2 && + j1 <= j && + j <= j2 + (m', n') = (i2-i1 + 1, j2-j1 + 1) + q = inBounds0 r i1 && + inBounds0 r i2 && + inBounds0 c j1 && + inBounds0 c j2 && + i2 >= i1 + +-- extract row / column +extractRowSM :: SpMatrix a -> Int -> SpMatrix a +extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1) + +extractColSM :: SpMatrix a -> Int -> SpMatrix a +extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j) + + + +-- demote (n x 1) or (1 x n) SpMatrix to SpVector +toSV :: SpMatrix a -> SpVector a +toSV (SM (m,n) im) = SV d $ snd . head $ IM.toList im where + d | m==1 && n==1 = 1 + | m==1 && n>1 = n + | n==1 && m>1 = m + | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n) + + +-- extract row or column and place into SpVector +extractCol :: SpMatrix a -> Int -> SpVector a +extractCol m i = toSV $ extractColSM m i + +extractRow :: SpMatrix a -> Int -> SpVector a +extractRow m j = toSV $ extractRowSM m j + + + + + + + +-- | ========= MATRIX STACKING + +vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a +vertStackSM mm1 mm2 = SM (m, n) $ IM.union u1 u2 where + nro1 = nrows mm1 + m = nro1 + nrows mm2 + n = max (ncols mm1) (ncols mm2) + u1 = immSM mm1 + u2 = IM.mapKeys (+ nro1) (immSM mm2) + +(-=-) = vertStackSM + + +horizStackSM, (-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a +horizStackSM mm1 mm2 = t (t mm1 -=- t mm2) where + t = transposeSM + +(-||-) = horizStackSM + + + + + + + + + +-- | ========= LOOKUP + +lookupSM :: SpMatrix a -> IM.Key -> IM.Key -> Maybe a +lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j + +-- | Looks up an element in the matrix (if not found, zero is returned) + +lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IM.Key, IM.Key) -> a +lookupWD_SM sm (i,j) = + fromMaybe 0 (lookupSM sm i j) + +lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IM.Key, IM.Key) -> a +lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) + +(@@) = lookupWD_SM + + + + +-- FIXME : to throw an exception or just ignore the out-of-bound access ? + + + + + + + +-- | ========= MISC SpMatrix OPERATIONS + +foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b +foldlSM f n (SM _ m)= foldlIM2 f n m + +ifoldlSM :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> SpMatrix a -> b +ifoldlSM f n (SM _ m) = ifoldlIM2' f n m + + + + + + +-- | mapping + + + + + +-- | folding + +-- count sub-diagonal nonzeros +countSubdiagonalNZSM :: SpMatrix a -> Int +countSubdiagonalNZSM (SM _ im) = countSubdiagonalNZ im + + +-- | filtering + +-- extractDiagonalSM :: (Num a, Eq a) => SpMatrix a -> SpVector a +-- extractDiagonalSM (SM (m,n) im) = mkSpVectorD m $ extractDiagonalIM2 im + +-- extract with default 0 +extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a +extractDiagonalDSM mm = mkSpVector1D n $ foldr ins [] ll where + ll = [0 .. n - 1] + n = nrows mm + ins i acc = mm@@(i,i) : acc + + + + + + + +-- filtering the index subset that lies below the diagonal + +subdiagIndicesSM :: SpMatrix a -> [(IM.Key, IM.Key)] +subdiagIndicesSM (SM _ im) = subdiagIndices im + + + + + +-- | sparsify : remove 0s (!!!) + +sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double) +sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps) + +sparsifySM :: SpMatrix Double -> SpMatrix Double +sparsifySM (SM d im) = SM d $ sparsifyIM2 im + + + +-- | ROUNDING operations (!!!) + +roundZeroOneSM :: SpMatrix Double -> SpMatrix Double +roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im + + + + + + + + + + + +-- | ========= ALGEBRAIC PRIMITIVE OPERATIONS + + +-- | transpose + + +transposeSM, (#^) :: SpMatrix a -> SpMatrix a +transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im) + +(#^) = transposeSM + + + +-- | A^T B +(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double +a #^# b = transposeSM a #~# b + + +-- | A B^T +(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double +a ##^ b = a #~# transposeSM b + + + + + + +-- | matrix action on a vector + +{- +FIXME : matVec is more general than SpVector's : + +\m v -> fmap (`dot` v) m + :: (Normed f1, Num b, Functor f) => f (f1 b) -> f1 b -> f b +-} + + + +-- matrix on vector +matVec, (#>) :: Num a => SpMatrix a -> SpVector a -> SpVector a +matVec (SM (nr, nc) mdata) (SV n sv) + | nc == n = SV nr $ fmap (`dot` sv) mdata + | otherwise = error $ "matVec : mismatching dimensions " ++ show (nc, n) + +(#>) = matVec + +-- vector on matrix (FIXME : transposes matrix: more costly than `matVec`) +vecMat, (<#) :: Num a => SpVector a -> SpMatrix a -> SpVector a +vecMat (SV n sv) (SM (nr, nc) mdata) + | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata) + | otherwise = error $ "vecMat : mismatching dimensions " ++ show (n, nr) + +(<#) = vecMat + + + + + + + + + + +-- matVec' mm vv = +-- withDim2 mm vv (\(nro, nco) nv _ _ -> nco == nv) matVecU "matVec : mismatching dimensions" +-- (\ m v -> unwords [show (dim m), show (dim v)]) + + +-- asdfm ll = unwords (map (show . dim) ll) + + + + +-- | matrix-matrix product + +-- unsafe matMat +matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a +matMatU m1 m2 = + SM (nrows m1, ncols m2) im where + im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1) + + +-- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a +-- matMat (SM (nr1,nc1) m1) (SM (nr2,nc2) m2) +-- | nc1 == nr2 = SM (nr1, nc2) $ +-- fmap (\vm1 -> fmap (`dot` vm1) (transposeIM2 m2)) m1 +-- | otherwise = error "matMat : incompatible matrix sizes" + +matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a +matMat m1 m2 + | c1 == r2 = matMatU m1 m2 + | otherwise = error $ "matMat : incompatible matrix sizes" ++ show (d1, d2) where + d1@(r1, c1) = dim m1 + d2@(r2, c2) = dim m2 + + +(##) = matMat + +-- matMat m1 m2 = +-- withDim2 m1 m2 +-- (\(r1,c1) (r2,c2) _ _ -> c1 == r2) +-- matMatU +-- "matMat : incompatible matrix sizes" +-- (\m1 m2 -> unwords [show (dim m1), show (dim m2)]) + + + + + +-- | sparsified matrix-matrix product (prunes all elements `x` for which `abs x <= eps`) +matMatSparsified, (#~#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double +matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2 + +(#~#) = matMatSparsified + + + + + + +-- | ========= predicates + +-- is the matrix orthogonal? i.e. Q^t ## Q == I +isOrthogonalSM :: SpMatrix Double -> Bool +isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where + rsm = roundZeroOneSM $ transposeSM sm ## sm + + + + + + + + + +-- | ========= condition number + +-- uses the R matrix from the QR factorization +conditionNumberSM :: SpMatrix Double -> Double +conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system" + | otherwise = kappa where + kappa = lmax / lmin + (_, r) = qr m + u = extractDiagonalDSM r -- FIXME : need to extract with default element 0 + lmax = abs (maximum u) + lmin = abs (minimum u) + + + + + + + +-- | ========= Householder transformation + +hhMat :: Num a => a -> SpVector a -> SpMatrix a +hhMat beta x = eye n ^-^ scale beta (x >< x) where + n = dim x + + +-- a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane: +-- v' = (I - 2 x >< x) v +hhRefl :: SpVector Double -> SpMatrix Double +hhRefl = hhMat 2.0 + + + + + + + + + + + +-- | ========= Givens rotation matrix + + +hypot :: Floating a => a -> a -> a +hypot x y = abs x * (sqrt (1 + y/x)**2) + +sign :: (Ord a, Num a) => a -> a +sign x + | x > 0 = 1 + | x == 0 = 0 + | otherwise = -1 + +givensCoef :: (Ord a, Floating a) => a -> a -> (a, a, a) +givensCoef a b -- returns (c, s, r) where r = norm (a, b) + | b==0 = (sign a, 0, abs a) + | a==0 = (0, sign b, abs b) + | abs a > abs b = let t = b/a + u = sign a * abs ( sqrt (1 + t**2)) + in (1/u, - t/u, a*u) + | otherwise = let t = a/b + u = sign b * abs ( sqrt (1 + t**2)) + in (t/u, - 1/u, b*u) + + +{- +Givens method, row version: choose other row index i' s.t. i' is : +* below the diagonal +* corresponding element is nonzero + +QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is +non-zero but A has zeros in row k for all columns less than j. +-} + +givens :: SpMatrix Double -> Int -> Int -> SpMatrix Double +givens mm i j + | validIxSM mm (i,j) && isSquareSM mm = + sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm)) + | otherwise = error "givens : indices out of bounds" + where + (c, s, _) = givensCoef a b + i' = head $ fromMaybe (error $ "givens: no compatible rows for entry " ++ show (i,j)) (candidateRows (immSM mm) i j) + a = mm @@ (i', j) + b = mm @@ (i, j) -- element to zero out + +-- is the `k`th the first nonzero column in the row? +firstNonZeroColumn :: IM.IntMap a -> IM.Key -> Bool +firstNonZeroColumn mm k = isJust (IM.lookup k mm) && + isNothing (IM.lookupLT k mm) + +-- returns a set of rows {k} that satisfy QR.C1 +candidateRows :: IM.IntMap (IM.IntMap a) -> IM.Key -> IM.Key -> Maybe [IM.Key] +candidateRows mm i j | IM.null u = Nothing + | otherwise = Just (IM.keys u) where + u = IM.filterWithKey (\irow row -> irow /= i && + firstNonZeroColumn row j) mm + + + + + +-- | ========= QR algorithm + +{- +applies Givens rotation iteratively to zero out sub-diagonal elements +-} + + +qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double) +qr mm = (transposeSM qmatt, rmat) where + qmatt = F.foldl' (#~#) ee $ gmats mm -- Q^T = (G_n * G_n-1 ... * G_1) + rmat = qmatt #~# mm -- R = Q^T A + ee = eye (nrows mm) + +-- Givens matrices in order [G1, G2, .. , G_N ] +gmats :: SpMatrix Double -> [SpMatrix Double] +gmats mm = gm mm (subdiagIndicesSM mm) where + gm m ((i,j):is) = let g = givens m i j + in g : gm (g #~# m) is + gm _ [] = [] + + + + + +-- -- | QR algorithm, state transformer version +-- gmatST0 (m, (i,j):is) = (m', is) where -- WRONG, possible access to [] +-- g = givens m i j +-- m' = g #~# m +-- gmatST0 (m, []) = (eye (nrows m), []) + +-- gmatST m = gmatST0 (m, subdiagIndicesSM m) + + + + + +-- | ========= Eigenvalues, using QR + + +eigsQR :: Int -> SpMatrix Double -> SpVector Double +eigsQR nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where + eigsStep m = r #~# q where (q, r) = qr m + convergtest g = modifyInspectN nitermax f g where + f [m1, m2] = let dm1 = extractDiagonalDSM m1 + dm2 = extractDiagonalDSM m2 + in norm2 (dm1 ^-^ dm2) <= eps + + + + + + +-- | ========= Eigenvalues, using Rayleigh iteration + +rayleighStep :: + SpMatrix Double -> + (SpVector Double, Double) -> + (SpVector Double, Double) -- updated estimate of (eigenvector, eigenvalue) +rayleighStep aa (b, mu) = (b', mu') where + ii = eye (nrows aa) + nom = (aa ^-^ (mu `matScale` ii)) <\> b + b' = normalize 2 nom + mu' = b' `dot` (aa #> b') / (b' `dot` b') + +eigRayleigh :: Int -- max # iterations + -> SpMatrix Double -- matrix + -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue) + -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue) +eigRayleigh nitermax m = execState (convergtest (rayleighStep m)) where + convergtest g = modifyInspectN nitermax f g where + f [(b1, _), (b2, _)] = norm2 (b2 ^-^ b1) <= eps + + + + + +-- | ========= Householder vector (G & VL Alg. 5.1.1, function `house`) + +-- hhV :: (Ord a, Floating a) => SpVector a -> (SpVector a, a) +hhV :: SpVector Double -> (SpVector Double, Double) +hhV x = (v, beta) where + n = dim x + tx = tailSV x + sigma = tx `dot` tx + vtemp = singletonSV 1 `concatSV` tx + (v, beta) | sigma <= eps = (vtemp, 0) + | otherwise = let mu = sqrt (headSV x**2 + sigma) + xh = headSV x + vh | xh <= 1 = xh - mu + | otherwise = - sigma / (xh + mu) + vnew = (1 / vh) .* insertSpVector 0 vh vtemp + in (vnew, 2 * xh**2 / (sigma + vh**2)) + + + + + + + + +-- | ========= SVD + +{- Golub & Van Loan, sec 8.6.2 (p 452 segg.) + +SVD of A : + +* reduce A to upper bidiagonal form B (Alg. 5.4.2) +* compute SVD of B (implicit-shift QR step, Alg. 8.3.2) + +-} + + + + + + + + + + + + + + + + + + +-- | ======================================================= + +-- | LINEAR SOLVERS : solve A x = b + +-- | numerical tolerance for e.g. solution convergence +eps :: Double +eps = 1e-8 + +-- | residual of candidate solution x0 +residual :: Num a => SpMatrix a -> SpVector a -> SpVector a -> SpVector a +residual aa b x0 = b ^-^ (aa #> x0) + +converged :: SpMatrix Double -> SpVector Double -> SpVector Double -> Bool +converged aa b x0 = normSq (residual aa b x0) <= eps + + + +-- | CGS + +-- | one step of CGS +cgsStep :: SpMatrix Double -> SpVector Double -> CGS -> CGS +cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1 + where + aap = aa #> p + alphaj = (r `dot` rhat) / (aap `dot` rhat) + q = u ^-^ (alphaj .* aap) + xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution + rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q)))-- updated residual + betaj = (rj1 `dot` rhat) / (r `dot` rhat) + uj1 = rj1 ^+^ (betaj .* q) + pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p))) + +data CGS = CGS { _x :: SpVector Double, + _r :: SpVector Double, + _p :: SpVector Double, + _u :: SpVector Double } deriving Eq + +cgs :: + SpMatrix Double -> + SpVector Double -> + SpVector Double -> + SpVector Double -> + -- Int -> + CGS +cgs aa b x0 rhat = + -- execState (replicateM n (modify (cgsStep aa rhat))) cgsInit where + execState (untilConverged _x (cgsStep aa rhat)) cgsInit where + r0 = b ^-^ (aa #> x0) -- residual of initial guess solution + p0 = r0 + u0 = r0 + cgsInit = CGS x0 r0 p0 u0 + + +instance Show CGS where + show (CGS x r p u) = "x = " ++ show x ++ "\n" ++ + "r = " ++ show r ++ "\n" ++ + "p = " ++ show p ++ "\n" ++ + "u = " ++ show u ++ "\n" + + + + + +-- | BiCSSTAB + +-- _aa :: SpMatrix Double, -- matrix +-- _b :: SpVector Double, -- rhs +-- _r0 :: SpVector Double, -- initial residual +-- _r0hat :: SpVector Double, -- candidate solution: r0hat `dot` r0 >= 0 + +-- | one step of BiCGSTAB +bicgstabStep :: SpMatrix Double -> SpVector Double -> BICGSTAB -> BICGSTAB +bicgstabStep aa r0hat (BICGSTAB x r p) = BICGSTAB xj1 rj1 pj1 where + aap = aa #> p + alphaj = (r `dot` r0hat) / (aap `dot` r0hat) + sj = r ^-^ (alphaj .* aap) + aasj = aa #> sj + omegaj = (aasj `dot` sj) / (aasj `dot` aasj) + xj1 = x ^+^ (alphaj .* p) ^+^ (omegaj .* sj) + rj1 = sj ^-^ (omegaj .* aasj) + betaj = (rj1 `dot` r0hat)/(r `dot` r0hat) * alphaj / omegaj + pj1 = rj1 ^+^ (betaj .* (p ^-^ (omegaj .* aap))) + +data BICGSTAB = BICGSTAB { _xBicgstab :: SpVector Double, + _rBicgstab :: SpVector Double, + _pBicgstab :: SpVector Double} deriving Eq + + +bicgstab + :: SpMatrix Double + -> SpVector Double + -> SpVector Double + -> SpVector Double + -- -> Int + -> BICGSTAB +bicgstab aa b x0 r0hat = + -- execState (replicateM n (modify (bicgstabStep aa r0hat))) bicgsInit where + execState (untilConverged _xBicgstab (bicgstabStep aa r0hat)) bicgsInit where + r0 = b ^-^ (aa #> x0) -- residual of initial guess solution + p0 = r0 + bicgsInit = BICGSTAB x0 r0 p0 + -- q (BICGSTAB xi _ _) = nor + +instance Show BICGSTAB where + show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" ++ + "r = " ++ show r ++ "\n" ++ + "p = " ++ show p ++ "\n" + + + + + + + + + + + +-- | ========= LINEAR SOLVERS INTERFACE + +data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show) + +-- random starting vector +linSolveM :: + PrimMonad m => + LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double) +linSolveM method aa b = do + let (m, n) = dim aa + nb = dim b + if n /= nb then error "linSolve : operand dimensions mismatch" else do + x0 <- randVec nb + case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0) + BICGSTAB_ -> return $ _x (cgs aa b x0 x0) + +-- deterministic starting vector (every component at 0.1) +linSolve :: + LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double +linSolve method aa b + | n /= nb = error "linSolve : operand dimensions mismatch" + | otherwise = solve aa b where + solve aa' b' | isDiagonalSM aa = (reciprocal aa') #> b' + | otherwise = solveWith aa' b' + solveWith aa' b' = case method of + CGS_ -> _xBicgstab (bicgstab aa' b' x0 x0) + BICGSTAB_ -> _x (cgs aa' b' x0 x0) + x0 = mkSpVectorD n $ replicate n 0.1 + (m, n) = dim aa + nb = dim b + +-- <\> : sets default solver method + +(<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double +(<\>) = linSolve BICGSTAB_ + + + + + + + + +-- | TODO : if system is poorly conditioned, is it better to warn the user or just switch solvers (e.g. via the pseudoinverse) ? + +-- linSolveQR aa b init f1 stepf +-- | isInfinite k = do +-- tell "linSolveQR : rank-deficient system" +-- | otherwise = do +-- solv aa b init +-- where +-- (q, r) = qr aa +-- k = conditionNumberSM r +-- solv aa b init = execState (untilConverged f1 stepf) init + + + + + + + + + + +-- | ========= PRETTY PRINTING + +-- | Show details and contents of sparse matrix + +sizeStr :: SpMatrix a -> String +sizeStr sm = + unwords ["(",show (nrows sm),"rows,",show (ncols sm),"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where + (SMInfo nz sy) = infoSM sm + + +showNonZero :: (Show a, Num a, Eq a) => a -> String +showNonZero x = if x == 0 then " " else show x + + +toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a] +toDenseRow (SM (_,ncol) im) irow = + fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1] + +toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String +toDenseRowClip sm irow ncomax + | ncols sm > ncomax = unwords (map show h) ++ " ... " ++ show t + | otherwise = show dr + where dr = toDenseRow sm irow + h = take (ncomax - 2) dr + t = last dr + +newline :: IO () +newline = putStrLn "" + +printDenseSM :: (Show t, Num t) => SpMatrix t -> IO () +printDenseSM sm = do + newline + putStrLn $ sizeStr sm + newline + printDenseSM' sm 5 5 + newline + where + printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO () + printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where + rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1] + rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_] + | otherwise = rr_ + + +toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String +toDenseListClip sv ncomax + | dim sv > ncomax = unwords (map show h) ++ " ... " ++ show t + | otherwise = show dr + where dr = toDenseListSV sv + h = take (ncomax - 2) dr + t = last dr + +printDenseSV :: (Show t, Num t) => SpVector t -> IO () +printDenseSV sv = do + newline + printDenseSV' sv 5 + newline where + printDenseSV' v nco = putStrLn rr_' where + rr_ = toDenseListClip v nco :: String + rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]] + | otherwise = rr_ + +class PrintDense a where + prd :: a -> IO () + +instance (Show a, Num a) => PrintDense (SpVector a) where + prd = printDenseSV + +instance (Show a, Num a) => PrintDense (SpMatrix a) where + prd = printDenseSM + + + + + + + + + + + + + +-- | rounding to 0 or 1 within some predefined numerical precision + + +almostZero, almostOne :: Double -> Bool +almostZero x = abs x <= eps +almostOne x = x >= (1-eps) && x < (1+eps) + +withDefault :: (t -> Bool) -> t -> t -> t +withDefault q d x | q x = d + | otherwise = x + +roundZero, roundOne :: Double -> Double +roundZero = withDefault almostZero 0 +roundOne = withDefault almostOne 1 + +with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t +with2Defaults q1 q2 d1 d2 x | q1 x = d1 + | q2 x = d2 + | otherwise = x + +roundZeroOne :: Double -> Double +roundZeroOne = with2Defaults almostZero almostOne 0 1 + + + + + +-- | transform state until a condition is met + +modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s +modifyUntil q f = do + x <- get + let y = f x + put y + if q y then return y + else modifyUntil q f + +loopUntilAcc :: Int -> ([t] -> Bool) -> (t -> t) -> t -> t +loopUntilAcc nitermax q f x = go 0 [] x where + go i ll xx | length ll < 2 = go (i + 1) (y : ll) y + | otherwise = if q ll || i == nitermax + then xx + else go (i + 1) (take 2 $ y:ll) y + where y = f xx + +-- modify state and append, until max # of iterations is reached +modifyInspectN :: + MonadState s m => Int -> ([s] -> Bool) -> (s -> s) -> m s +modifyInspectN nitermax q f + | nitermax > 0 = go 0 [] + | otherwise = error "modifyInspectN : n must be > 0" where + go i ll = do + x <- get + let y = f x + if length ll < 2 + then do put y + go (i + 1) (y : ll) + else if q ll || i == nitermax + then do put y + return y + else do put y + go (i + 1) (take 2 $ y : ll) + + +meanl :: (Foldable t, Fractional a) => t a -> a +meanl xx = 1/fromIntegral (length xx) * sum xx + +norm2l :: (Foldable t, Functor t, Floating a) => t a -> a +norm2l xx = sqrt $ sum (fmap (**2) xx) + + +diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1] + + + + + + + + + + +untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a +untilConverged fproj = modifyInspectN 100 (normDiffConverged fproj) + +-- convergence check (FIXME) +normDiffConverged :: (Foldable t, Functor t) => + (a -> SpVector Double) -> t a -> Bool +normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps + + + +-- run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0 + +runAppendN :: ([t] -> Bool) -> (t -> t) -> Int -> t -> [t] +runAppendN qq ff niter x0 | niter<0 = error "runAppendN : niter must be > 0" + | otherwise = go qq ff niter x0 [] where + go q f n z xs = + let x = f z in + if n <= 0 || q xs then xs + else go q f (n-1) x (x : xs) + +-- ", NO convergence check +runAppendN' :: (t -> t) -> Int -> t -> [t] +runAppendN' ff niter x0 | niter<0 = error "runAppendN : niter must be > 0" + | otherwise = go ff niter x0 [] where + go f n z xs = + let x = f z in + if n <= 0 then xs + else go f (n-1) x (x : xs) + +-- runN :: Int -> (a -> a) -> a -> a +-- runN n stepf x0 = runAppendN' stepf n x0 + + + + + + + + + + +-- | random matrices and vectors + +-- dense + +randMat :: PrimMonad m => Int -> m (SpMatrix Double) +randMat n = do + g <- MWC.create + aav <- replicateM (n^2) (MWC.normal 0 1 g) + let ii_ = [0 .. n-1] + (ix_,iy_) = unzip $ concatMap (zip ii_ . replicate n) ii_ + return $ fromListSM (n,n) $ zip3 ix_ iy_ aav + +randVec :: PrimMonad m => Int -> m (SpVector Double) +randVec n = do + g <- MWC.create + bv <- replicateM n (MWC.normal 0 1 g) + let ii_ = [0..n-1] + return $ fromListSV n $ zip ii_ bv + + + +-- sparse + +randSpMat :: Int -> Int -> IO (SpMatrix Double) +randSpMat n nsp | nsp > n^2 = error "randSpMat : nsp must be < n^2 " + | otherwise = do + g <- MWC.create + aav <- replicateM nsp (MWC.normal 0 1 g) + ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int) + jj <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int) + return $ fromListSM (n,n) $ zip3 ii jj aav + + +randSpVec :: Int -> Int -> IO (SpVector Double) +randSpVec n nsp | nsp > n = error "randSpVec : nsp must be < n" + | otherwise = do + g <- MWC.create + aav <- replicateM nsp (MWC.normal 0 1 g) + ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int) + return $ fromListSV n $ zip ii aav + + + + + + +-- | misc utils + + + +-- | integer-indexed ziplist +denseIxArray :: [b] -> [(Int, b)] +denseIxArray xs = zip [0..length xs-1] xs + +-- ", 2d arrays +denseIxArray2 :: Int -> [c] -> [(Int, Int, c)] +denseIxArray2 m xs = zip3 (concat $ replicate n ii_) jj_ xs where + ii_ = [0 .. m-1] + jj_ = concatMap (replicate m) [0 .. n-1] + ln = length xs + n = ln `div` m + + +-- folds + +foldrMap :: (Foldable t, Functor t) => (a -> b) -> (b -> c -> c) -> c -> t a -> c +foldrMap ff gg x0 = foldr gg x0 . fmap ff + +foldlStrict :: (a -> b -> a) -> a -> [b] -> a +foldlStrict f = go + where + go z [] = z + go z (x:xs) = let z' = f z x in z' `seq` go z' xs + +-- indexed fold? +ifoldr :: Num i => + (a -> b -> b) -> b -> (i -> c -> d -> a) -> c -> [d] -> b +ifoldr mjoin mneutral f = go 0 where + go i z (x:xs) = mjoin (f i z x) (go (i+1) z xs) + go _ _ [] = mneutral + + +-- bounds checking + +type LB = Int +type UB = Int + +inBounds :: LB -> UB -> Int -> Bool +inBounds ibl ibu i = i>= ibl && i (Int, Int) -> Bool +inBounds2 (ibl,ibu) (ix,iy) = inBounds ibl ibu ix && inBounds ibl ibu iy + + +-- ", lower bound = 0 +inBounds0 :: UB -> Int -> Bool +inBounds0 = inBounds 0 + +inBounds02 :: (UB, UB) -> (Int, Int) -> Bool +inBounds02 (bx,by) (i,j) = inBounds0 bx i && inBounds0 by j + + + + + + + + + +-- + +tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double +tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)] + +tv0, tv1 :: SpVector Double +tv0 = mkSpVectorD 2 [5, 6] + + +tv1 = SV 2 $ IM.singleton 0 1 + +-- wikipedia test matrix for Givens rotation + +tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3] + +tm1g1 = givens tm1 1 0 +tm1a2 = tm1g1 ## tm1 + +tm1g2 = givens tm1a2 2 1 +tm1a3 = tm1g2 ## tm1a2 + +tm1q = transposeSM (tm1g2 ## tm1g1) + + +-- wp test matrix for QR decomposition via Givens rotation + +tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41] + + + + +tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9] + +tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c] + where c= 0.4961 + s = 0.8682 + + +-- + +tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12] + + +-- playground + +-- | terminate after n iterations or when q becomes true, whichever comes first +untilC :: (a -> Bool) -> Int -> (a -> a) -> a -> a +untilC p n f = go n + where + go m x | p x || m <= 0 = x + | otherwise = True `seq` go (m-1) (f x) + + + + + +-- testing State + + +-- data T0 = T0 {unT :: Int} deriving Eq +-- instance Show T0 where +-- show (T0 x) = show x + +-- -- modifyT :: MonadState T0 m => (Int -> Int) -> m String +-- modifyT f = state (\(T0 i) -> (i, T0 (f i))) + + +-- t00 = T0 0 + +-- testT n = execState $ replicateM n (modifyT (+1)) + + +-- testT2 = execState $ when + + +-- replicateSwitch p m f = loop m where +-- loop n | n <= 0 || p = pure (#) +-- | otherwise = f *> loop (n-1) + + + +-- testing Writer + +-- asdfw n = runWriter $ do +-- tell $ "potato " ++ show n +-- tell "jam" +-- return (n+1) + + +-- -- + + +-- testing State and Writer + + + +-- runMyApp runA k maxDepth = +-- let config = maxDepth +-- state = 0 +-- in runStateT (runWriterT (runA k) config) state + + + diff --git a/src/Math/Linear/Sparse/IntMap.hs b/src/Math/Linear/Sparse/IntMap.hs new file mode 100644 --- /dev/null +++ b/src/Math/Linear/Sparse/IntMap.hs @@ -0,0 +1,149 @@ +module Math.Linear.Sparse.IntMap where + +import qualified Data.IntMap.Strict as IM + + + + +-- | ========= IntMap-of-IntMap (IM2) stuff + + +-- insert an element +insertIM2 :: + IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) +insertIM2 i j x imm = IM.insert i ro imm where + ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm) +{-# inline insertIM2 #-} + +-- lookup a key +lookupIM2 :: + IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a +lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j +{-# inline lookupIM2 #-} + +-- populate an IM2 from a list of (row index, column index, value) +fromListIM2 :: + Foldable t => + t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) +fromListIM2 iix sm = foldl ins sm iix where + ins t (i,j,x) = insertIM2 i j x t + + +-- | folding + +-- indexed fold over an IM2 +ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b +ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where + accRow acc i r = IM.foldlWithKey' (accElem i) acc r + accElem i acc j x = f i j x acc +{-# inline ifoldlIM2' #-} + +ifoldlIM2 :: + (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) -> + IM.IntMap (IM.IntMap t) -> + IM.IntMap a +ifoldlIM2 f m = IM.foldlWithKey' accRow IM.empty m where + accRow acc i row = IM.foldlWithKey' (accElem i) acc row + accElem i acc j x = f i j x acc +{-# inline ifoldlIM2 #-} + +foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b +foldlIM2 f empty mm = IM.foldl accRow empty mm where + accRow acc r = IM.foldl accElem acc r + accElem acc x = f x acc +{-# inline foldlIM2 #-} + + +-- transposeIM2 : inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed. +transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) +transposeIM2 = ifoldlIM2 (flip insertIM2) +{-# inline transposeIM2 #-} + +-- specialized folds + +-- -- extract diagonal elements +-- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a] +-- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) [] + + + + +-- | filtering + +-- map over outer IM and filter all inner IM's +ifilterIM2 :: + (IM.Key -> IM.Key -> a -> Bool) -> + IM.IntMap (IM.IntMap a) -> + IM.IntMap (IM.IntMap a) +ifilterIM2 f = + IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) +{-# inline ifilterIM2 #-} + +-- specialized filtering function + +-- keep only sub-diagonal elements +filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) +filterSubdiag = ifilterIM2 (\i j _ -> i>j) + +countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int +countSubdiagonalNZ im = + IM.size $ IM.filter (not . IM.null) (filterSubdiag im) + +-- list of (row, col) indices of (nonzero) subdiagonal elements +subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)] +subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where + im' = filterSubdiag im + +rpairs :: (a, [b]) -> [(a, b)] +rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj +rpairs (_, []) = [] + +-- -- list of (row, col) indices of elements that satisfy a criterion +-- indicesThatIM2 :: +-- (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)] +-- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where +-- im' = IM.filterWithKey f im + + + + +-- | mapping + +-- map over IM2 + +mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b) +mapIM2 = IM.map . IM.map -- imapIM2 (\_ _ x -> f x) + + +-- indexed map over IM2 +imapIM2 :: + (IM.Key -> IM.Key -> a -> b) -> + IM.IntMap (IM.IntMap a) -> + IM.IntMap (IM.IntMap b) +imapIM2 f im = IM.mapWithKey ff im where + ff j x = IM.mapWithKey (`f` j) x + + + +-- mapping keys + +mapKeysIM2 :: + (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) +mapKeysIM2 fi fj im = IM.map adjCols adjRows where + adjRows = IM.mapKeys fi im + adjCols = IM.mapKeys fj + + + + +-- map over a single `column` + +mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b) +mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im + + + + + +-- sparsification : + diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-}