diff --git a/Changelog.md b/Changelog.md new file mode 100644 --- /dev/null +++ b/Changelog.md @@ -0,0 +1,8 @@ +# Changelog + +## 0.0.1 + - raising upper bounds of optics and attoparsec + +## 0.0.0 +Initial version, providing PCA, DBScan and hierarchical clustering. +Command line interface to CREST XYZ trajectories. diff --git a/ConClusion.cabal b/ConClusion.cabal new file mode 100644 --- /dev/null +++ b/ConClusion.cabal @@ -0,0 +1,119 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.34.4. +-- +-- see: https://github.com/sol/hpack + +name: ConClusion +version: 0.0.1 +synopsis: Cluster algorithms, PCA, and chemical conformere analysis +description: Please see the README on GitLab at +category: Statistics, Chemistry +bug-reports: https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster/-/issues +author: Phillip Seeber +maintainer: phillip.seeber@googlemail.com +copyright: 2021 Phillip Seeber +license: AGPL-3 +license-file: LICENSE.md +build-type: Simple +extra-source-files: + README.md + LICENSE.md + Changelog.md + +source-repository head + type: git + location: https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster + +library + exposed-modules: + ConClusion.Chemistry.Topology + ConClusion.Numeric.Data + ConClusion.Numeric.Statistics + other-modules: + Paths_ConClusion + hs-source-dirs: + src + default-extensions: + BangPatterns + OverloadedStrings + NoImplicitPrelude + FlexibleContexts + ScopedTypeVariables + OverloadedLabels + DataKinds + FlexibleInstances + MultiParamTypeClasses + UndecidableInstances + TypeFamilies + DuplicateRecordFields + ScopedTypeVariables + DataKinds + DeriveAnyClass + DeriveDataTypeable + DeriveGeneric + DeriveTraversable + FlexibleInstances + GeneralizedNewtypeDeriving + TypeApplications + RecordWildCards + NamedFieldPuns + ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints + build-depends: + PSQueue >=1.1.0.1 && <1.2 + , aeson ==1.5.* + , attoparsec >=0.13.0.0 && <0.15 + , base >=4.7 && <4.15 + , containers >=0.6.0.0 && <0.7 + , formatting >=7.1.0 && <7.2 + , hmatrix >=0.20.0 && <0.21 + , massiv >=0.6.0.0 && <0.7 + , rio >=0.1.13.0 && <0.2 + default-language: Haskell2010 + +executable conclusion + main-is: ConClusion.hs + other-modules: + Paths_ConClusion + hs-source-dirs: + app + default-extensions: + BangPatterns + OverloadedStrings + NoImplicitPrelude + FlexibleContexts + ScopedTypeVariables + OverloadedLabels + DataKinds + FlexibleInstances + MultiParamTypeClasses + UndecidableInstances + TypeFamilies + DuplicateRecordFields + ScopedTypeVariables + DataKinds + DeriveAnyClass + DeriveDataTypeable + DeriveGeneric + DeriveTraversable + FlexibleInstances + GeneralizedNewtypeDeriving + TypeApplications + RecordWildCards + NamedFieldPuns + ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N + build-depends: + ConClusion + , PSQueue >=1.1.0.1 && <1.2 + , aeson ==1.5.* + , attoparsec >=0.13.0.0 && <0.15 + , base >=4.7 && <4.15 + , cmdargs >=0.10.0 && <0.11 + , containers >=0.6.0.0 && <0.7 + , formatting >=7.1.0 && <7.2 + , hmatrix >=0.20.0 && <0.21 + , massiv >=0.6.0.0 && <0.7 + , optics >=0.3 && <0.5 + , rio >=0.1.13.0 && <0.2 + , text >=1.2.0.0 && <1.3 + default-language: Haskell2010 diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,662 @@ + + GNU AFFERO GENERAL PUBLIC LICENSE + Version 3, 19 November 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 Affero General Public License is a free, copyleft license for +software and other kinds of works, specifically designed to ensure +cooperation with the community in the case of network server software. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +our General Public Licenses are 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. + + 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. + + Developers that use our General Public Licenses protect your rights +with two steps: (1) assert copyright on the software, and (2) offer +you this License which gives you legal permission to copy, distribute +and/or modify the software. + + A secondary benefit of defending all users' freedom is that +improvements made in alternate versions of the program, if they +receive widespread use, become available for other developers to +incorporate. Many developers of free software are heartened and +encouraged by the resulting cooperation. However, in the case of +software used on network servers, this result may fail to come about. +The GNU General Public License permits making a modified version and +letting the public access it on a server without ever releasing its +source code to the public. + + The GNU Affero General Public License is designed specifically to +ensure that, in such cases, the modified source code becomes available +to the community. It requires the operator of a network server to +provide the source code of the modified version running there to the +users of that server. Therefore, public use of a modified version, on +a publicly accessible server, gives the public access to the source +code of the modified version. + + An older license, called the Affero General Public License and +published by Affero, was designed to accomplish similar goals. This is +a different license, not a version of the Affero GPL, but Affero has +released a new version of the Affero GPL which permits relicensing under +this license. + + 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 Affero 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. Remote Network Interaction; Use with the GNU General Public License. + + Notwithstanding any other provision of this License, if you modify the +Program, your modified version must prominently offer all users +interacting with it remotely through a computer network (if your version +supports such interaction) an opportunity to receive the Corresponding +Source of your version by providing access to the Corresponding Source +from a network server at no charge, through some standard or customary +means of facilitating copying of software. This Corresponding Source +shall include the Corresponding Source for any work covered by version 3 +of the GNU General Public License that is incorporated pursuant to the +following paragraph. + + 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 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 work with which it is combined will remain governed by version +3 of the GNU General Public License. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU Affero 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 Affero 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 Affero 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 Affero 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 Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If your software can interact with users remotely through a computer +network, you should also make sure that it provides a way for users to +get its source. For example, if your program is a web application, its +interface could display a "Source" link that leads users to an archive +of the code. There are many ways you could offer source, and different +solutions will be better for different programs; see section 13 for the +specific requirements. + + 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 AGPL, see +. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,123 @@ +# ConClusion +ConClusion provides principal component analysis, hierarchical clustering and DBScan in Haskell. +There is also a command line interface for processing of [CREST](https://github.com/grimme-lab/crest) conformere trajectories. +Hence the name: CONformere CLUStering. +The procedure to analyse conformere data has three steps: + + 1. Read the trajectory and calculate a set of features for each conformere. + The features can include the energy, a set of bond lengths, a set of bond angles, and a set of dihedral angles in arbitrary combination. + Those descriptors form a feature matrix. + 2. A [principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) of the feature matrix might be perfomed to reduce the number of dimensions and remove redundancies. + 3. The (potentially PCA-processed) feature matrix is being clustered. + Different distance measures are available. + Either [DBScan](https://en.wikipedia.org/wiki/DBSCAN) or [hierarchical clustering](https://en.wikipedia.org/wiki/Hierarchical_clustering) can be used to group different conformeres. + +While the command line interface only fits the work flow described above, the underlying clustering algorithms and PCA are implemented in a general way and can be utilised independently as library. + +## Installation +### Bundled Archive +A self-contained executable archive is build for the main branch and for releases. +This can be executed directly on any Linux and has just to be downloaded. +Go to the [page of releaes](https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster/-/releases) and download an archive. +Make it executable (e.g. `chmod +x conclusion`) and you are done! + +### From Source +If you have the Haskell toolchain intalled and therefore working Cabal and GHC, you may build ConClusion from source. +This also requires working BLAS and LAPACK libraries on your system. + +```bash +git clone https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster.git ConClusion +cd ConClusion +cabal install --installdir=$(PREFIX) +``` +Choose a `PREFIX` where to install the executable. +`$HOME/.local/bin/` is often a good choice. + +If you would like to use ConClusion on systems where Nix is not available (Windows, BSD, ...) this is the way to go. + + +### With Nix +When you have Nix available on your system, everything can be build by Nix: + +```bash +git clone https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster.git ConClusion +cd ConClusion/nix +nix-build -A ConClusion.components.exes.conclusion +``` + + +## Usage +The command line interface to the `conclusion` executable offers full control about all three steps described above. + + - `-x --xyz` takes the path to the XYZ trajectory that is to be processed. + It must contain the energy as comment line (which is the case for CREST trajectories). + + - `--dim` specifies the characterising features; therefore a set of internal coordinates and maybe the energy. + A comma separated list of an arbitrary amount of features can be specified with the following syntax: + + - `e` for the energy + - `b m n` for a bond length, where `m` and `n` are the 0-based atom indices + - `a m n o` for an angle, where `m`, `n` and `o` are 0-based atom indices. + Calculates the angle around `n` + - `d m n o p` for a dihedral angle, where `m`, `n`, `o` and `p` are 0-based atom indices. + Calculates the rotation of the bond between `n` and `o`. + Dihedrals use a metric, that respects periodicity and direction of the rotation. + For each dihedral there will be two rows in the feature matrix, therefore. + See [this paper](https://doi.org/10.1002/prot.20310). + - (indices are 0-based) + + - `-p --pca` activates dimensionalty reduction by principal component analysis. + Give an integer to specify how many principal components are kept. + During the execution of ConClusion the error introduced by PCA will be printed. + + - `-c --cluster` activates clustering of the results. + The clustering algorithm can be selected by giving either `dbscan` or `hca` + + - `--measure` specifies the distance measure between the conformeres. + By default an euclidean distance is used, but Manhattan and Mahalanobis distances are also available, as well as a general L_p norm. + If Mahalanobis distances are used, it might be worth a try to disable PCA. + + - `--joinstrat` controls how inter-cluster distances are calculated in hierarchical clustering. + `single` might be the best choice to get dense groups of conformeres. + + - `--distance` gives the search radius in DBScan or the dendrogram cut distance in hierarchical clustering. + + - `--minsize` is the minimum size of a cluster in DBScan. + If `--forcemin` is given, the clusters obtained by HCA are also filtered for a minimum size. + + - `--forcemin` forces filtering of HCA clusters for their minimum size as given by `--minsize`. + Disabled by default. + +Each processing step will produce a Gnuplot compatible file (space separated columns). The pure feature matrix will be `features.dat`, the results of the PCA will be in `pca.dat` and the clustering results will be in `cluster.dat`. The first column in `cluster.dat` will be an integer giving the cluster number this point belongs to, that can be used for colour-coding in Gnuplot. + +## Example +A perylene dye with four phenoxy groups has different conformeres, that have different spectral properties. +For solubility the dye has also some alkyl groups. +Crest finds about 1400 conformeres, most of them being different only in the alkyl side-chains, that do not influence spectral properties. +Therefore, a much smaller group of different conformeres with respect to different positions of phenoxy groups exist. +From each of those groups the lowest energy conformere shall be obtained. +We therefore select eight dihedral angles and the energy as features; 2 dihedrals for each phenoxy group. +One dihedral per phenoxy group describing the rotation of the perylene-O bond, the second one describing the rotation around the O-Ph bond. +As the dihedral angles are not independent from each other, as some orientations of phenoxy groups are not possible, we use a PCA to reduce dimensionalty and remove redundancies. +After the PCA, DBScan is used to obtain clusters of similar conformeres. +The lowest index in each cluster is also the lowest energy conformere in each group, as CREST sorts conformeres by energy. + +```bash +conclusion \ + --xyz=crest_conformers.xyz \ + --pca=3 \ + --dim="e, d 19 18 2 56, d 18 2 56 65, d 16 15 1 45, d 15 1 45 46, d 11 13 0 34, d 13 0 34 43, d 29 31 3 67, d31 3 67 76" \ + --measure=manhattan \ + --cluster=dbscan \ + --distance=0.3 \ + --minsize=5 +``` + +## Library/Haskell Package +ConClusion provides principal components analysis and the clustering algorithms DBScan and hierarchical clustering. +The algorithms are implemented in efficient parallel arrays and perform quite well. +For the API see the haddock documentation, which can be generated by: + +```bash +cabal haddock +``` diff --git a/app/ConClusion.hs b/app/ConClusion.hs new file mode 100644 --- /dev/null +++ b/app/ConClusion.hs @@ -0,0 +1,742 @@ +-- | +-- Module : Main +-- Description : Command Line Interface to analyse CREST results. +-- Copyright : Phillip Seeber, 2021 +-- License : AGPL-3 +-- Maintainer : phillip.seeber@uni-jena.de +-- Stability : experimental +-- Portability : POSIX, Windows +module Main (main) where + +import ConClusion.Chemistry.Topology hiding (xyz) +import qualified ConClusion.Numeric.Statistics as Statistics +import Data.Aeson +import Data.Attoparsec.Text hiding (D) +import Data.Foldable +import qualified Data.IntSet as IntSet +import Data.Massiv.Array as Massiv hiding (B, D) +import qualified Data.Massiv.Array as Massiv +import Data.Text.Lazy (toStrict) +import Data.Text.Lazy.Builder as TB +import Data.Version (showVersion) +import Formatting hiding (char, (%)) +import qualified Formatting as F +import Optics hiding (view) +import Paths_ConClusion (version) +import RIO hiding (Lens', Vector, lens, view, (^.)) +import qualified RIO.Text as Text +import System.Console.CmdArgs +import System.IO.Unsafe (unsafePerformIO) + +main :: IO () +main = do + app <- cmdArgs conClusionCmd >>= cmdToEnv + runRIO app $ do + logInfo $ "ConClusion, version " <> (displayShow . showVersion $ version) <> "\n" + + -- Processing. + (featureMat, _trj) <- processTrajectory + pcaMat <- dimReduction featureMat + _clusters <- clusterFeatures pcaMat + + logInfo "Done. Good Luck with the results!" + +---------------------------------------------------------------------------------------------------- +-- Processing Logic + +-- | Read the trajectory and calculate the specified features. +{-# SCC processTrajectory #-} +processTrajectory :: (HasXYZ env, HasDim env, HasLogFunc env) => RIO env (Matrix DL Double, Trajectory) +processTrajectory = do + -- Obtain environment + file <- view xyzL + dims <- view dimL + let outputFile = "features.dat" + + -- Logging + {- ORMOLU_DISABLE -} + logInfo $ + "+-----------------------+\n\ + \| Trajectory Processing |\n\ + \+-----------------------+\n\n\ + \ Trajectory File : " <> displayShow file <> "\n\ + \ Feature Matrix File : " <> displayShow outputFile <> "\n\ + \ Features :\n" <> displayFeatures dims + {- ORMOLU_ENABLE -} + + -- Processing. + raw <- readFileUtf8 file + trj <- handleFailure $ parse' trajectory raw + features <- handleFailure $ getFeatures dims trj + + -- Output. + toGnuPlot Nothing (Just . Massiv.fromList Seq $ dims) (compute features) + >>= writeFileUtf8 outputFile + + logInfo "\n\n\n" + return (features, trj) + +-- | Possibly perform a PCA. +{-# SCC dimReduction #-} +dimReduction :: (HasPrincipalComponentAnalysis env, HasLogFunc env) => Matrix DL Double -> RIO env (Matrix U Double) +dimReduction mat = do + doPCA <- view pcaL + case doPCA of + Nothing -> return . compute $ featureMat + Just PrincipalComponentAnalysis {keep} -> do + let outputFile = "pca.dat" + + -- Logging. + {- ORMOLU_DISABLE -} + logInfo + "+------------------------------+\n\ + \| Principal Component Analysis |\n\ + \+------------------------------+\n" + {- ORMOLU_ENABLE -} + + -- Processing + pcaData <- Statistics.pca keep featureMat + let behaviourPercent = Statistics.remaining pcaData + + -- Logging. + {- ORMOLU_DISABLE -} + logInfo $ + " Dimensions : " <> display keep <> "\n\ + \ Eigenvalues : " <> (displayEigenvalues . Statistics.allEigenValues $ pcaData) <> "\n\ + \ Mean Squared Error : " <> (displayMSE . Statistics.mse $ pcaData) <> "\n\ + \ Behaviour Captured : " <> displayPercent behaviourPercent <> "\n\ + \ PCA Matrix File : " <> displayShow outputFile + {- ORMOLU_ENABLE -} + + -- Output + toGnuPlot Nothing Nothing (Statistics.y pcaData) >>= writeFileUtf8 outputFile + + logInfo "\n\n\n" + return . Statistics.y $ pcaData + where + featureMat = compute @U mat + +-- | Cluster analysis of the (dimensionalty reduced) feature matrix. +{-# SCC clusterFeatures #-} +clusterFeatures :: + (HasClustering env, HasLogFunc env) => Matrix U Double -> RIO env Statistics.Clusters +clusterFeatures mat = do + doClustering <- view clusteringL + case doClustering of + Nothing -> let Sz (_ :. n) = size mat in return $ makeArray @Massiv.B @Ix1 @IntSet Par (Sz n) $ \i -> IntSet.singleton i + Just cl -> do + let outputFile = "cluster.dat" + -- Logging + {- ORMOLU_DISABLE -} + logInfo + "+------------+\n\ + \| Clustering |\n\ + \+------------+\n" + {- ORMOLU_ENABLE -} + + clusters <- case cl of + DBScan {distance, measure, minSize} -> clusterDB distance measure minSize + HCA {distance, measure, joinstrat, forcemin} -> clusterHCA distance measure joinstrat forcemin + + {- ORMOLU_DISABLE -} + logInfo $ + " Number of Clusters : " <> (display . Massiv.elemsCount $ clusters) <> "\n\ + \ Cluster Data File : " <> displayShow outputFile <> "\n\ + \ Cluster : " <> displayClusters clusters + {- ORMOLU_ENABLE -} + + -- Output + toGnuPlot (Just clusters) Nothing mat >>= writeFileUtf8 outputFile + + return clusters + where + -- Make the distance function for clustering. + distFn m = case m of + Lr i -> Statistics.lpNorm i + Manhattan -> Statistics.manhattan + Euclidean -> Statistics.euclidean + Mahalanobis -> Statistics.mahalanobis + + -- Cluster data with DBScan algorithm. + clusterDB d m s = do + clusters <- Statistics.dbscan (distFn m) s d mat + {- ORMOLU_DISABLE -} + logInfo $ + " Algorithm : DBScan\n\ + \ Distance Measure : " <> displayShow m <> "\n\ + \ Search Distance ε : " <> display d <> "\n\ + \ Minimal Cluster Size : " <> display s + {- ORMOLU_ENABLE -} + return clusters + + -- Cluster data with Hierarchical Cluster Analysis. + clusterHCA d m j fm = do + dendrogram <- Statistics.hca (distFn m) j mat + let allClusters = Statistics.cutDendroAt dendrogram d + dendroFile = "dendrogram.json" + clusters = case fm of + Nothing -> allClusters + Just minSize -> compute . sfilter (\c -> IntSet.size c >= minSize) $ allClusters + + -- Write dendrogram file to disk. + liftIO $ encodeFile dendroFile dendrogram + + {- ORMOLU_DISABLE -} + logInfo $ + " Algorithm : Hierarchical Cluster Analysis\n\ + \ Distance Measure : " <> displayShow m <> "\n\ + \ Search Distance : " <> display d <> "\n\ + \ Cluster Join Strategy : " <> displayShow j <> "\n" <> + (case fm of + Nothing -> mempty + Just minSize -> " Minimal Cluster Size : " <> display minSize <> "\n") <> + " Dendrogram File : " <> displayShow dendroFile + {- ORMOLU_ENABLE -} + return clusters + +---------------------------------------------------------------------------------------------------- +-- Command line interfaces + +-- | Command line arguments for ConClusion. +data ConClusionCmd = ConClusionCmd + { -- | Filepath to the CREST xyz trajectory file. + xyz :: FilePath, + -- | Dimensions to process. + dim :: String, + -- | If PCA should be used to reduce dimensionalty of the problem. If so, the number of + -- dimensions to keep. + pca :: Maybe Int, + -- | Cluster algorithm to be used for clustering, if any. + cluster :: Maybe String, + -- | Distance measure of points. + measure :: String, + -- | Cluster joining strategy for hierarchical clustering. + joinstrat :: String, + -- | Search distance in case of DBScan or distance for the Dendogram cut in case of hierarchical + -- clustering. + distance :: Double, + -- | Mininmal size of a cluster in case of DBScan. + minsize :: Int, + -- | If to force a minimum size of clusters in HCA. + forcemin :: Bool, + -- | Do verbose logging. + verbose :: Bool + } + deriving (Show, Data, Typeable) + +conClusionCmd :: ConClusionCmd +conClusionCmd = + ConClusionCmd + { xyz = + "crest_conformers.xyz" + &= typFile + &= help "Path to CREST trajectory.", + dim = + "e" + &= help + "Expression of dimensions to analyse:\n\ + \ - \"e\" Energy\n\ + \ - \"b m n\" Bond length between the atoms with\n\ + \ indices m and n\n\ + \ - \"a m n o\" Angle between the atoms with indices\n\ + \ m, n and o\n\ + \ - \"d m n o p\" Sinus of the dihedral between the atoms\n\ + \ m, n, o and p\n\ + \Multiple features may be combined by giving a comma separated list of them.\n\ + \Default: e", + pca = + Nothing + &= typ "INT" + &= help + "Perform a principal component analysis on the features.\ + \Optionally give the number of dimensions you want to keep.\n\ + \Default: 2" + &= opt (2 :: Int), + cluster = + Nothing + &= help + "Cluster algorithm to cluster conformeres.\n\ + \ - \"dbscan\" DBScan algorithm. Depends on \n\ + \ \"measure\", \"distance\" and \"minSize\"\n\ + \ - \"hca\" Hierarchical Cluster Analysis. Depends on\n\ + \ \"measure\", \"joinstrat\" and \"distance\"\n\ + \Default: dbscan" + &= opt ("dbscan" :: String), + measure = + "mahalanobis" + &= help + "Distance measure for distance between conformeres.\n\ + \ - \"lrX\" General L_r norm, where X is a value\n\ + \ between 3 and 9\n\ + \ - \"manhattan\" Manhattan distances\n\ + \ - \"euclidean\" Euclidean distances\n\ + \ - \"mahalanobis\" Mahalanobis distances\n\ + \Default: mahalanobis", + joinstrat = + "ward" + &= help + "Strategy for inter-cluster distance in hierarchical clustering.\n\ + \ - \"single\" Single-Linkage (minimal distance of\n\ + \ points between clusters)\n\ + \ - \"complete\" Complete-Linkage (maximal distance of\n\ + \ points between clusters)\n\ + \ - \"median\" Median distance linkage between\n\ + \ all points between clusters)\n\ + \ - \"upgma\" Average-Group-Linkage (average distance\n\ + \ of all points of joined clusters)\n\ + \ - \"wpgma\" Average-Linkage (average distance\n\ + \ of all points between clusters)\n\ + \ - \"centroid\" Distance of cluster centroids\n\ + \ - \"ward\" Ward distance of minimal residues\n\ + \Default: ward", + distance = + 0.1 + &= typ "FLOAT" + &= help + "Search distance in DBScan or Dendogram cut distance in hierarchical clustering.\n\ + \Default: 0.1", + minsize = + 5 + &= typ "INT" + &= help + "Minimal size of a cluster in DBScan.\n\ + \Default: 5", + forcemin = + False + &= help + "Force a minimal cluster size in hierarchical clustering.\n\ + \Default: False", + verbose = + False + &= help "Verbose logging." + } + &= summary + "Analysis of conformere ensembles.\n\ + \Processes an annotated XYZ trajectory (as generated by CREST, for example), where the energy\ + \is provided in the comment line. The analysis is a three-step procedure:\n\ + \ 1. Calculate selected features from the ensemble, that are meant to\n\ + \ characterise the conformeres.\n\ + \ Relevant settings: \"--dim\"\n\ + \ 2. (Optional) Perform a PCA analysis on the feature matrix to reduce the\n\ + \ number of dimensions.\n\ + \ Relevant settings: \"--pca\"\n\ + \ 3. (Optional) Cluster the feature matrix into groups of conformeres. If a PCA\n\ + \ has been performed, use the feature matrix with reduced dimensionalty.\n\ + \ Relevant settins: \"--cluster\", \"--measure\", \"joinstrat\", \"--distance\",\n\ + \ \"--minsize\"" + +---------------------------------------------------------------------------------------------------- +-- Runtime Environment. + +-- Reader classes +class HasXYZ a where + xyzL :: Lens' a FilePath + +class HasDim a where + dimL :: Lens' a [Feature] + +class HasPrincipalComponentAnalysis a where + pcaL :: Lens' a (Maybe PrincipalComponentAnalysis) + +class HasClustering a where + clusteringL :: Lens' a (Maybe Clustering) + +-- | Full runtime environment of ConClusion with all settings combined. +data App = App + { -- | Filepath to the CREST xyz trajectory file. + xyz :: FilePath, + -- | Features that are to be analysed. + dim :: [Feature], + -- | Settings for principal component analysis. + pca :: Maybe PrincipalComponentAnalysis, + -- | Settings for Clustering. + clustering :: Maybe Clustering, + -- | RIO's 'logFunc'. + logFunc :: LogFunc + } + +-- Lenses +instance (k ~ A_Lens, a ~ FilePath, b ~ a) => LabelOptic "xyz" k App App a b where + labelOptic = lens (xyz :: App -> FilePath) $ \s b -> (s {xyz = b} :: App) + +instance (k ~ A_Lens, a ~ [Feature], b ~ a) => LabelOptic "dim" k App App a b where + labelOptic = lens (dim :: App -> [Feature]) $ \s b -> (s {dim = b} :: App) + +instance (k ~ A_Lens, a ~ Maybe PrincipalComponentAnalysis, b ~ a) => LabelOptic "pca" k App App a b where + labelOptic = lens (pca :: App -> Maybe PrincipalComponentAnalysis) $ \s b -> (s {pca = b} :: App) + +instance (k ~ A_Lens, a ~ Maybe Clustering, b ~ a) => LabelOptic "clustering" k App App a b where + labelOptic = lens clustering $ \s b -> s {clustering = b} + +instance (k ~ A_Lens, a ~ LogFunc, b ~ a) => LabelOptic "logFunc" k App App a b where + labelOptic = lens logFunc $ \s b -> s {logFunc = b} + +-- Reader Classes +instance HasXYZ App where + xyzL = castOptic #xyz + +instance HasDim App where + dimL = castOptic #dim + +instance HasPrincipalComponentAnalysis App where + pcaL = castOptic #pca + +instance HasClustering App where + clusteringL = castOptic #clustering + +instance HasLogFunc App where + logFuncL = toLensVL #logFunc + +-- | Settings for principal component analysis. +newtype PrincipalComponentAnalysis = PrincipalComponentAnalysis + { -- | Number of dimensions to keep from PCA. + keep :: Int + } + +-- Lenses +instance (k ~ A_Lens, a ~ Int, b ~ a) => LabelOptic "keep" k PrincipalComponentAnalysis PrincipalComponentAnalysis a b where + labelOptic = lens keep $ \s b -> s {keep = b} + +-- | Settings for clustering. +data Clustering + = -- | DBScan. + DBScan + { -- | Search distance around a point. + distance :: Double, + -- | Distance measure in clustering. + measure :: Measure, + -- | Minimal size of a cluster. + minSize :: Int + } + | -- | Hierarchical clustering. + HCA + { -- | Cut distance in the dendogram to select the number of clusters. + distance :: Double, + -- | Distance measure in clustering. + measure :: Measure, + -- | Join strategy for clusters. + joinstrat :: Statistics.JoinStrat Double, + -- | A minimum cluster size might be forced in post processing for HCA. + forcemin :: Maybe Int + } + +-- | Distance measure between clusters. +data Measure + = Lr Int + | Manhattan + | Euclidean + | Mahalanobis + deriving (Show) + +---------------------------------------------------------------------------------------------------- +-- Utilities + +-- | Conversion of command line arguments into runtime environment. +cmdToEnv :: ConClusionCmd -> IO App +cmdToEnv ConClusionCmd {xyz, dim, pca, cluster, measure, joinstrat, distance, minsize, forcemin, verbose} = do + -- Dimensionalty/Feature construction. + dimApp <- case parseOnly dimParser (Text.pack dim) of + Left err -> throwM $ ArgumentException "dim" err "a comma separated list of dimension specifications" + Right res -> return res + + -- PCA construction. + pcaSettings <- case pca of + Nothing -> return Nothing + Just x -> + if x >= 1 + then return . Just $ PrincipalComponentAnalysis {keep = x} + else throwM $ ArgumentException "pca" (show x) "a positive integer smaller than the number of clusters." + + -- Clustering construction. + clDistance <- + if distance <= 0 + then throwM $ ArgumentException "distance" (show distance) "a positive number" + else return distance + clMeasure <- case measure of + "lr3" -> return $ Lr 3 + "lr4" -> return $ Lr 4 + "lr5" -> return $ Lr 5 + "lr6" -> return $ Lr 6 + "lr7" -> return $ Lr 7 + "lr8" -> return $ Lr 8 + "lr9" -> return $ Lr 9 + "manhattan" -> return Manhattan + "euclidean" -> return Euclidean + "mahalanobis" -> return Mahalanobis + _ -> throwM $ ArgumentException "measure" measure "one of \"manhattan | euclidean | mahalanobis | lrX\" with X = [ 3 .. 9]" + clMinsize <- + if minsize <= 0 + then throwM $ ArgumentException "minsize" (show minsize) ">= 1" + else return minsize + clJoinstrat <- case joinstrat of + "single" -> return Statistics.SingleLinkage + "complete" -> return Statistics.CompleteLinkage + "median" -> return Statistics.Median + "upgma" -> return Statistics.UPGMA + "wpgma" -> return Statistics.WPGMA + "centroid" -> return Statistics.Centroid + "ward" -> return Statistics.Ward + _ -> throwM $ ArgumentException "joinstrat" joinstrat "one of \"single | complete | upgma | wpgma | centroid | ward\"" + clustering <- case cluster of + Just "dbscan" -> + return . Just $ + DBScan + { distance = clDistance, + measure = clMeasure, + minSize = clMinsize + } + Just "hca" -> + return . Just $ + HCA + { distance = clDistance, + measure = clMeasure, + joinstrat = clJoinstrat, + forcemin = if forcemin then Just clMinsize else Nothing + } + Nothing -> return Nothing + Just c -> throwM $ ArgumentException "cluster" c "one of \"dbscan | hca\"" + + -- Log function setup. + logOptions <- + id + . setLogUseTime False + . setLogUseColor True + . setLogVerboseFormat False + . setLogUseLoc False + <$> logOptionsHandle stdout verbose + (lf, _ :: IO ()) <- newLogFunc logOptions + return + App + { xyz = xyz, + dim = dimApp, + pca = pcaSettings, + clustering = clustering, + logFunc = lf + } + +-- | Parser for dimensionalty expressions. +dimParser :: Parser [Feature] +dimParser = many1 $ do + skipMany (char ' ') + fType <- eParser <|> bParser <|> aParser <|> dParser + skipMany (char ' ') + _ <- option ',' (char ',') + return fType + where + eParser = char 'e' *> return Energy + bParser = do + _ <- char 'b' + skipMany (char ' ') + a <- decimal + skipMany (char ' ') + b <- decimal + return . Bond $ B a b + aParser = do + _ <- char 'a' + skipMany (char ' ') + a <- decimal + skipMany (char ' ') + b <- decimal + skipMany (char ' ') + c <- decimal + return . Angle $ A a b c + dParser = do + _ <- char 'd' + skipMany (char ' ') + a <- decimal + skipMany (char ' ') + b <- decimal + skipMany (char ' ') + c <- decimal + skipMany (char ' ') + d <- decimal + return . Dihedral $ D a b c d + +-- | Substitution for RIO's @view@ for optics. +view :: (Is k A_Getter, MonadReader s f) => Optic' k is s b -> f b +view l = (^. l) <$> ask + +-- | Parsers lifted to 'MonadThrow'. +parse' :: MonadThrow m => Parser x -> Text -> m x +parse' p t = case parseOnly p t of + Left err -> throwM $ ParseException err + Right res -> return res + +-- | Logging a problem as an error and then rethrow the underlying exception. +handleFailure :: (Exception e, HasLogFunc env) => Either e r -> RIO env r +handleFailure (Right r) = return r +handleFailure (Left e) = do + logError . displayShow $ e + throwM e + +-- | Display features for a trajectory in a nice way. +displayFeatures :: Foldable f => f Feature -> Utf8Builder +displayFeatures features = + foldl + ( \acc f -> + acc + <> ( case f of + Energy -> " - E" + Bond (B a b) -> " - B " <> display a <> " " <> display b + Angle (A a b c) -> " - A " <> display a <> " " <> display b <> " " <> display c + Dihedral (D a b c d) -> " - D " <> display a <> " " <> display b <> " " <> display c <> " " <> display d + ) + <> "\n" + ) + mempty + features + +-- | Display a percent value. +displayPercent :: Double -> Utf8Builder +displayPercent v = display $ sformat (fixed 1 F.% " %") v + +-- | Display a mean squared error. +displayMSE :: Double -> Utf8Builder +displayMSE v = display $ sformat (fixed 2) v + +-- | Display eigenvalues +displayEigenvalues :: Vector U Double -> Utf8Builder +displayEigenvalues vec = display . Massiv.fold . Massiv.map fmtD $ vec + where + fmtD d = format (right 6 ' ' F.%. fixed 2) d + +-- | Display best cluster members. +displayBestMember :: Vector U Int -> Utf8Builder +displayBestMember vec = + display . toLazyText . Massiv.fold . Massiv.map (\i -> bformat (" - " F.% int F.% "\n") i) $ vec + +-- | Converts a matrix to a Gnuplot matrix. Prints a header line as a comment with the features, if +-- available. Adds the assigned cluster number in the first column if available. +{-# SCC toGnuPlot #-} +toGnuPlot :: + ( MonadThrow m + ) => + Maybe Statistics.Clusters -> + Maybe (Vector Massiv.B Feature) -> + Matrix U Double -> + m Text +toGnuPlot clusters features mat = do + -- Build the lines representing the points. If cluster information is available, the first column + -- is the cluster, that this point is assigned to. + valueLines <- case clusters of + Nothing -> return . Massiv.fold . Massiv.map (writeLine . Left) . innerSlices $ mat + Just cl -> do + let points = innerSlices mat :: Vector Massiv.D (Vector M Double) + labeledCl = compute . Massiv.imap (,) $ cl :: Vector Massiv.B (Int, IntSet) + + -- Anotate all points with their assignment to a cluster number. If a point is from DBScan, it + -- is not necessarily assigned to a cluster. + annoPoints <- Massiv.iforM @Massiv.B points $ \pointIx vec -> do + let filteredCl = Massiv.dropWhile (\(_, clPoints) -> not $ pointIx `IntSet.member` clPoints) labeledCl + let assignedCl = filteredCl !? 0 :: Maybe (Int, IntSet) + clNumber = fst <$> assignedCl :: Maybe Int + return (clNumber, vec) + + -- return annoClusters + return . Massiv.fold . Massiv.map (writeLine . Right) $ annoPoints + + featureLine <- case features of + Nothing -> return mempty + Just f -> do + let Sz nFeatureLabels = actualFeatureSize f + featureLine = writeFeatures f + if nFeatures == nFeatureLabels + then return featureLine + else throwM $ SizeMismatchException (Sz nFeatureLabels) (Sz nFeatures) + return . toStrict . TB.toLazyText $ featureLine <> valueLines + where + Sz (nFeatures :. _) = size mat + + -- The column width in the Gnuplot file. + cw = 14 + + -- Formatter for doubles. + dForm a = bformat (left cw ' ' F.%. fixed 4) a + + -- Formatter for Integers. + iForm = fitLeft 3 F.%. left 3 ' ' F.%. int + + -- Formatter for Features + fForm f = case f of + Energy -> bformat (left cw ' ' F.%. builder) "E" + Bond (B a b) -> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm)) "B" a b + Angle (A a b c) -> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A" a b c + Dihedral (D a b c d) -> + bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A sin" a b c d + <> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A cos" a b c d + + -- Print a vector of numeric values. + doubleVecB :: Source r Ix1 Double => Vector r Double -> TB.Builder + doubleVecB v = Massiv.fold . Massiv.map dForm $ v + + -- Writer for a value line. + writeLine :: + (Source r Ix1 Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder + writeLine (Left vec) = doubleVecB vec <> "\n" + writeLine (Right (cl, vec)) = bformat iForm (fromMaybe (-1) cl) <> doubleVecB vec <> "\n" + + -- Writer for the Feature line. + writeFeatures :: Vector Massiv.B Feature -> TB.Builder + writeFeatures vec = "#" <> (Massiv.fold . Massiv.map fForm $ vec) <> "\n" + + -- Actual feature size. + actualFeatureSize vec = + unsafePerformIO $ + Massiv.foldlP + ( \acc f -> case f of + Dihedral _ -> 2 + acc + _ -> 1 + acc + ) + 0 + (+) + 0 + vec + +-- | Pretty printer for clusters. +displayClusters :: Statistics.Clusters -> Utf8Builder +displayClusters clusters = display . toLazyText . Massiv.fold . Massiv.map printCl $ clusters + where + -- Print clusters as list of indices. + iForm = fitLeft 5 F.%. left 5 ' ' F.%. int + + -- Printer of a single cluster. + printCl :: IntSet -> TB.Builder + printCl cl = + fst $ + IntSet.foldl + ( \(bAcc, counter) i -> + let this = bformat iForm i <> if counter `mod` (16 :: Int) == 0 then "\n " else mempty + in (bAcc <> this, counter + 1) + ) + ("\n\n - ", 1) + cl + +---------------------------------------------------------------------------------------------------- +-- Exceptions + +-- | Exception thrown on unexcpected command line arguments. +data ArgumentException = ArgumentException + { argName :: String, + is :: String, + should :: String + } + +instance Show ArgumentException where + show ArgumentException {argName, is, should} = + "ArgumentException: Field \"" + <> argName + <> "\" got \"" + <> is + <> "\" but must be " + <> should + <> "." + +instance Exception ArgumentException + +-- | Exceptions thrown on failed parsing. +newtype ParseException = ParseException String + +instance Show ParseException where + show (ParseException err) = "ParseException: \"" <> err <> "\"" + +instance Exception ParseException diff --git a/src/ConClusion/Chemistry/Topology.hs b/src/ConClusion/Chemistry/Topology.hs new file mode 100644 --- /dev/null +++ b/src/ConClusion/Chemistry/Topology.hs @@ -0,0 +1,210 @@ +-- | +-- Module : ConClusion.Chemistry.Topology +-- Description : Principal Component Analysis +-- Copyright : Phillip Seeber, 2021 +-- License : AGPL-3 +-- Maintainer : phillip.seeber@googlemail.com +-- Stability : experimental +-- Portability : POSIX, Windows +-- +-- This module implements routines to work with simple molden style XYZ trajectories. This includes +-- parsers as well as functions to obtain structural features in internal coordinates. +-- +-- For an introduction into PCA see . +-- +-- Diherdrals require a special metric, see . +module ConClusion.Chemistry.Topology + ( Molecule, + Trajectory, + xyz, + trajectory, + B (..), + A (..), + D (..), + Feature (..), + getFeatures, + ) +where + +import ConClusion.Numeric.Data hiding (angle) +import qualified ConClusion.Numeric.Data as ND +import Data.Attoparsec.Text hiding (D) +import Data.Foldable +import Data.Massiv.Array as Massiv hiding (B, D, forM) +import qualified Numeric.LinearAlgebra as LA +import RIO +import qualified RIO.Seq as Seq +import qualified RIO.Text as Text +import qualified RIO.Vector.Boxed as VectorB + +-- | A Molecule in cartesian coordinates. +data Molecule = Molecule + { -- | The energy of the molecule. + energy :: !Double, + -- | Chemical symbols or names of the atoms. \(N\) vector. + atoms :: !(VectorB.Vector Text), + -- | Cartesian coordinates. Atoms as rows, xyz as columns. \(N \times 3\) matrix. + coordinates :: !(Massiv.Matrix S Double) + } + +type Trajectory = Seq Molecule + +-- | Parser for molecules in Molden XYZ format. +xyz :: Parser Molecule +xyz = do + -- Starting line with number of atoms. + nAtoms <- do + _ <- skipMany . char $ ' ' + n <- decimal + _ <- skipMany . char $ ' ' + endOfLine + return n + + -- Comment line. Must contain the energy. + energy <- do + _ <- skipMany . char $ ' ' + e <- double + _ <- skipMany . char $ ' ' + endOfLine + return e + + -- Multiple lines of atoms. + atoms <- count nAtoms $ do + _ <- skipMany . char $ ' ' + name <- many1 $ letter <|> digit + _ <- skipMany1 . char $ ' ' + x <- double + _ <- skipMany1 . char $ ' ' + y <- double + _ <- skipMany1 . char $ ' ' + z <- double + _ <- skipMany . char $ ' ' + endOfLine + let row = Massiv.expandOuter (Sz 1) const . Massiv.fromList @S Seq $ [x, y, z] + return (Text.pack name, row) + + return + Molecule + { energy = energy, + atoms = VectorB.fromListN nAtoms . fmap fst $ atoms, + coordinates = Massiv.compute . Massiv.concat' (Dim 2) . fmap snd $ atoms + } + +-- | Parser for trajectories in XYZ format as produced by CREST. +{-# SCC trajectory #-} +trajectory :: Parser Trajectory +trajectory = do + mols <- Seq.fromList <$> many1 xyz + skipSpace + endOfInput + return mols + +-- | Selection of a bond between two atoms. +data B = B Int Int + +-- | Selection of an angle between three atoms. +data A = A Int Int Int + +-- | Selection of a dihedral angle between four atoms. Rotation around the central two. +data D = D Int Int Int Int + +-- | Selections +data Feature + = Energy + | Bond B + | Angle A + | Dihedral D + +-- | Calculate a bond distance. +bond :: MonadThrow m => B -> Molecule -> m Double +bond (B a b) Molecule {coordinates} + | a == b = throwM . IndexException $ "selected atoms are identical" + | otherwise = do + coordA <- compute @U <$> (coordinates !?> a) + coordB <- compute @U <$> (coordinates !?> b) + vecAB <- coordA .-. coordB + return . sqrt . Massiv.sum . Massiv.map (^ (2 :: Int)) $ vecAB + +-- | Calculates the sinus of an angle defined by three atoms. +angle :: MonadThrow m => A -> Molecule -> m Double +angle (A a b c) Molecule {coordinates} + | a == b || b == c || a == c = throwM . IndexException $ "selected atoms are identical" + | otherwise = do + coordA <- compute @U <$> (coordinates !?> a) + coordB <- compute @U <$> (coordinates !?> b) + coordC <- compute @U <$> (coordinates !?> c) + vecAB <- coordA .-. coordB + vecCB <- coordC .-. coordB + return $ ND.angle vecAB vecCB + +-- | Calculates the dihedral angle defined by four atoms. Respects rotation direction. Obtains the +-- result in radian. +dihedral' :: MonadThrow m => D -> Molecule -> m Double +dihedral' (D a b c d) Molecule {coordinates} + | a == b || a == c || a == d || b == c || b == d || c == d = throwM . IndexException $ "selected atoms are indentical" + | otherwise = do + coordA <- compute @U <$> (coordinates !?> a) + coordB <- compute @U <$> (coordinates !?> b) + coordC <- compute @U <$> (coordinates !?> c) + coordD <- compute @U <$> (coordinates !?> d) + vecAB <- coordA .-. coordB + vecBC <- coordB .-. coordC + vecCD <- coordC .-. coordD + let planeABC = vecH2M $ LA.cross (vecM2H vecAB) (vecM2H vecBC) :: Massiv.Vector U Double + planeBCD = vecH2M $ LA.cross (vecM2H vecBC) (vecM2H vecCD) :: Massiv.Vector U Double + normVecRot = LA.cross (vecM2H vecCD) (vecM2H vecBC) :: LA.Vector Double + rotDir = + if vecH2M normVecRot !.! vecAB < 0 + then -1 + else 1 + return $ rotDir * ND.angle planeABC planeBCD + +-- | Calculates the a metric value of the dihedral angle defined by four atoms. This must create 2 +-- values in the feature matrix, instead of one. +-- See +dihedral :: MonadThrow m => D -> Molecule -> m (Double, Double) +dihedral d mol = do + dihedRad <- dihedral' d mol + return (sin dihedRad, cos dihedRad) + +-- | Get all selected features from a molecule. +getFeature :: + ( Traversable t, + MonadThrow m + ) => + -- | Selection of multiple features. + t Feature -> + -- | A Molecule where to apply to. + Molecule -> + m (Massiv.Vector Massiv.DL Double) +getFeature sel mol = do + features <- for sel $ \s -> case s of + Energy -> return . Left . energy $ mol + Bond b -> Left <$> bond b mol + Angle a -> Left <$> angle a mol + Dihedral d -> Right <$> dihedral d mol + + let featureVec = + foldl + ( \acc val -> case val of + Left s -> acc `Massiv.snoc` s + Right (a, b) -> acc `Massiv.snoc` a `Massiv.snoc` b + ) + mempty + features + + return featureVec + +-- | Obtains the feature matrix \(\mathbf{X}\) for a principal component analysis. Given \(m\) +-- features to analyse in \(n\) measurements, \(\mathbf{X}\) will be a \(m \times n\) matrix. +{-# SCC getFeatures #-} +getFeatures :: + ( MonadThrow m, + Traversable f + ) => + f Feature -> + Trajectory -> + m (Massiv.Matrix DL Double) +getFeatures sel trj = traverse toCols trj >>= concatM (Dim 1) + where + toCols v = expandInner @Ix2 (Sz 1) const . compute @U <$> getFeature sel v diff --git a/src/ConClusion/Numeric/Data.hs b/src/ConClusion/Numeric/Data.hs new file mode 100644 --- /dev/null +++ b/src/ConClusion/Numeric/Data.hs @@ -0,0 +1,197 @@ +-- | +-- Module : ConClusion.Numeric.Data +-- Description : Type castings, conversions and utilities +-- Copyright : Phillip Seeber, 2021 +-- License : AGPL-3 +-- Maintainer : phillip.seeber@googlemail.com +-- Stability : experimental +-- Portability : POSIX, Windows +module ConClusion.Numeric.Data + ( -- * Conversion of array types. + + -- Conversion between HMatrix and Massiv's array types + IndexException (..), + vecH2M, + vecM2H, + matH2M, + matM2H, + + -- * Array Processing + magnitude, + normalise, + angle, + minDistAt, + minDistAtVec, + iMinimumM, + + -- * Utilities + printMat, + + -- * Binary Trees + BinTree (..), + root, + takeBranchesWhile, + takeLeafyBranchesWhile, + ) +where + +import Data.Aeson hiding (Array) +import Data.Massiv.Array as Massiv hiding (IndexException) +import Data.Massiv.Array.Manifest.Vector as Massiv +import Formatting +import Numeric.LinearAlgebra as LA hiding (magnitude, (<>)) +import RIO +import qualified RIO.Vector.Storable as VectorS +import System.IO.Unsafe (unsafePerformIO) + +-- | Exception regarding indexing in some kind of aaray. +newtype IndexException = IndexException String deriving (Show) + +instance Exception IndexException + +-- | Converts a vector from the HMatrix package to the Massiv representation. +{-# SCC vecH2M #-} +vecH2M :: (Element e, Mutable r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e +vecH2M hVec = fromVector' Seq (Sz $ VectorS.length hVec) hVec + +-- | Converts a vector from the Massiv representation to the HMatrix representation. +{-# SCC vecM2H #-} +vecM2H :: (Manifest r Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e +vecM2H = Massiv.toVector + +-- | Converts a matrix from the HMatrix representation to the Massiv representation. +{-# SCC matH2M #-} +matH2M :: (Mutable r Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e +matH2M hMat = Massiv.resize' (Sz $ nRows :. nCols) . vecH2M . LA.flatten $ hMat + where + nRows = LA.rows hMat + nCols = LA.cols hMat + +-- | Converts a matrix from Massiv to HMatrix representation. +{-# SCC matM2H #-} +matM2H :: (Manifest r Ix1 e, Element e, Resize r Ix2, Load r Ix2 e) => Massiv.Matrix r e -> LA.Matrix e +matM2H mMat = LA.reshape nCols . vecM2H . Massiv.flatten $ mMat + where + Sz (_nRows :. nCols) = Massiv.size mMat + +-- | Magnitude of a vector (length). +magnitude :: (Massiv.Numeric r e, Source r Ix1 e, Floating e) => Massiv.Vector r e -> e +magnitude v = sqrt $ v !.! v + +-- | Normalise a vector. +normalise :: (Massiv.Numeric r e, Source r Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e +normalise v = v .* (1 / magnitude v) + +-- | Angle between two vectors. +angle :: (Massiv.Numeric r e, Source r Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e +angle a b = acos $ a !.! b / (magnitude a * magnitude b) + +-- | Find the minimal distance in a distance matrix, which is not the main diagonal. +{-# SCC minDistAt #-} +minDistAt :: + ( Manifest r Ix2 e, + MonadThrow m, + Ord e + ) => + Massiv.Matrix r e -> + m (e, Ix2) +minDistAt arr + | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr) + | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr + where + ix0 = pureIndex 0 + e0 = arr Massiv.! ix0 + start = (e0, ix0) + minFold acc@(eA, _) ix@(m :. n) e = if e < eA && m > n then (e, ix) else acc + chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc + +-- | Find the minimal element of a vector, which is at a larger than the supplied index. +minDistAtVec :: + ( Manifest r Ix1 e, + MonadThrow m, + Ord e + ) => + Ix1 -> + Massiv.Vector r e -> + m (e, Ix1) +minDistAtVec ixStart vec + | isEmpty vec = throwM $ SizeEmptyException (Massiv.size vec) + | ixStart >= nElems = throwM $ IndexOutOfBoundsException (Sz nElems) ixStart + | otherwise = do + let (minE, minIx) = unsafePerformIO $ ifoldlP minFold startAcc chFold startAcc searchVec + return (minE, minIx + ixStart + 1) + where + Sz nElems = Massiv.size vec + searchVec = Massiv.drop (Sz $ ixStart + 1) vec + ix0 = 0 + e0 = searchVec Massiv.! ix0 + startAcc = (e0, ix0) + minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc + chFold acc ch = min acc ch + +-- | Like 'Massiv.minimumM' but also returns the index of the minimal element. +iMinimumM :: + ( Manifest r ix a, + MonadThrow m, + Ord a + ) => + Array r ix a -> + m (a, ix) +iMinimumM arr + | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr) + | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr + where + ix0 = pureIndex 0 + e0 = arr Massiv.! ix0 + start = (e0, ix0) + + minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc + chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc + +-- | Quickly print a matrix with numerical values +printMat :: (Source r Ix2 e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text +printMat mat = Massiv.map (sformat (left 4 ' ' %. fixed 2)) mat + +---------------------------------------------------------------------------------------------------- +-- Binary Trees. + +-- | A binary tree. +data BinTree e = Leaf e | Node e (BinTree e) (BinTree e) + deriving (Eq, Show, Generic) + +instance (FromJSON e) => FromJSON (BinTree e) + +instance (ToJSON e) => ToJSON (BinTree e) + +instance Functor BinTree where + fmap f (Leaf a) = Leaf (f a) + fmap f (Node a l r) = Node (f a) (fmap f l) (fmap f r) + +-- | Look at the root of a binary tree. +root :: BinTree e -> e +root (Leaf e) = e +root (Node e _ _) = e + +-- | Steps down each branch of a tree until some criterion is satisfied or the end of the branch is +-- reached. Each end of the branch is added to a result. +takeBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a +takeBranchesWhile chk tree = go tree (Massiv.empty @DL) + where + go (Leaf v) acc = if chk v then acc `snoc` v else acc + go (Node v l r) acc = + let vAcc = if chk v then acc `snoc` v else acc + lAcc = go l vAcc + rAcc = go r lAcc + in if chk v then rAcc else vAcc + +-- | Takes the first value in each branch, that does not fullfill the criterion anymore and adds it +-- to the result. Terminal leafes of the branches are always taken. +takeLeafyBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a +takeLeafyBranchesWhile chk tree = go tree (Massiv.empty @DL) + where + go (Leaf v) acc = acc `snoc` v + go (Node v l r) acc = + let vAcc = if chk v then acc else acc `snoc` v + lAcc = go l vAcc + rAcc = go r lAcc + in if chk v then rAcc else vAcc diff --git a/src/ConClusion/Numeric/Statistics.hs b/src/ConClusion/Numeric/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/ConClusion/Numeric/Statistics.hs @@ -0,0 +1,947 @@ +-- | +-- Module : ConClusion.Numeric.Statistics +-- Description : Statistical Functions +-- Copyright : Phillip Seeber, 2021 +-- License : AGPL-3 +-- Maintainer : phillip.seeber@googlemail.com +-- Stability : experimental +-- Portability : POSIX, Windows +module ConClusion.Numeric.Statistics + ( -- * PCA + PCA (..), + pca, + + -- * Variance + normalise, + meanDeviation, + covariance, + + -- * Distance Metrics + DistFn, + lpNorm, + manhattan, + euclidean, + mahalanobis, + + -- * Cluster Algorithms + Clusters, + + -- ** DBScan + DistanceInvalidException (..), + dbscan, + + -- ** Hierarchical Cluster Analysis + Dendrogram, + JoinStrat (..), + hca, + cutDendroAt, + ) +where + +import ConClusion.Numeric.Data hiding (normalise) +import Data.Aeson hiding (Array) +import Data.Complex +import qualified Data.IntSet as IntSet +import Data.Massiv.Array as Massiv +import Data.Massiv.Array.Unsafe as Massiv +import qualified Data.PSQueue as PQ +import qualified Numeric.LinearAlgebra as LA +import RIO hiding (Vector) +import System.IO.Unsafe (unsafePerformIO) + +---------------------------------------------------------------------------------------------------- +-- Others/Helpers + +-- | Solves eigenvalue problem of a square matrix and obtains its eigenvalues and eigenvectors. +{-# SCC eig #-} +eig :: + ( Mutable r1 Ix1 (Complex Double), + Mutable r2 Ix1 (Complex Double), + LA.Field e, + Manifest r3 Ix1 e, + Resize r3 Ix2, + Load r3 Ix2 e, + MonadThrow m + ) => + Matrix r3 e -> + m (Vector r1 (Complex Double), Matrix r2 (Complex Double)) +eig covM + | m /= n = throwM $ IndexException "eigenvalue problems can only be solved for square matrix" + | otherwise = return . bimap vecH2M matH2M . LA.eig $ cov + where + Sz (m :. n) = size covM + cov = matM2H covM + +-- | Sort eigenvalues and eigenvectors by magnitude of the eigenvalues in descending order (largest +-- eigenvalues first). Eigenvectors are the columns of the input matrix. +{-# SCC eigSort #-} +eigSort :: + ( Load r2 Ix2 e, + MonadThrow m, + Source r1 Ix1 e, + Source r2 Ix2 e, + Mutable r1 Ix1 e, + Mutable r2 Ix2 e, + Unbox e, + Ord e + ) => + (Vector r1 e, Matrix r2 e) -> + m (Vector r1 e, Matrix r2 e) +eigSort (vec, mat) + | m /= n = throwM $ IndexException "matrix of the eigenvectors is not a square matrix" + | n /= n' = throwM $ IndexException "different number of eigenvalues and eigenvectors" + | otherwise = do + let ixedEigenvalues = Massiv.zip vec ixVec + (eigenValSortAsc, ixSort) = (\a -> (get fst a, get snd a)) . quicksort . compute @U $ ixedEigenvalues + eigenVecSortAsc = backpermute' (Sz $ m :. n) (\(r :. c) -> r :. (ixSort ! c)) mat + eigenValSort = reverse' (Dim 1) eigenValSortAsc + eigenVecSort = reverse' (Dim 1) eigenVecSortAsc + return (compute eigenValSort, compute eigenVecSort) + where + Sz (m :. n) = size mat + Sz n' = size vec + ixVec = makeArrayLinear @D Seq (Sz n') id + get acc = compute @U . Massiv.map acc + +---------------------------------------------------------------------------------------------------- +-- Principal Component Analysis + +data PCA = PCA + { -- | Original feature matrix. + x :: Matrix U Double, + -- | Feature matrix in mean deviation form. + x' :: Matrix U Double, + -- | Transformed data. + y :: Matrix U Double, + -- | Transformation matrix to transform feature matrix into PCA result matrix. + a :: Matrix U Double, + -- | Mean squared error introduced by PCA. + mse :: Double, + -- | Percentage of the behaviour captured in the remaining dimensions. + remaining :: Double, + -- | All eigenvalues from the diagonalisation of the covariance matrix. + allEigenValues :: Vector U Double, + -- | Eigenvalues that were kept for PCA. + pcaEigenValues :: Vector U Double, + -- | All eigenvectors from the diagonalisation of the covariance matrix. + allEigenVecs :: Matrix U Double, + -- | Eigenvectors that were kept for PCA. + pcaEigenVecs :: Matrix U Double + } + +-- | Transform the input values with a transformation matrix \(\mathbf{A}\), where \(\mathbf{A}\) is +-- constructed from the eigenvectors associated to the largest eigenvalues. +{-# SCC transformToPCABasis #-} +transformToPCABasis :: + ( Source (R r) Ix2 e, + Extract r Ix2 e, + Mutable r Ix2 e, + Numeric r e, + MonadThrow m + ) => + -- | Number of dimensions to keep from PCA. + Int -> + -- | Matrix of the eigenvectors, sorted descendingly by eigenvalues, where the eigenvectors are + -- the columns of the matrix. + Matrix r e -> + -- | Feature matrix in mean deviation form. + Matrix r e -> + -- | Input data transformed by PCA to lower dimensions, and the transformation matrix + -- \(\mathbf{A}\). + m (Matrix r e, Matrix r e) +transformToPCABasis nDim eigenVecMat featureMat + | mE /= nE = throwM $ IndexException "the matrix of the eigenvectors must be a quadratic matrix" + | nDim <= 0 = throwM $ IndexException "the number of dimensions of the PCA is smaller than or zero" + | nDim >= nE = throwM $ IndexException "more than the possible amount of dimensions has been selected" + | mE /= mF = throwM $ IndexException "eigenvector matrix and feature matrix have mismatching dimensions" + | otherwise = do + matA <- compute . transpose <$> extractM (0 :. 0) (Sz $ mE :. nDim) eigenVecMat + pcaData <- matA .><. featureMat + return (pcaData, matA) + where + Sz (mE :. nE) = size eigenVecMat + Sz (mF :. _nF) = size featureMat + +-- | Performs a PCA on the feature matrix \(\mathbf{X}\) by solving the eigenproblem of the +-- covariance matrix. The function takes the feature matrix directly and perfoms the conversion +-- to mean deviation form, the calculation of the covariance matrix and the eigenvalue problem +-- automatically. +{-# SCC pca #-} +pca :: + ( Numeric r Double, + Mutable r Ix2 Double, + Manifest r Ix1 Double, + Source (R r) Ix2 Double, + Extract r Ix2 Double, + MonadThrow m + ) => + -- | Dimensionalty after PCA transformation. + Int -> + -- | \(m \times n\) Feaute matrix \(\mathbf{X}\), with \(m\) different measurements (rows) in + -- \(n\) different trials (columns). + Matrix r Double -> + m PCA +pca dim x = do + -- Calculate the mean deviation form of the feature matrix and the covariance matrix from it. + let x' = normalise . meanDeviation $ x + cov = covariance x' + + -- Obtain eigenvalues and eigenvectors of the covariance matrix and sort them. + (eigValsC :: Vector U (Complex Double), eigVecsC :: Matrix U (Complex Double)) <- eig cov + let eigValsR = compute @U . Massiv.map realPart $ eigValsC + eigVecsR = compute . Massiv.map realPart $ eigVecsC + (eValS, eVecS) <- eigSort (eigValsR, eigVecsR) + + -- Use the subset of the eigenvectors with the largest eigenvalues to transform the features in + -- mean deviation form into the result matrix Y. + (pcaData, matA) <- transformToPCABasis dim eVecS x' + + -- Reconstuct the original data from lower dimensions and calculate the mean squared deviation. + reconstructX <- (compute . transpose $ matA) .><. pcaData + mse <- (/ fromIntegral n) . Massiv.sum . Massiv.map (** 2) <$> (x' .-. reconstructX) + + -- For output give the eigenvalues and eigenvectors that were kept. + pcaEigenValues <- extractM 0 (Sz dim) eValS + pcaEigenVecs <- extractM (0 :. 0) (Sz $ m :. dim) eVecS + + -- Calculate the amount of behaviour that could be kept. + let remaining = (Massiv.sum pcaEigenValues / Massiv.sum eValS) * 100 + + return + PCA + { x = compute x, + x' = compute x', + y = compute pcaData, + a = compute matA, + mse = mse, + remaining = remaining, + allEigenValues = eValS, + pcaEigenValues = compute pcaEigenValues, + allEigenVecs = compute eVecS, + pcaEigenVecs = compute pcaEigenVecs + } + where + Sz (m :. n) = size x + +---------------------------------------------------------------------------------------------------- +-- Variance + +-- | Subtract the mean value of all columns from the feature matrix. Brings the feature matrix to +-- mean deviation form. +{-# SCC meanDeviation #-} +meanDeviation :: + ( Source r Ix2 e, + Fractional e, + Unbox e, + Numeric r e, + Mutable r Ix2 e + ) => + Matrix r e -> + Matrix r e +meanDeviation mat = mat !-! compute meanMat + where + Sz (_ :. n) = Massiv.size mat + featueMean = Massiv.foldlInner (+) 0 mat .* (1 / fromIntegral n) + meanMat = expandInner (Sz n) const . compute @U $ featueMean + +-- | Obtains the covariance matrix \(\mathbf{C_X}\) from the feature matrix \(\mathbf{X}\). +-- \[ +-- \mathbf{C_X} \equiv \frac{1}{n - 1} \mathbf{X} \mathbf{X}^T +-- \] +-- where \(n\) is the number of columns in the matrix. +-- +-- The feature matrix should be in mean deviation form, see 'meanDeviation'. +{-# SCC covariance #-} +covariance :: (Numeric r e, Mutable r Ix2 e, Fractional e) => Matrix r e -> Matrix r e +covariance x = (1 / (fromIntegral n - 1)) *. (x !> + Array r Ix2 e -> + Array r Ix2 e +normalise mat = + let absMat = Massiv.map abs mat + maxPerRow = compute @U . foldlInner max 0 $ absMat + divMat = compute . Massiv.map (1 /) . expandInner @Ix2 (Sz n) const $ maxPerRow + in divMat !*! mat + where + Sz (_ :. n) = size mat + +---------------------------------------------------------------------------------------------------- +-- Distance Measures + +-- | Distance matrix generator functions. +type DistFn r e = Matrix r e -> Matrix r e + +-- | Builds the distance measures in a permutation matrix/distance matrix. +buildDistMat :: + (Mutable r Ix2 e) => + -- | Zip function to combine the elements of vectors \(\mathbf{a}\) \(\mathbf{b}\). Usually @(-)@. + -- \( f(\mathbf{a}_i, \mathbf{b}_i) = \mathbf{c} \) + (e -> e -> a) -> + -- | Fold the vector \(\mathbf{c}\) elementwise to a distance \(d\). + (a -> a -> a) -> + -- | Accumulator of the fold function. + a -> + -- | \(m \times n\) matrix, with \(n\) \(m\)-dimensional points (column vectors of the matrix). + Matrix r e -> + -- | Resulting distance matrix. + Matrix D a +buildDistMat zipFn foldFn acc mat = + let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ mat + b = transposeInner a + ab = Massiv.zipWith zipFn a b + d = foldlInner foldFn acc ab + in d + where + Sz (_ :. n) = size mat + +-- | The \(L_p\) norm between two vectors. Generalisation of Manhattan and Euclidean distances. +-- \[ +-- d(\mathbf{a}, \mathbf{b}) = \left( \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert ^p \right) ^ \frac{1}{p} +-- \] +{-# SCC lpNorm #-} +lpNorm :: (Mutable r Ix2 e, Floating e) => Int -> DistFn r e +lpNorm p = compute . buildDistMat zipFn foldFn 0 + where + zipFn a b = (^ p) . abs $ a - b + foldFn a b = (** (1 / fromIntegral p)) $ a + b + +-- | The Manhattan distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 1\). +-- \[ +-- d(\mathbf{a}, \mathbf{b}) = \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert +-- \] +{-# SCC manhattan #-} +manhattan :: (Mutable r Ix2 e, Floating e) => DistFn r e +manhattan = lpNorm 1 + +-- | The Euclidean distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 2\). +-- \[ +-- d(\mathbf{a}, \mathbf{b}) = \sqrt{\sum \limits_{i=1}^n (\mathbf{a}_i - \mathbf{b}_i)^2} +-- \] +{-# SCC euclidean #-} +euclidean :: (Mutable r Ix2 e, Floating e) => DistFn r e +euclidean = lpNorm 2 + +-- | Mahalanobis distance between points. Suitable for non correlated axes. +-- \[ +-- d(\mathbf{a}, \mathbf{b}) = \sqrt{(\mathbf{a} - \mathbf{b})^T \mathbf{S}^{-1} (\mathbf{a} - \mathbf{b})} +-- \] +-- where \(\mathbf{S}\) is the covariance matrix. +{-# SCC mahalanobis #-} +mahalanobis :: (Unbox e, Numeric r e, Mutable r Ix2 e, Mutable r Ix1 e, LA.Field e) => DistFn r e +mahalanobis points = + let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ points + b = transposeInner a + abDiff = compute @U $ a !-! b + ixArray = makeArray @U @Ix2 @Ix2 Par (Sz $ n :. n) id + distMat = + Massiv.map + ( \(x :. y) -> + let ab = compute @U $ abDiff !> x !> y + in ab > Exception (DistanceInvalidException e) + +-- | Representation of clusters. +type Clusters = Vector B IntSet + +-- | DBScan algorithm. +{-# SCC dbscan #-} +dbscan :: + ( MonadThrow m, + Ord e, + Num e, + Typeable e, + Show e, + Source r Ix2 e + ) => + -- | Distance measure to build the distance matrix of all points. + DistFn r e -> + -- | Minimal number of members in a cluster. + Int -> + -- | Search radius \(\epsilon\) + e -> + -- | \(n\) \(m\)-dimensional data points as column vectors of a \(m \times n\) matrix. + Matrix r e -> + -- | Resulting clusters. + m Clusters +dbscan distFn nPoints epsilon points + | isEmpty points = throwM $ SizeEmptyException (Sz 0 :: Sz1) + | nPoints < 1 = throwM $ SizeNegativeException (Sz nPoints) + | epsilon <= 0 = throwM $ DistanceInvalidException epsilon + | otherwise = + let pointNeighbours = ifoldlInner collectNeighbours mempty distMat + allClusters = joinOverlapping . compute @B $ pointNeighbours + largeClusters = sfilter (\s -> IntSet.size s >= nPoints) allClusters + in return $ compute largeClusters + where + -- The distance matrix in the measure chosen by the distance function. + distMat = distFn points + + -- Function to collect the neighbours of a point within the search radius epsilon. + {-# SCC collectNeighbours #-} + collectNeighbours (_ :. n) acc d = if d <= epsilon then IntSet.insert n acc else acc + + -- Construct the overlap matrix of all clusters. + compareSets :: (IntSet -> IntSet -> Bool) -> Vector B IntSet -> Matrix D Bool + compareSets fn clVec = + let a = expandOuter @Ix2 sz const clVec + b = transpose a + compareMat = Massiv.zipWith fn a b + in compareMat + where + sz = size clVec + + -- Overlap matrix. Checks if two sets have any overlap. Sets do overlap with themself. + overlap :: Vector B IntSet -> Matrix D Bool + overlap = compareSets (\a b -> not $ IntSet.disjoint a b) + + -- Check if any set overlaps wiht **any** other set. + anyOtherOverlap :: Vector B IntSet -> Bool + anyOtherOverlap = Massiv.any (== True) . imap (\(m :. n) v -> if m == n then False else v) . overlap + + -- Check if two sets are identical. Sets are identical to themself. + same :: Vector B IntSet -> Matrix D Bool + same = compareSets (==) + + -- Join all overlapping clusters recursively. + {-# SCC joinOverlapping #-} + joinOverlapping :: Vector B IntSet -> Vector B IntSet + joinOverlapping clVec = + let -- The overlap matrix of the clusters. + ovlpMat = compute @U . overlap $ clVec + anyOvlp = anyOtherOverlap clVec + + -- Join all sets that have overlap but keep them redundantly (no reduction of the amount + -- of clusters). + joined = + ifoldlInner + (\(_ :. n) acc ovlp -> if ovlp then (clVec ! n) <> acc else acc) + mempty + ovlpMat + + -- Find all sets at different indices that are the same. This is an upper triangular + -- matrix with the main diagonal being False. + sameMat = + compute @U + . imap (\(m :. n) v -> if m >= n then False else v) + . same + . compute @B + $ joined + + -- Remove all sets that are redundant. Redundancy is checked by two criteria: + -- 1. Is this cluster the same set of points as any other cluster? If yes, it is + -- redundant. + -- 2. Is this cluster isolated it is not redundant. + nonRed = + compute @B + . sifilter + ( \ix _ -> + let sameAsAnyOther = Massiv.any (== True) $ sameMat !> ix + in not sameAsAnyOther + ) + $ joined + in if anyOvlp then joinOverlapping nonRed else clVec + +---------------------------------------------------------------------------------------------------- +-- Hierarchical Cluster Analysis + +-- | Nodes of a dendrogram. +data DendroNode e = DendroNode + { distance :: e, + cluster :: IntSet + } + deriving (Eq, Show, Generic) + +instance (FromJSON e) => FromJSON (DendroNode e) + +instance (ToJSON e) => ToJSON (DendroNode e) + +-- | A dendrogram as a binary tree. +newtype Dendrogram e = Dendrogram {unDendro :: BinTree (DendroNode e)} + deriving (Show, Eq, Generic) + +instance ToJSON e => ToJSON (Dendrogram e) + +instance FromJSON e => FromJSON (Dendrogram e) + +-- | An accumulator to finally build a dendrogram by a bottom-up algorithm. Not to be exposed in the +-- API. +type DendroAcc e = Vector B (Dendrogram e) + +-- | Mutable version of the dendrogram accumulator. +type DendroAccM m e = MArray (PrimState m) B Ix1 (Dendrogram e) + +-- | Cut a 'Dendrogram' at a given distance and obtain all clusters from it. +cutDendroAt :: Ord e => Dendrogram e -> e -> Clusters +cutDendroAt dendro dist = + let nodes = takeLeafyBranchesWhile (\DendroNode {distance} -> distance >= dist) . unDendro $ dendro + clusters = compute @B . Massiv.map cluster . compute @B $ nodes + in clusters + +-- | A strategy/distance measure for clusters. +data JoinStrat e + = SingleLinkage + | CompleteLinkage + | Median + | UPGMA + | WPGMA + | Centroid + | Ward + | LWFB e + | LW e e e e + deriving (Eq, Show) + +-- | Lance Williams formula to update distances. +{-# SCC lanceWilliams #-} +lanceWilliams :: + Fractional e => + -- | How to calculate distance between clusters of points. + JoinStrat e -> + -- | Number of points in cluster \(A\). + Int -> + -- | Number of points in cluster \(B\) + Int -> + -- | Number of points in cluster \(C\) + Int -> + -- | \(d(A, B)\) + e -> + -- | \(d(A, C)\) + e -> + -- | \(d(B, C)\) + e -> + -- | Updated distance \(D \(A \cup B, C\) + e +lanceWilliams js nA nB nC dAB dAC dBC = alpha1 * dAC + alpha2 * dBC + beta * dAB + gamma * abs (dAC - dBC) + where + nA' = fromIntegral nA + nB' = fromIntegral nB + nC' = fromIntegral nC + (alpha1, alpha2, beta, gamma) = case js of + SingleLinkage -> (1 / 2, 1 / 2, 0, - 1 / 2) + CompleteLinkage -> (1 / 2, 1 / 2, 0, 1 / 2) + Median -> (1 / 2, 1 / 2, - 1 / 4, 0) + UPGMA -> (nA' / (nA' + nB'), nB' / (nA' + nB'), 0, 0) + WPGMA -> (1 / 2, 1 / 2, 0, 0) + Centroid -> (nA' / (nA' + nB'), nB' / (nA' + nB'), - (nA' * nB') / ((nA' + nB') ^ (2 :: Int)), 0) + Ward -> ((nA' + nC') / (nA' + nB' + nC'), (nA' + nC') / (nA' + nB' + nC'), - (nA' + nC') / (nA' + nB' + nC'), 0) + LWFB b -> ((1 - b) / 2, (1 - b) / 2, b, 0) + LW a b c d -> (a, b, c, d) + +---------------------------------------------------------------------------------------------------- +-- Müllner Generic Hierarchical Clustering + +-- | A neighbourlist. At index @i@ of the vector it contains a tuple with the minimal distance of +-- this cluster to any other cluster and the index of the other cluster. +type Neighbourlist r e = Vector r (e, Ix1) + +-- | A distance matrix. +type DistanceMatrix r e = Matrix r e + +-- | Performance improved hierarchical clustering algorithm. @GENERIC_LINKAGE@ from figure 3, +-- . +{-# SCC hca #-} +hca :: + ( MonadThrow m, + Mutable r Ix1 e, + Mutable r Ix2 e, + Mutable r Ix1 (e, Ix1), + Manifest (R r) Ix1 e, + OuterSlice r Ix2 e, + Ord e, + Unbox e, + Fractional e + ) => + DistFn r e -> + JoinStrat e -> + Matrix r e -> + m (Dendrogram e) +hca distFn joinStrat points + | Massiv.isEmpty points = throw $ SizeEmptyException (Sz nPoints) + | otherwise = do + let -- The distance matrix from the points. + distMat = distFn points + + -- Initial vector of nearest neighbour to each point. + nNghbr <- nearestNeighbours distMat + + let -- Initial priority queue of points. Has the minimum distance of all points. + pq = PQ.fromList . Massiv.toList . Massiv.imap (\k (d, _) -> k PQ.:-> d) $ nNghbr + -- Set of points not joined yet. Initially all points. + s = IntSet.fromDistinctAscList [0 .. nPoints - 1] + -- Initial dendrogram accumulator. The vector of all points as their own cluster. + dendroAcc = + makeArray @B @Ix1 + Par + (Sz nPoints) + (\p -> Dendrogram . Leaf $ DendroNode {distance = 0, cluster = IntSet.singleton p}) + + distMatM <- return . unsafePerformIO . thaw $ distMat + nNghbrM <- return . unsafePerformIO . thaw $ nNghbr + dendroAccM <- return . unsafePerformIO . thaw $ dendroAcc + + return . unsafePerformIO $ agglomerate joinStrat distMatM nNghbrM pq s dendroAccM + where + Sz (_mFeatures :. nPoints) = size points + +-- | Agglomerative clustering by the improved generic linkage algorithm. This is the main loop +-- recursion L 10-43. +{-# SCC agglomerate #-} +agglomerate :: + ( MonadThrow m, + PrimMonad m, + MonadUnliftIO m, + PrimState m ~ RealWorld, + Mutable r Ix2 e, + OuterSlice r Ix2 e, + Manifest (R r) Ix1 e, + Mutable r Ix1 (e, Ix1), + Fractional e, + Ord e + ) => + -- | Join strategy for clusters and therefore how to calculate cluster-cluster distances. + JoinStrat e -> + -- | Distance matrix. + MArray (PrimState m) r Ix2 e -> + -- | List of nearest neighbours for each point. + MArray (PrimState m) r Ix1 (e, Ix1) -> + -- | Priority queue with the distances as priorities and the cluster index as keys. + PQ.PSQ Ix1 e -> + -- | A set \(S\), that keeps track which clusters have already been joined. + IntSet -> + -- | Accumulator of the dendrogram. Should collapse to a singleton vector. + DendroAccM m e -> + -- | The final dendrogram, after all clusters have been joined. + m (Dendrogram e) +agglomerate joinStrat distMat nNghbr pq s dendroAcc + | IntSet.null s = throw $ IndexException "No clusters left. This must never happen." + | otherwise = do + -- Obtain candidates for the two clusters to join and the minimal distance in the priority queue. + candidates <- getJoinCandidates nNghbr pq + + -- If the distance between a b is not the minimal distance that the priority queue has found, the + -- neighbour list must be wrong and recalculated. + (a, b, delta, nNghbrU1, pqU1) <- recalculateNghbr candidates s distMat nNghbr pq + + -- Remove the minimal element from the priority queue and join clusters a and b. The cluster + -- accumulator is reduced in its size: a is removed and b is updated with the joined cluster. + (newS, pqU2, newAcc) <- joinClusters a b delta s pqU1 dendroAcc + + -- Update the distance matrix in the row and column of b but not at (b,b) and not at (a,b) and + -- (b,a). + newDistMat <- updateDistMat joinStrat a b newS distMat newAcc + + -- Redirect neighbours to b, if they previously pointed to a. + nNghbrU2 <- redirectNeighbours a b newS newDistMat nNghbrU1 + + -- Update the neighbourlist and priority queue with the new distances to b. + (newNNghbr, newPQ) <- updateBNeighbour b s newDistMat nNghbrU2 pqU2 + + -- If the problem has been reduced to a single cluster the algorithm is done and the final + -- dendrogram can be obtained from the accumulator at index b. Otherwise join further. + if IntSet.size newS == 1 + then newAcc `readM` b + else agglomerate joinStrat newDistMat newNNghbr newPQ newS newAcc + +-- | Obtain candidates for the clusters to join by looking at the minimal distance in the priority +-- queue and the neighbourlist. L 11-13 +{-# SCC getJoinCandidates #-} +getJoinCandidates :: + ( MonadThrow m, + PrimMonad m, + Mutable r Ix1 (e, Ix1), + Ord e + ) => + MArray (PrimState m) r Ix1 (e, Ix1) -> + PQ.PSQ Ix1 e -> + m (Ix1, Ix1, e) +getJoinCandidates nNghbr pq = do + (a PQ.:-> d) <- case PQ.findMin pq of + Nothing -> throwM $ IndexException "Empty priority queue" + Just v -> return v + (_, b) <- nNghbr `readM` a + return (a, b, d) + +-- | If the minimal distance @d@ found is not the distance between @a@ and @b@ recalculate the +-- neighbour list, update the priority queue and obtain a new set of a,b and a distance between them. +-- L 14-20. +{-# SCC recalculateNghbr #-} +recalculateNghbr :: + ( MonadThrow m, + PrimMonad m, + MonadUnliftIO m, + PrimState m ~ RealWorld, + OuterSlice r Ix2 e, + Manifest (R r) Ix1 e, + Mutable r Ix1 (e, Ix1), + Mutable r Ix2 e, + Ord e + ) => + (Ix1, Ix1, e) -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + MArray (PrimState m) r Ix1 (e, Ix1) -> + PQ.PSQ Ix1 e -> + m (Ix1, Ix1, e, MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e) +recalculateNghbr (cA, cB, d) s distMat nNghbr pq = do + dAB <- distMat `readM` (cA :. cB) + if d == dAB + then return (cA, cB, d, nNghbr, pq) + else do + -- Recalculate the nearest neighbours just on index cA. Consider only clusters, that were not + -- merged yet. + dmRowA <- searchRow cA s distMat >>= unsafeFreeze Par + newNeighbourA@(minDistA, _) <- minimumM dmRowA + writeM nNghbr cA newNeighbourA + + -- Update the priority queue at key cA with the new distance. + let newPQ = PQ.adjust (const minDistA) cA pq + + -- Determine new a, b and d from the updated neighbour list and priority queue. + (a PQ.:-> newD) <- case PQ.findMin newPQ of + Nothing -> throwM $ IndexException "Empty priority queue" + Just v -> return v + (_, b) <- nNghbr `readM` a + recalculateNghbr (a, b, newD) s distMat nNghbr newPQ + +-- | Joins the selected clusters \(A\) and \(B\) and updates the dendrogram accumulator at index b. +-- A will not be removed so that the accumulator never shrinks. +-- L 21-24 +{-# SCC joinClusters #-} +joinClusters :: + ( MonadThrow m, + PrimMonad m, + Ord e + ) => + Ix1 -> + Ix1 -> + e -> + IntSet -> + PQ.PSQ Ix1 e -> + DendroAccM m e -> + m (IntSet, PQ.PSQ Ix1 e, DendroAccM m e) +joinClusters a b d s pq acc = do + clA <- acc `readM` a + let newPQ = PQ.deleteMin pq + modifyM_ + acc + ( \clB -> + return + . Dendrogram + $ Node + ( DendroNode + { distance = d, + cluster = (cluster . root . unDendro $ clA) <> (cluster . root . unDendro $ clB) + } + ) + (unDendro clA) + (unDendro clB) + ) + b + let newS = IntSet.delete a s + return (newS, newPQ, acc) + +-- | Update the distance matrix with a Lance-Williams update in the rows and columns of cluster b. +-- L 25-27 +{-# SCC updateDistMat #-} +updateDistMat :: + ( MonadThrow m, + PrimMonad m, + MonadUnliftIO m, + Mutable r Ix2 e, + Fractional e + ) => + JoinStrat e -> + Ix1 -> + Ix1 -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + DendroAccM m e -> + m (MArray (PrimState m) r Ix2 e) +updateDistMat js a b s distMat dendroAcc + | nDM /= nDM = throwM $ SizeMismatchException (Sz nDM) (Sz nCl) + | mDM /= nDM = throwM $ SizeMismatchException (Sz mDM) (Sz nDM) + | otherwise = do + dAB <- distMat `readM` (a :. b) + nA <- clSize a + nB <- clSize b + forIO_ ixV $ \ix -> do + dAX <- distMat `readM` (a :. ix) + nX <- clSize ix + modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (ix :. b) + modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (b :. ix) + return distMat + where + Sz (mDM :. nDM) = msize distMat + Sz nCl = msize dendroAcc + ixV = Massiv.fromList @U Par . IntSet.toAscList . IntSet.delete b $ s + clSize i = IntSet.size . cluster . root . unDendro <$> dendroAcc `readM` i + +-- | Updates the neighbourlist. All elements with a smaller index than a, that had a as a nearest +-- neighbour are blindly redirected to the union of a and b, now at index b. +-- L 28-32 +{-# SCC redirectNeighbours #-} +redirectNeighbours :: + ( MonadThrow m, + PrimMonad m, + MonadUnliftIO m, + Mutable r Ix1 (e, Ix1), + Mutable r Ix2 e + ) => + Ix1 -> + Ix1 -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + MArray (PrimState m) r Ix1 (e, Ix1) -> + m (MArray (PrimState m) r Ix1 (e, Ix1)) +redirectNeighbours a b s distMat nNghbr = do + forIO_ ixV $ \ix -> + modifyM_ + nNghbr + ( \old@(_, nghbrX) -> + if nghbrX == a + then distMat `readM` (ix :. b) >>= \dXB -> return (dXB, b) + else return old + ) + ix + return nNghbr + where + ixV = compute @U . sfilter (< a) . Massiv.fromList @U Par . IntSet.toAscList $ s + +-- | Updates the list of nearest neighbours for all combinations that might have changed by +-- recalculation with the joined cluster AB at index b. +-- L +{-# SCC updateWithNewBDists #-} +updateWithNewBDists :: + ( MonadThrow m, + MonadUnliftIO m, + PrimMonad m, + Mutable r Ix2 e, + Mutable r Ix1 (e, Ix1), + Ord e + ) => + Ix1 -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + MArray (PrimState m) r Ix1 (e, Ix1) -> + PQ.PSQ Ix1 e -> + m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e) +updateWithNewBDists b s distMat nNghbr pq = do + pqT <- newTVarIO pq + forIO_ ixV $ \ix -> do + dBX <- distMat `readM` (ix :. b) + currentPQ <- readTVarIO pqT + minDistX <- case PQ.lookup ix currentPQ of + Nothing -> throwM $ IndexException "Empty priority queue." + Just v -> return v + if dBX < minDistX + then do + writeM nNghbr ix (dBX, b) + atomically . writeTVar pqT . PQ.adjust (const dBX) ix $ currentPQ + else atomically . writeTVar pqT $ currentPQ + + newPQ <- readTVarIO pqT + return (nNghbr, newPQ) + where + ixV = compute @U . Massiv.sfilter (< b) . Massiv.fromList @U Par . IntSet.toAscList $ s + +-- | Updates the list of nearest neighbours and the priority queue at key b. +-- L 39-40 +{-# SCC updateBNeighbour #-} +updateBNeighbour :: + ( MonadThrow m, + PrimMonad m, + MonadUnliftIO m, + Mutable r Ix1 (e, Ix1), + Mutable r Ix2 e, + Ord e + ) => + Ix1 -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + MArray (PrimState m) r Ix1 (e, Ix1) -> + PQ.PSQ Ix1 e -> + m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e) +updateBNeighbour b s distMat nNghbr pq = + if b >= nNeighbours + then return (nNghbr, pq) + else do + rowAB <- searchRow b s distMat >>= unsafeFreeze Par + newNeighbourB@(distB, neighbourB) <- minimumM rowAB + writeM nNghbr b newNeighbourB + let newPQ = PQ.adjust (const distB) neighbourB pq + return (nNghbr, newPQ) + where + Sz nNeighbours = msize nNghbr + +-- | Find the nearest neighbour for each point from a distance matrix. For each point it stores the +-- minimum distance and the index of the other point, that is the nearest neighbour but at a higher +-- index. +{-# SCC nearestNeighbours #-} +nearestNeighbours :: + ( MonadThrow m, + Mutable r Ix1 e, + Mutable r Ix1 (e, Ix1), + OuterSlice r Ix2 e, + Source (R r) Ix1 e, + Ord e, + Unbox e + ) => + Matrix r e -> + m (Vector r (e, Ix1)) +nearestNeighbours distMat + | m /= n = throwM $ IndexException "Distance matrix is not square" + | m == 0 = throwM $ IndexException "Distance matrix is empty" + | otherwise = + let rows = compute @B . outerSlices $ distMat + minDistIx = + Massiv.imap (\i v -> unsafePerformIO . minDistAtVec i . compute @U $ v) . init $ rows + in return . compute $ minDistIx + where + Sz (m :. n) = size distMat + +-- Make a search row for distances. Takes row x from a distance matrix and zips them with their +-- column index. Then keeps only the valid elements of the row, that are still part of the available +-- points. A minimum or maximum search can be performed on the resulting vector and a valid pair of +-- distance and index can be obtained. +searchRow :: + ( PrimMonad m, + MonadThrow m, + MonadUnliftIO m, + Mutable r Ix2 e, + Mutable r Ix1 (e, Ix1) + ) => + Ix1 -> + IntSet -> + MArray (PrimState m) r Ix2 e -> + m (MArray (PrimState m) r Ix1 (e, Ix1)) +searchRow x s dm = + makeMArray Par (size ixV) $ \ix -> do + dmIx <- ixV !? ix + (dm `readM` (x :. dmIx)) >>= \dist -> return (dist, dmIx) + where + ixV :: Vector U Ix1 + ixV = compute @U . sfilter (> x) . Massiv.fromList @U Par . IntSet.toAscList $ s