[med-svn] [spaced] 01/06: Imported Upstream version 1.0.2
Fabian Klötzl
kloetzl-guest at moszumanska.debian.org
Wed Mar 23 12:26:06 UTC 2016
This is an automated email from the git hooks/post-receive script.
kloetzl-guest pushed a commit to branch master
in repository spaced.
commit 6fe59c5d336af6ac42625ee76626915eae607459
Author: Fabian Klötzl <fabian at kloetzl.info>
Date: Wed Mar 23 11:26:53 2016 +0100
Imported Upstream version 1.0.2
---
COPYING | 674 ++++++++++++++++++++++++++++++++++++
Makefile | 4 +
README | 37 ++
extkey.cpp | 73 ++++
extkey.h | 41 +++
patternset.cpp | 1050 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
patternset.h | 103 ++++++
sort.h | 694 +++++++++++++++++++++++++++++++++++++
spaced | Bin 0 -> 158712 bytes
spaced.cc | 248 +++++++++++++
variance.cpp | 788 ++++++++++++++++++++++++++++++++++++++++++
variance.h | 112 ++++++
12 files changed, 3824 insertions(+)
diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000..94a9ed0
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ <program> Copyright (C) <year> <name of author>
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..0515a61
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,4 @@
+
+all:
+ g++ -O3 -std=c++11 spaced.cc variance.cpp patternset.cpp extkey.cpp -o spaced -fopenmp
+
diff --git a/README b/README
new file mode 100644
index 0000000..ad2c70a
--- /dev/null
+++ b/README
@@ -0,0 +1,37 @@
+to compile type:
+ make
+
+run with:
+ ./spaced [options] <sequence file>
+
+options:
+ Options:
+ -h: print this help and exit"
+ -o <file>: output file name (default: DMat)"
+ -k <integer>: pattern weight (default 14)"
+ -l <integer>: pattern don't care positions (default 15)"
+ -n <integer>: number of patterns (default: 5)"
+ -f <file>: use patterns in <file> instead of random patterns"
+ -t <integer>: numer of threads (default: 25 threads)"
+ -r: don't consider the reverse complement"
+ -d EU | JS | EV: change distance type to Euclidean, Jensen-Shannon, evolutionary distance (default: EV)"
+_______________________________________________________________________________________________________________
+Scientific publications using Spaced Words should cite:
+
+ C.-A. Leimeister, M. Boden, S. Horwege, S. Lindner, B. Morgenstern (2014)
+ Fast alignment-free sequence comparison using spaced-word frequencies
+ Bioinformatics, DOI: 10.1093/bioinformatics/btu177 (http://bioinformatics.oxfordjournals.org/content/early/2014/04/03/bioinformatics.btu177)
+
+
+ S. Horwege, S. Linder, M. Boden, K. Hatje, M. Kollmar, C.-A. Leimeister, B. Morgenstern (2014)
+ Spaced words and kmacs: fast alignment-free sequence comparison based on inexact word matches
+ Nuc. Acids Research 42, W7-W11 (http://nar.oxfordjournals.org/content/42/W1/W7.abstract)
+
+
+ B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ Estimating evolutionary distances between genomic sequences from spaced-word matches
+ Algorithms for Molecular Biology 10, 5
+_______________________________________________________________________________________________________________
+
+contact:
+chris.leimeister at stud.uni-goettingen.de
diff --git a/extkey.cpp b/extkey.cpp
new file mode 100644
index 0000000..5395aad
--- /dev/null
+++ b/extkey.cpp
@@ -0,0 +1,73 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * extended key object file
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#include "extkey.h"
+
+/*---Variables---------------------------------------------------------------*/
+
+
+/*===Main-Part===============================================================*/
+/*---Constructor-& Init------------------------------------------------------*/
+extkey::extkey(){
+ extkey(1,1);
+}
+
+extkey::extkey(int pos, double value){
+ this->pos = pos;
+ this->value = value;
+}
+
+
+/*---Functions---------------------------------------------------------------*/
+int extkey::GetPos(){
+ return pos;
+}
+
+double extkey::GetValue(){
+ return value;
+}
+
+void extkey::SetValue(double value){
+ this->value = value;
+}
+
+bool extkey::operator < (extkey const& p) const{
+ return (value < p.value);
+}
+
+bool extkey::operator <= (extkey const& p) const{
+ return (value <= p.value);
+}
+
+bool extkey::operator > (extkey const& p) const{
+ return (value > p.value);
+}
+
+bool extkey::operator >= (extkey const& p) const{
+ return (value >= p.value);
+}
+
+bool extkey::operator == (extkey const& p) const{
+ return (value == p.value);
+}
+
+bool extkey::operator || (extkey const& p) const{
+ return (value || p.value);
+}
+
+bool extkey::operator != (extkey const& p) const{
+ return (value != p.value);
+}
diff --git a/extkey.h b/extkey.h
new file mode 100644
index 0000000..463a560
--- /dev/null
+++ b/extkey.h
@@ -0,0 +1,41 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * extended key object header
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#ifndef EXTKEY_H_
+#define EXTKEY_H_
+
+class extkey{
+ public:
+ extkey();
+ extkey(int pos, double value);
+
+ int GetPos();
+ double GetValue();
+ void SetValue(double value);
+
+ bool operator < (extkey const& p) const;
+ bool operator <= (extkey const& p) const;
+ bool operator > (extkey const& p) const;
+ bool operator >= (extkey const& p) const;
+ bool operator == (extkey const& p) const;
+ bool operator != (extkey const& p) const;
+ bool operator || (extkey const& p) const;
+
+ private:
+ int pos;
+ double value;
+};
+#endif
diff --git a/patternset.cpp b/patternset.cpp
new file mode 100644
index 0000000..5568559
--- /dev/null
+++ b/patternset.cpp
@@ -0,0 +1,1050 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * patternset object file
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#include "patternset.h"
+
+
+/*---Variables---------------------------------------------------------------*/
+
+std::default_random_engine generator(std::random_device{}());
+
+/*===Main-Part===============================================================*/
+/*---Constructor-------------------------------------------------------------*/
+/**
+ * Default constructor, sets the default vaulues, pattern will be generated automatically.
+ */
+patternset::patternset() {
+ this->size = 10;
+ this->length = new int[2];
+ this->length[0] = 14;
+ this->length[1] = 14;
+ this->weight = 8;
+ patternset(NULL, size, length, weight);
+}
+
+
+/**
+ * File constructor, sets values only from files; resets automatically if there are problems.
+ */
+patternset::patternset(char* pattern_file) {
+ this->pattern_file = pattern_file;
+ this->length = new int[2];
+ this->length[0] = 14;
+ this->length[1] = 14;
+ patternset(pattern_file, 10, length, 8);
+}
+
+/**
+ * Short constructor, sets some default vaulues, just pattern dimension is set.
+ *
+ * @param size
+ * The amount of patterns; pattern number.
+ *
+ * @param length
+ * The pattern length for each pattern of the pattern set.
+ *
+ * @param weigth
+ * The weight (match positions; '1') for each pattern of the pattern set.
+ */
+patternset::patternset(int size, int *length, int weight) {
+ patternset(NULL, size, length, weight);
+}
+
+
+/**
+ * Long constructor, sets the values; resets automatically, if there are problems.
+ *
+ * @param pattern_file File, that may contains submitted pattern
+ *
+ * @param align_file File, that may contains an alignment file to estimate p, q, l_hom, li and lj
+ *
+ * @param size The amount of patterns; pattern number.
+ *
+ * @param length The pattern length for each pattern of the pattern set.
+ *
+ * @param weigth The weight (match positions; '1') for each pattern of the pattern set.
+ *
+ * @param l_hom In theory, the amount of homologous positions of two sequences in an multiple alignment.
+ *
+ * @param l1 In theory, the first(represents each sequence i) sequence length of two observed sequences.
+ *
+ * @param l2 In theory, the second(represents each sequence j) sequence length of two observed sequences.
+ *
+ * @param p The match probability ( = #matches / #l_hom)
+ *
+ * @param q The background probability for each nucleotide A,C,G,T
+ */
+patternset::patternset(char* pattern_file, int size, int *length, int weight) {
+ this->pattern_file = pattern_file;
+ this->size = size;
+ this->length = length;
+ this->weight = weight;
+ update = false;
+ improve = true;
+ silent = false;
+ randpatleng = false;
+ ReInitPattern();
+}
+
+/**
+ * Default destructor, deletes all vectors and matrices in the object.
+ */
+patternset::~patternset() {
+ delete[] length;
+ Clear();
+}
+
+/*---Init--------------------------------------------------------------------*/
+
+/**
+ * Creates for the submitted or default values a set of pattern and calculates the first variance.
+ * If possible estimates p, q, all sequences lengths and all combinations of homologous sequence positions
+ *
+ * Also reset Pattern if needed, not necessary to create new object
+ */
+void patternset::ReInitPattern() { /*====Main-TODO=====*/
+ uint64_t tmp;
+ int lgth, lgt;
+ Clear();
+
+ lgth = 0;
+
+ if (pattern_file != NULL) {
+ TestPattern();
+ }
+ if(string_pat.size()== 0){
+ VerifyConditions();
+ CreateRandomPattern();
+ }
+ else{
+ lengths.clear();
+ for(int i = 0; i < size; i++){
+ tmp = ToBit(string_pat[i]);
+ pattern_set.push_back(tmp);
+ lgt = (int)(log2(tmp))+1;
+ lengths.push_back(lgt);
+ }
+ length[0] = lengths[0];
+ length[1] = length[0];
+ for(int i = 1; i < size; i++){
+ if(length[0] > lengths[i]){
+ length[0] = lengths[i];
+ }
+ if(length[1] < lengths[i]){
+ length[1] = lengths[i];
+ }
+ }
+ }
+ for(int i = 0; i < size; i++){
+ lgth += lengths[i];
+ }
+ length_mean = lgth/size;
+
+ return;
+}
+/**
+ * Due to some cases, it is possible necessary to
+ * reset some pattern parameters.
+ *
+ * This is not a complete delete, so therefore it
+ * it is only a part of the destructor.
+ */
+void patternset::Clear(){
+ lengths.clear();
+ pattern_set.clear();
+ string_pat.clear();
+}
+
+
+/*---Functions---------------------------------------------------------------*/
+/**
+ * If there is an submitted pattern_file, it will check, if it is containing
+ * a possible pattern, and if this pattern is in the correct format.
+ * Pattern can be parsed in a different formats
+ */
+void patternset::TestPattern() {
+ std::ifstream patternfile;
+ std::vector<std::string> pattern_tmp;
+ std::string tmp;
+ char tokens[4] = { '.',' ',',',';' }; /*These tokens are allowed to seperate patterns*/
+ size_t f_size;
+ bool start;
+
+ start = false;
+
+ patternfile.open(pattern_file);
+ patternfile.seekg(0, std::ios::end);
+ f_size = (size_t)patternfile.tellg();
+ if (!patternfile) {
+ SecureMessage("file", -1);
+ }
+ else if (f_size == 0) { /*[FILE].eof() does not recognize empty files...*/
+ SecureMessage("empty", -1);
+ }
+ else {
+ patternfile.close(); /*..therefore it has also to be closed and opened -.-** */
+ patternfile.open(pattern_file);
+ if(!silent){
+ std::cout << "Reading pattern from submitted patternfile ...\n" << std::endl;
+ }
+ patternfile >> tmp;
+ while (!patternfile.eof()) {
+ if (!ValidatePatternsFormat(tmp)) {
+ patternfile >> tmp; /*Ignoring each incorrect pattern. It is easier to calculate with all the rests*/
+ }
+ else {
+ pattern_tmp = SplitString(tmp, tokens);
+ for (unsigned int i = 0; i < pattern_tmp.size(); i++) {
+ string_pat.push_back(pattern_tmp[i]); /*Each pattern needs to be saved in its own std::string for comparison*/
+ }
+ patternfile >> tmp;
+ }
+ }
+ if (size > 0) { /*For an empty set we do not habe to validate*/
+ start = ValidatePatternConditions();
+ }
+ if (!start) { /*Some conditions, like changed weight is too much much amount of work*/
+ SecureMessage("pattern", -1); /*Therefore we just use our default values*/
+ size = 10;
+ weight = 8;
+ length = new int[2];
+ length[0] = 14;
+ length[1] = 14;
+ string_pat.clear(); /*Do not forget to reset, or the pattern set will not replaces, just increased*/
+ }
+ else {
+ size = (int )string_pat.size();
+ length = PatternLength(string_pat);
+ weight = PatternWeight(string_pat[0]);
+ if(!silent){
+ for (int i = 0; i < size - 1; i++) {
+ std::cout << string_pat[i] << "\n";
+ }
+ std::cout << string_pat[size - 1] << std::endl;
+ std::cout << "\n... Done!\n" << std::endl;
+ }
+ }
+ }
+ patternfile.close();
+}
+
+
+/**
+ * It will Verify the pattern conditions: is the weight correct,
+ * not negative values, weight not above length and so on..
+ */
+void patternset::VerifyConditions(){
+ std::vector<int> leng_new;
+ int tmp, diff, leng_old;
+
+ if(size <= 0){
+ SecureMessage("size", -1);
+ size = 10;
+ update = true;
+ }
+ if(length[1] < length[0]){
+ tmp = length[1];
+ length[1] = length[0];
+ length[0] = tmp;
+ }
+ if(weight < 2){
+ SecureMessage("weight_pat", -1);
+ weight = 8;
+ update = true;
+ }
+
+ if(length[0] < 3 || length[1] < 3){
+ SecureMessage("length",-1);
+ diff = length[1] - length[0];
+ length[0] = weight*2;
+ length[1] = length[0]+diff;
+ update = true;
+ }
+
+ diff = length[1]-length[0];
+ leng_old = length[0];
+ if(length[1] > 63){ /*max length of 64-Bit Integer*/
+ SecureMessage("length",-1);
+ if(diff >= 63-weight){
+ diff = 63-weight-1;
+ }
+ length[1] = 63;
+ length[0] = length[1]-diff;
+ update = true;
+ }
+ if(length[0] < weight){
+ SecureMessage("length",-1);
+ if(length[0] != leng_old){
+ length[0] = leng_old;
+ }
+ else{
+ length[0] = weight+2;
+ }
+ if(length[0] > length[1]){
+ length[0] = length[1];
+ }
+ update = true;
+ }
+
+ if(length[1] < weight || length[0] < weight || length[0] > 63 || length[1] > 63){
+ SecureMessage("nkl", -1); /*...if everything is fucked up and rescue does not work...*/
+ length[0] = 14;
+ length[1] = 14;
+ weight = 7;
+ update = true;
+ }
+
+ if(length[0] == weight && length[1] == weight){
+ if(size > 1){
+ SecureMessage("max_number_pattern",-1);
+ size = 1;
+ update = true;
+ }
+ if(improve){
+ improve = false;
+ SecureMessage("noimprove",-1);
+ update = true;
+ }
+ }
+
+ CreateLengths();
+
+ tmp = 1;
+ diff = 0;
+ if(lengths[0] == weight && size > 1){
+ diff++;
+ if(lengths[1] == lengths[0]){
+ lengths[1]++;
+ }
+ }
+
+ for(int i = 1+diff; i < size; i++){
+ if(lengths[i] < lengths[i-1]){
+ lengths[i] = lengths[i-1];
+ }
+ if(lengths[i] == lengths[i-1]){
+ tmp++;
+ if(tmp > MaxNumberPattern(weight-2, lengths[i]-2)){
+ lengths[i]++;
+ if(lengths[i] > length[1]){
+ lengths[i] = length[1];
+ size = i;
+ if(size <= 0){
+ size = 1;
+ }
+ }
+ tmp = 1;
+ }
+ }
+ else{
+ tmp = 1;
+ }
+ }
+
+
+ if(size != (int)lengths.size()){
+ SecureMessage("max_number_pattern",-1);
+ for(int i = 0; i < size; i++){
+ leng_new.push_back(lengths[i]);
+ }
+ lengths.clear();
+ lengths = leng_new;
+ update=true;
+ }
+ return;
+}
+
+
+
+/*---Functions---------------------------------------------------------------*/
+/*---Func-Create-------------------------------------------------------------*/
+/**
+ * Splits a string, read by a pattern file. Mayby each pattern does not get
+ * a new line, it has to be parsed, when a pattern starts and ends
+ *
+ * @param pattern_split The string containing a few pattern
+ *
+ * @param tokens The allowed tokens which can be used to seperate patterns in a line
+ */
+std::vector<std::string> patternset::SplitString(std::string pattern_split, char* tokens) {
+ std::vector<std::string> patternset;
+ std::string tmp = "";
+ bool flag_token = false;
+
+ for (unsigned int i = 0; i < pattern_split.length(); i++) {
+ for (unsigned int j = 0; j < strlen(tokens); j++) {
+ if (pattern_split[i] == tokens[j]) {
+ flag_token = true;
+ }
+ }
+ if (flag_token) { /*token found, which means in one line more patterns --> start new pattern, save last pattern*/
+ patternset.push_back(tmp);
+ tmp = "";
+ }
+ else {
+ tmp = tmp + pattern_split[i]; /*concatenating patternparts*/
+ }
+ flag_token = false;
+ }
+ patternset.push_back(tmp);
+ return patternset;
+}
+
+/**
+ * Validates a pattern, if it contains only pattern symbols and seperation tokens
+ *
+ * @param pattern_form The pattern which has to be investigate for symbols and tokens
+ *
+ * @return Returns true if this one pattern is in right format, false else
+ */
+bool patternset::ValidatePatternsFormat(std::string pattern_form) {
+ bool flag = true;
+
+ for (unsigned int i = 0; i < pattern_form.length(); i++) { /*allowed tokens in patternformat, also separating tokens*/
+ if (pattern_form[i] != '1' && pattern_form[i] != '0' && pattern_form[i] != ' ' && pattern_form[i] != ',' && pattern_form[i] != '.' && pattern_form[i] != ';') {
+ flag = false;
+ SecureMessage("format", -1);
+ break;
+ }
+ }
+ if (pattern_form[0] != '1' || pattern_form[pattern_form.length() - 1] != '1') {
+ flag = false;
+ SecureMessage("startend", -1);
+ }
+
+ return flag;
+}
+
+/**
+ * Validates a pattern set, if all patterns have the same length and weight
+ *
+ * @return Returns true if this one pattern is in right format, false else
+ */
+bool patternset::ValidatePatternConditions() {
+ int com_weight, com_size;
+ bool condition;
+
+ com_size = (int) string_pat.size();
+ com_weight = PatternWeight(string_pat[0]); /*as fixpoint saving first patternweight*/
+ condition = true;
+
+ for (int i = 1; i < com_size; i++) {
+ if (PatternWeight(string_pat[i]) != com_weight) {
+ SecureMessage("weight", i);
+ condition = false;
+ }
+ }
+ return condition;
+}
+
+/**
+ * Estimates the weight of a pattern
+ *
+ * @param pattern_str The pattern which has to be investigate for the weight
+ *
+ * @return Returns true if this one pattern is in right format, false else
+ */
+int patternset::PatternWeight(std::string pattern_wght) {
+ int str_weight;
+ int str_length;
+
+ str_weight = 0;
+ str_length = (int) pattern_wght.length();
+
+ for (int i = 0; i < str_length; i++) {
+ if (pattern_wght[i] == '1') {
+ str_weight++;
+ }
+ }
+ return str_weight;
+}
+
+/**
+ * Estimates the weight of a pattern
+ *
+ * @param pattern_str The pattern which has to be investigate for the weight
+ *
+ * @return Returns true if this one pattern is in right format, false else
+ */
+int* patternset::PatternLength(std::vector<std::string> pattern_length) {
+ int* pat_length;
+ int minp, maxp;
+
+ pat_length = new int[2];
+ minp = (int) pattern_length[0].length();
+ maxp = (int) pattern_length[0].length();
+
+ for (unsigned int i = 1; i < pattern_length.size(); i++) {
+ if ((int)pattern_length[i].length() > maxp) {
+ maxp = (int) pattern_length[i].length();
+ }
+ if ((int)pattern_length[i].length() < minp) {
+ minp = (int) pattern_length[i].length();
+ }
+ }
+ pat_length[0] = minp;
+ pat_length[1] = maxp;
+ return pat_length;
+}
+
+/**
+ * This function will create for our lenghts intervall
+ * for each pattern a specific length!.
+ *
+ */
+void patternset::CreateLengths(){
+ double step, leng2;
+ int leng, one;
+
+ one = 0;
+ leng = length[1];
+ if(length[0] != length[1]){
+ leng--;
+ }
+ if(length[0] == weight){
+ one = 1;
+ }
+ std::uniform_int_distribution<int> pat_leng(length[0]+one, leng);
+
+ if(randpatleng){
+ for (int i = 0; i < size-1; i++) {
+ if (i < (size / 4)) {
+ lengths.push_back(length[1]);
+ }
+ else {
+ lengths.push_back(pat_leng(generator));
+ }
+ }
+ lengths.push_back(length[0]);
+ }
+ else{
+ step = ((double)(length[1]-length[0]+1)/(size-1));
+ leng = length[0]+one;
+ leng2 = (double) leng;
+ if(step == 0){
+ step++;
+ }
+ for(int i = one; i < size; i++){
+ if(leng > length[1]){
+ leng = length[1];
+ }
+ lengths.push_back(leng);
+ leng2 += step;
+ leng = (int)leng2;
+ }
+ if(one==1){ /*if one pattern has to be without don't care positions*/
+ lengths.push_back(length[0]);
+ }
+ }
+ std::sort(lengths.begin(),lengths.end());
+}
+
+/**
+ * Creates random a set of pattern. For convention a pattern has to start and end with '1'
+ * Reason 10010 ~ 1001
+ *
+ * @return Returns a randomly created pattern set
+ */
+void patternset::CreateRandomPattern() {
+ std::vector<uint64_t> tmp;
+ uint64_t prototype, position, one;
+ int counter;
+
+ bool flag;
+
+ one = 1;
+
+ std::sort(lengths.rbegin(),lengths.rend());
+
+ for (int i = 0; i < size; i++) {
+ flag = true;
+ while(flag){
+ flag = false;
+ prototype = 1;
+ prototype = ((prototype << (lengths[i] - 1))| 1); /*Prototype, first and last position is a match!*/
+ std::uniform_int_distribution<int> distribution(1, lengths[i] - 2); /*better random generater than time, uses likelihood for evenly random distribution*/
+ counter = 2;
+ while (counter < weight) {
+ position = (uint64_t) distribution(generator);
+ if ((prototype & (one << position)) == 0) { /*have fun and fill with '1' randomly :D */
+ counter++;
+ prototype |= (one << position);
+ }
+ }
+ for(int j = 0; j < i; j++){
+ if(prototype == tmp[j]){
+ flag = true;
+ }
+ }
+ if(!flag){
+ tmp.push_back(prototype);
+ }
+ }
+ }
+ for(int i = 0; i < size; i++){
+ pattern_set.push_back(tmp[i]);
+ string_pat.push_back(BitString(tmp[i]));
+ }
+ tmp.clear();
+}
+
+
+/**
+ * Scans if there is another pattern in the same format
+ *
+ * @return returns boolean if there is another same pattern
+ */
+bool patternset::UniqPattern(int number) {
+ bool uniq = true;
+ for (int i = 0; i < size; i++) {
+ if (number != i) {
+ if (pattern_set[i] == pattern_set[number]) {
+ uniq = false;
+ }
+ }
+ }
+ return uniq;
+}
+
+/*--Func-Change--------------------------------------------------------------*/
+/**
+ * Changes two different positions ('1' and '0') in a specific bit-pattern
+ * Start and end are excluded
+ *
+ * @param number
+ * The pattern which has to be modified
+ */
+void patternset::ChangeBits(int number) {
+ int pos1, pos0;
+ bool flag = true;
+
+ if(lengths[number] == weight){ /*only possible, if lenghts[size-1] == weight*/
+ number--;
+ }
+
+ while (flag && improve) {
+ flag = false;
+ pos1 = SymbolRandPos(number, '1');
+ pos0 = SymbolRandPos(number, '0');
+ pattern_set[number] = ChangeBitPos(number, pos1, pos0);
+ if (!UniqPattern(number)) {
+ flag = true;
+ }
+ }
+}
+
+/**
+ * Exchange of two positions; changing match to don't care and vice versa.
+ *
+ * @param pos The pattern index of the modifying pattern
+ *
+ * @param pos_one Matchposition which has to be a don't care
+ *
+ * @param pos_zero Don'tCare position which has to be a match
+ *
+ * @return The modified pattern as an 64-Bit integer
+ */
+uint64_t patternset::ChangeBitPos(int pos, int pos_one, int pos_zero) {
+ uint64_t diff, one_mask, mask, p;
+
+ p = pattern_set[pos];
+ diff = 1;
+ diff = diff << pos_one;
+ one_mask = 1;
+
+ one_mask = one_mask << ((int)log2(p)+1);
+ one_mask--;
+
+ mask = one_mask ^ diff;
+ p = p & mask;
+
+ diff = 1;
+ diff = diff << pos_zero;
+ p = p | diff;
+
+ return p;
+}
+
+/**
+ * Returns for a specific symbol ('1', '0') a random chosen position
+ *
+ * @param number The pattern index of the modifying pattern
+ *
+ * @param symb The specific symbol, either '1' or '0'
+ *
+ * @return A random chosen position in the pattern of the specific symbol
+ */
+int patternset::SymbolRandPos(int number, char symb) {
+ std::vector<int> positions;
+ int pos;
+
+ positions = GetSymbol(number, symb);
+
+ std::uniform_int_distribution<int> distribution(0, (int) positions.size() - 1);
+
+ pos = distribution(generator);
+ pos = positions[pos];
+
+ positions.clear();
+ return pos;
+}
+
+/**
+ * Investigates for a specific symbol all positions of a pattern.
+ * For a match, the first and last positions are excluded!
+ *
+ * @param number The pattern index of the modifying pattern
+ *
+ * @param symb The specific symbol, either '1' or '0'
+ *
+ * @return All positions of the pattern of the specific symbol
+ */
+std::vector<int> patternset::GetSymbol(int number, char symb){
+ std::vector<int> positions;
+ int lgth;
+ uint64_t tmp, chr;
+
+ tmp = pattern_set[number];
+ lgth = (int)(log2(tmp));
+
+ chr = 0;
+ if (symb == '1') {
+ chr = 1;
+ }
+
+ for (int i = 1; i < lgth; i++) {
+ if (((tmp >> i) & 1) == chr) {
+ positions.push_back(i);
+ }
+ }
+ return positions;
+}
+
+
+/*---stuff-------------------------------------------------------------------*/
+/**
+* Method to print the current pattern
+*/
+void patternset::Print() {
+ ToString();
+ for (int i = 0; i < size; i++) {
+ std::cout << string_pat[i] << std::endl;;
+ }
+}
+
+/**
+ * Saves the current 64-Bit integer in the string pat as printable strings
+ */
+void patternset::ToString() {
+ for (int i = 0; i < size; i++) {
+ string_pat[i] = BitString(pattern_set[i]);
+ }
+}
+
+/**
+ * Turns a string into a 64-Bit Integer as binary pattern
+ *
+ * @param p The number of the pattern in the string pattern set
+ *
+ * @return The binary representation of the input string
+ */
+uint64_t patternset::ToBit(std::string p) {
+ std::string tmp;
+ uint64_t val, d;
+ val = 0;
+ d = 1;
+
+ for (unsigned int i = 0; i < p.length(); i++) {
+ if (p[i] == '1') {
+ val += d << (p.length() - (i + 1));
+ }
+ }
+ return val;
+}
+
+/**
+ * Changes a pattern in 64-Bit integer format into a printable string
+ *
+ * @param p The pattern in 64-bit integer format
+ *
+ * @return The printable string of the input pattern
+ */
+std::string patternset::BitString(uint64_t p) {
+ std::string bit_string;
+ int lgth;
+ lgth = (int) log2(p);
+
+ bit_string = "";
+
+ for (int i = lgth; i >= 0; i--) {
+ if (((p >> i) & 1) == 1) {
+ bit_string += "1";
+ }
+ else {
+ bit_string += "0";
+ }
+ }
+ bit_string += "\0";
+ return bit_string;
+}
+
+
+/**
+ * Determines the maximum number of patterns with weight and length
+ *
+ *@param weight Complete pattern weight-2, start and end have to be match and do not change
+ *
+ *@param length Complete pattern lengt-2, start and end have to be match and do not change
+ *
+ *@return max number of possible pattern
+ */
+double patternset::MaxNumberPattern(int p_weight, int p_length) {
+ double tmpa = Faculty(p_length);
+ double tmpb = Faculty(p_length - p_weight);
+ double tmpc = Faculty(p_weight);
+ double tmp = tmpa / (tmpb*tmpc);
+ return tmp;
+}
+
+/**
+ * Calculates the faculty
+ *
+ *@param value Calculating the faculty for value
+ *
+ *@return faculty(value) = 'value!'
+ */
+double patternset::Faculty(int value) {
+ double tmp;
+ tmp = 1;
+ for (int i = 1; i <= value; i++) {
+ tmp *= i;
+ }
+ return tmp;
+}
+
+/**
+ * A Method to collect all errormessages. Just easier for programmer to change
+ * the text or extend.
+ *
+ * @param errmsg Due to a few possible errormessages, this is the option, which has to be printed.
+ *
+ * @param pos The position of the incorrect patterns.
+ */
+void patternset::SecureMessage(std::string errmsg, int pos) {
+ /*if (errmsg == "wrongindex") {
+ std::cerr << "ERROR! Pattern " << pos << " does not exist... do nothing\n" << std::endl;
+ return;
+ }
+ if (errmsg == "file") {
+ std::cerr << "ERROR! Pattern file \'" << pattern_file << "\' could not be found!" << std::endl;
+ std::cerr << "Return to submitted or default values\n" << std::endl;
+ return;
+ }
+ if (errmsg == "empty") {
+ std::cerr << "ERROR! File \'" << pattern_file << "\' is an empty file!" << std::endl;
+ std::cerr << "Return to submitted or default values\n" << std::endl;
+ return;
+ }
+ if (errmsg == "pattern") {
+ std::cerr << "ERROR! Patternconditions from pattern file were not correct (different weight or length)!" << std::endl;
+ std::cerr << "Return to submitted or default values\n" << std::endl;
+ return;
+ }
+ if (errmsg == "nkl") {
+ std::cerr << "ERROR! Wrong values for weight, pattern number or pattern length!" << std::endl;
+ std::cerr << "Return to default values\n" << std::endl;
+ return;
+ }
+ if (errmsg == "weight_pat") {
+ std::cerr << "ERROR! Weight of a pattern cannot under 2!" << std::endl;
+ std::cerr << "Return to submitted or default values\n" << std::endl;
+ return;
+ }
+ if (errmsg == "format") {
+ std::cerr << "FORMAT-ERROR: Pattern containes illegal characters!" << std::endl;
+ std::cerr << "Allowed characters: '0','1' for pattern; ','|'.'|';'|' ' to seperate patterns" << std::endl;
+ std::cerr << "Go on to next pattern.\n" << std::endl;
+ return;
+ }
+ if (errmsg == "max_number_pattern") {
+ std::cerr << "The number of patterns is too high for your configuration and will be reset!" << std::endl;
+ return;
+ }
+ if (errmsg == "size") {
+ std::cerr << "ERROR! The number of patterns has to be greater than 0!\n" << std::endl;
+ return;
+ }
+ if (errmsg == "weight") {
+ std::cerr << "ERROR! By comparing with the first pattern, the " << pos + 1 << ". pattern has a different weight!\n" << std::endl;
+ return;
+ }
+ if(errmsg == "length"){
+ std::cerr << "ERROR! The patternlengths has to be greater or equal the weight, greater than 2 and below 63!" << std::endl;
+ return;
+ }
+ if (errmsg == "startend") {
+ std::cerr << "FORMAT-ERROR: Pattern has to start and end with a match position '1' !\n" << std::endl;
+ return;
+ }
+ if (errmsg == "noimprove") {
+ std::cerr << "Using your pattern conditions it is not sensible to improve your pattern, sorry!" << std::endl;
+ std::cerr << "Deactivating improve mode\n" << std::endl;
+ return;
+ }
+ if (errmsg == "update"){
+ std::cerr << "\n####IMPROTANT####\nDue to some configuration errors, your submitted parameters have been updated!" << std::endl;
+ return;
+ }*/
+}
+
+
+/*---Set-&-GetFunc-----------------------------------------------------------*/
+/**
+ * Changes to silent output, nothing will be printed automatically!
+ * Errors are going to be printed!
+ */
+void patternset::Silent(){
+ silent = true;
+}
+
+/**
+ * Changes the patternlengths from evenly distributed to random pattern lengths.
+ * One quarter has maximum length, the rest randomly distributed between
+ * maximum and minimum length.
+ */
+void patternset::RandPatLength(){
+ if(!randpatleng){
+ randpatleng = true;
+ ReInitPattern();
+ }
+}
+
+/**
+ * Returns the weight of each pattern, the match positions
+ *
+ * @return returns weight
+ */
+int patternset::GetWeight() {
+ return weight;
+}
+
+/**
+ * Returns the amount of patters; number of patterns
+ *
+ * @return returns size
+ */
+int patternset::GetSize() {
+ return size;
+}
+
+/**
+ * Returns the length of each pattern
+ *
+ * @return returns length
+ */
+int* patternset::GetLength() {
+ return length;
+}
+
+/**
+ * Returns the mean of all pattern lengths
+ *
+ * @return The mean of all lengths;
+ */
+int patternset::GetLengthMean(){
+ return length_mean;
+}
+
+/**
+ * Returns a boolean, if an improvement is possible
+ * Imagine, a improvment for pattern '101' is not possible!
+ *
+ * @return The boolean, if this improvement is possible
+ */
+bool patternset::GetImprove() {
+ return improve;
+}
+
+/**
+ * Returns a boolean, if due to some configuration errors
+ * values set by the user had to be updated to better configuration.
+ * Imagine, the weight was above the pattern length.
+ *
+ * @return The boolean, if something had been updated.
+ */
+bool patternset::GetUpdate() {
+ return update;
+}
+
+
+/**
+ * Returns a specific pattern of the pattern set in 64-Bit integer format
+ *
+ * @return specific string in 64-Bit integer format
+ */
+uint64_t patternset::GetPattern(int number) {
+ if (number >= size) {
+ SecureMessage("wrongindex", number);
+ return 0;
+ }
+ else {
+ return pattern_set[number];
+ }
+}
+
+/**
+ * Returns complete pattern set
+ *
+ * @return complete pattern set
+ */
+std::vector<std::string> patternset::GetStringPattern() {
+ ToString();
+ return string_pat;
+}
+
+/**
+ * Returns only one pattern as a string from the current patternset.
+ *
+ * @param number The index of the chosen pattern
+ *
+ * @return The pattern as a printable string
+ */
+std::string patternset::GetString(int number){
+ if (number >= size) {
+ SecureMessage("wrongindex", number);
+ return NULL;
+ }
+ else {
+ ToString();
+ return string_pat[number];
+ }
+}
+
+/**
+ * Changes a pattern in 64-Bit integer format
+ *
+ * @param number The index of the pattern which will be changed
+ *
+ * @param patt The new pattern
+ */
+void patternset::SetPattern(int number, uint64_t patt){
+ if (number >= size) {
+ SecureMessage("wrongindex", number);
+ }
+ else {
+ pattern_set[number] = patt;
+ }
+}
diff --git a/patternset.h b/patternset.h
new file mode 100644
index 0000000..12a60b2
--- /dev/null
+++ b/patternset.h
@@ -0,0 +1,103 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * patternset object header
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#ifndef PATTERNSET_H_
+#define PATTERNSET_H_
+
+
+#include <iostream>
+#include <fstream>
+#include <random>
+#include <vector>
+#include <string>
+#include <algorithm>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <inttypes.h>
+#include "extkey.h"
+
+class patternset {
+public:
+ patternset();
+ patternset(char* pattern_file);
+ patternset(int size, int *length, int weight);
+ patternset(char* pattern_file, int size, int *length, int weight);
+ ~patternset();
+
+ void ReInitPattern();
+
+ void ChangeBits(int number);
+ uint64_t ChangeBitPos(int pos, int pos_one, int pos_zero);
+ bool UniqPattern(int number);
+
+ std::vector<std::string> GetStringPattern();
+ std::string GetString(int number);
+ uint64_t GetPattern(int number);
+ int GetWeight();
+ int GetSize();
+ int* GetLength();
+ int GetLengthMean();
+ bool GetImprove();
+ bool GetUpdate();
+
+ double GetValue(int number);
+ void SetPattern(int number, uint64_t patt);
+ void Print();
+ void Silent();
+ void RandPatLength();
+
+protected:
+ void Clear();
+ void TestPattern();
+ void VerifyConditions();
+
+ std::vector<std::string> SplitString(std::string pattern, char* tokens);
+ bool ValidatePatternsFormat(std::string pattern_form);
+ bool ValidatePatternConditions();
+ int PatternWeight(std::string pattern_wght);
+ int* PatternLength(std::vector<std::string> pattern_length);
+
+ void CreateLengths();
+ void CreateRandomPattern();
+ bool IsSetScore(int pattern);
+ int SymbolRandPos(int number, char symb);
+ int SymbolCalcPos(int number, char symb);
+ std::vector<int> GetSymbol(int number, char symb);
+
+ void ToString();
+ uint64_t ToBit(std::string p);
+ std::string BitString(uint64_t p);
+
+ double MaxNumberPattern(int p_weight, int p_length);
+ double Faculty(int value);
+ void SecureMessage(std::string errmsg, int pos);
+
+private:
+ std::vector<uint64_t> pattern_set;
+ std::vector<std::string> string_pat;
+ std::vector<int> lengths;
+ int size;
+ int *length;
+ int length_mean;
+ int weight;
+ char* pattern_file = NULL;
+ bool improve;
+ bool update;
+ bool silent;
+ bool randpatleng;
+};
+#endif
diff --git a/sort.h b/sort.h
new file mode 100644
index 0000000..bef900e
--- /dev/null
+++ b/sort.h
@@ -0,0 +1,694 @@
+/**
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details at
+ * http://www.gnu.org/copyleft/gpl.html
+ *
+ */
+
+#ifndef SORT_H_
+#define SORT_H_
+
+#include <stdio.h>
+ #include <string.h>
+#include <getopt.h>
+#include <string>
+#include <stdlib.h>
+#include <iostream>
+#include <ctype.h>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <fstream>
+#include <tr1/unordered_map>
+#include <map>
+#include <limits>
+#include <unordered_map>
+#include <unistd.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+#include <iomanip>
+#include <random>
+#include <bitset>
+
+#define DNA_BUCKETS 256
+#define PROT_BUCKETS 26
+#define EU 0
+#define JS 1
+#define EV 2
+
+using namespace std;
+
+template<typename uint>
+struct word{
+ unsigned long long key;
+ unsigned short seq;
+ bool operator<( const word& val ) const {
+ return key < val.key;
+ }
+};
+
+template<typename uint>
+struct sequence{
+ uint headerStart;
+ uint headerEnd;
+ uint seqStart;
+ uint seqEnd;
+ uint seqWordEnd;
+ uint length=0;
+ uint revCompStart;
+ uint revCompEnd;
+ uint revCompWordEnd;
+ vector<double> frequencies={0,0,0,0};
+ vector<double> frequenciesRevComp={0,0,0,0};
+};
+
+void getFirstBits(vector<unsigned char>&pattern, unsigned char* firstBits, unsigned long long pos, unsigned long length, unsigned char* seqData, int weight){
+ unsigned char alphabet[256]={};
+ alphabet['A']=0;alphabet['C']=1;alphabet['G']=2;alphabet['T']=3; // mapping A -> 00, C -> 01 etc. without overhead
+ unsigned char w;
+ // mapping like:
+ // ACAG -> w = 00 01 00 10
+ // uses one byte for first bits, thus max weight = 4
+ unsigned char p0,p1,p2,p3;
+ p0=pattern[0];p1=pattern[1];p2=pattern[2];p3=pattern[3];
+ if(weight>=4){
+ for(int i=0; i<length;i++){
+ w=0;
+ unsigned char* tmp=&seqData[i+pos];
+ w|=alphabet[*(tmp+p0)] << 6;
+ w|=alphabet[*(tmp+p1)] << 4;
+ w|=alphabet[*(tmp+p2)] << 2;
+ w|=alphabet[*(tmp+p3)];
+ firstBits[i]=w;
+ }
+ }
+ else if (weight==3){
+ unsigned char p0,p1,p2;
+ p0=pattern[0];p1=pattern[1];p2=pattern[2];
+ for(int i=0; i<length;i++){
+ w=0;
+ unsigned char* tmp=&seqData[i+pos];
+ w|=alphabet[*(tmp+p0)] << 4;
+ w|=alphabet[*(tmp+p1)] << 2;
+ w|=alphabet[*(tmp+p2)];
+ firstBits[i]=w;
+ }
+ }
+ else if(weight==2){
+ unsigned char p0,p1;
+ p0=pattern[0];p1=pattern[1];
+ for(int i=0; i<length;i++){
+ w=0;
+ unsigned char* tmp=&seqData[i+pos];
+ w|=alphabet[*(tmp+p0)] << 2;
+ w|=alphabet[*(tmp+p1)];
+ firstBits[i]=w;
+ }
+
+ }
+ else{
+ unsigned char p0;
+ p0=pattern[0];
+ for(int i=0; i<length;i++){
+ unsigned char* tmp=&seqData[i+pos];
+ w=alphabet[*(tmp+p0)];
+ firstBits[i]=w;
+ }
+ }
+}
+
+template<typename uint>
+void writeDmat(vector<vector<double> > dmat, vector<sequence<uint> >& sequences, string filename, unsigned char* seqData){
+ ofstream outfile;
+ outfile.open(filename.c_str());
+ outfile << sequences.size() << endl;
+ for (int i = 0; i < sequences.size(); i++) {
+ for(uint k=0; k<10;k++){
+ if( k+sequences[i].headerStart <= sequences[i].headerEnd )
+ outfile << seqData[k+sequences[i].headerStart];
+ else
+ outfile << " ";
+ }
+ outfile << " ";
+ for (int j = 0; j < sequences.size(); j++) {
+ if (i > j)
+ outfile << setprecision(12) << dmat[i][j] << " ";
+ else if(j>i)
+ outfile << setprecision(12) << dmat[j][i] << " ";
+ else
+ outfile << setprecision(12) << "0" << " ";
+ }
+ outfile << endl;
+ }
+ //outfile << dmat[1][0] << endl;
+}
+
+template<typename uint>
+void spacedDNA(vector<string>& patternSet, unsigned char* seqData, uint n, int distance, int threads, int weight, int dontCare, bool revComp, string output){
+ #ifdef _OPENMP
+ omp_set_dynamic(0);
+ int ell=dontCare+weight-1;
+ omp_set_num_threads(threads);
+ #endif
+ unsigned char alphabet[256]={};
+ alphabet['A']=0;alphabet['C']=1;alphabet['G']=2;alphabet['T']=3;
+ vector<sequence<uint> > sequences;
+ int seqNum;
+ uint totalSeqLength=0;
+ for(uint i=0,j=0; i<n;i++){
+ if(seqData[i]==';'){ //remove comments
+ while(seqData[++i]!='\n');
+ }else if(seqData[i]=='>'){//read header
+ sequence<uint> tmpSeq={};
+ seqData[j++]=seqData[i];
+ tmpSeq.headerStart=j;
+ totalSeqLength++;
+ if(sequences.size()>0){
+ sequences.back().seqEnd=j-1;
+ sequences.back().seqWordEnd=j-(weight+dontCare);
+ }
+ while(seqData[++i]!='\n'){
+ seqData[j++]=seqData[i];
+ totalSeqLength++;
+ }
+ tmpSeq.headerEnd=j-1;
+ tmpSeq.seqStart=j;
+ sequences.push_back(tmpSeq);
+ }else if(isalpha(seqData[i])){
+ char c=toupper(seqData[i]);
+ switch(c){
+ case 'A':
+ seqData[j++]='A';
+ sequences.back().frequencies[0]++;
+ sequences.back().length++;
+ break;
+ case 'C':
+ seqData[j++]='C';
+ sequences.back().frequencies[1]++;
+ sequences.back().length++;
+ break;
+ case 'G':
+ seqData[j++]='G';
+ sequences.back().frequencies[2]++;
+ sequences.back().length++;
+ break;
+ case 'T':
+ seqData[j++]='T';
+ sequences.back().frequencies[3]++;
+ sequences.back().length++;
+ break;
+ default:
+ seqData[j++]='N';
+ sequences.back().length++;
+ break;
+ }
+ totalSeqLength++;
+ }
+ }
+ sequences.back().seqEnd=totalSeqLength;
+ sequences.back().seqWordEnd=totalSeqLength-(weight+dontCare)+1;
+ seqNum=sequences.size();
+ uint cnt=0;
+ for(uint i=0; i<seqNum;i++){
+ uint start=sequences[i].seqStart;
+ uint end=sequences[i].seqEnd-1;
+ sequences[i].revCompStart=cnt+totalSeqLength;
+ sequences[i].revCompEnd=cnt+totalSeqLength+sequences[i].length;
+ sequences[i].revCompWordEnd=sequences[i].revCompEnd-(weight+dontCare)+1;
+ sequences[i].frequenciesRevComp[0]=sequences[i].frequencies[3];
+ sequences[i].frequenciesRevComp[1]=sequences[i].frequencies[2];
+ sequences[i].frequenciesRevComp[2]=sequences[i].frequencies[1];
+ sequences[i].frequenciesRevComp[3]=sequences[i].frequencies[0];
+ for(uint j=end; j>=start ;j--){
+ switch(seqData[j]){
+ case 'A': seqData[cnt++ + totalSeqLength]='T';break;
+ case 'C': seqData[cnt++ + totalSeqLength]='G';break;
+ case 'G': seqData[cnt++ + totalSeqLength]='C';break;
+ case 'T': seqData[cnt++ + totalSeqLength]='A';break;
+ case 'N': seqData[cnt++ + totalSeqLength]='N';break;
+ default: exit(-1);
+ }
+ }
+ }
+ for(int i=0; i< sequences.size();i++){
+ for(int j=0; j<4;j++){
+ sequences[i].frequencies[j]=sequences[i].frequencies[j]/sequences[i].length;
+ sequences[i].frequenciesRevComp[j]=sequences[i].frequenciesRevComp[j]/sequences[i].length;
+ }
+ }
+ totalSeqLength*=2;
+ vector<unsigned char> firstBits(totalSeqLength);
+ vector<uint> prefixSum(256,0);
+ vector<unsigned char> matchPos;
+ vector<unsigned char> dCarePos;
+ vector<vector <vector <double> > > dmat(threads, vector< vector<double> >(seqNum, vector<double>(seqNum,0)));
+ vector<vector <double> > dmatFinal(seqNum, vector<double>(seqNum,0));
+ vector<vector<double> > row(threads, vector<double>(seqNum*2));
+ for(int p = 0; p < patternSet.size(); p++ ){
+ for(int k = 0; k < patternSet[p].length(); k++){
+ if (patternSet[p][k] == '1')
+ matchPos.push_back(k);
+ else
+ dCarePos.push_back(k);
+ }
+ getFirstBits(matchPos, &firstBits[0], 0, totalSeqLength, seqData, weight);
+ fill(prefixSum.begin(), prefixSum.end(), 0);
+
+ for(int i=0;i<seqNum;i++){
+ for(uint j=sequences[i].seqStart; j<sequences[i].seqWordEnd;j++){
+ ++prefixSum[(firstBits[j])];
+ }
+ for(uint j=sequences[i].revCompStart; j<sequences[i].revCompWordEnd;j++){
+ ++prefixSum[(firstBits[j])];
+ }
+ }
+ #pragma omp parallel for schedule(runtime)
+ for(unsigned int k=0; k<DNA_BUCKETS; k++){
+ vector<word<uint> > words(prefixSum[k]);
+ double cnt=0;
+ unsigned long long w;
+ unsigned char bits;
+ for(int i=0;i<seqNum;i++){
+ for(uint j=sequences[i].seqStart; j<sequences[i].seqWordEnd;j++){
+ if(firstBits[j]==k){
+ w=0;
+ bits=weight*2-2;
+ unsigned char c;
+ bool correctWord=true;
+ for(int o=0; o<matchPos.size();o++){
+ uint pos=matchPos[o];
+ c= seqData[j+pos];
+ if(c=='N')
+ correctWord=false;
+ w|=(unsigned long long) alphabet[c] << bits;
+ bits-=2;
+ }
+ if(correctWord){
+ words[cnt].key=w;
+ words[cnt++].seq=i;
+ }
+ }
+ }
+ if(revComp){
+ for(uint j=sequences[i].revCompStart; j<sequences[i].revCompWordEnd;j++){
+ if(firstBits[j]==k){
+ w=0;
+ bits=weight*2-2;
+ unsigned char c;
+ bool correctWord=true;
+ for(int o=0; o<matchPos.size();o++){
+ uint pos=matchPos[o];
+ c= seqData[j+pos];
+ if(c=='N')
+ correctWord=false;
+ w|=(unsigned long long) alphabet[c] << bits;
+ bits-=2;
+ }
+ if(correctWord){
+ words[cnt].key=w;
+ words[cnt++].seq=i+seqNum;
+ }
+ }
+ }
+ }
+ }
+ std::sort(words.begin(),words.begin()+cnt);
+ int thread=0;
+ #ifdef _OPENMP
+ thread=omp_get_thread_num();
+ #endif
+ if(distance==EU){
+ fill(row[thread].begin(), row[thread].end(), 0);
+ for(uint i=0; i<cnt;){
+ unsigned long long key_t=words[i].key;
+ row[thread][words[i].seq]++;
+ while(++i<cnt &&key_t==words[i].key){
+ row[thread][words[i].seq]++;
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ double t1;
+ if(revComp)
+ t1=abs(row[thread][i]-(row[thread][j]+row[thread][j+seqNum]));
+ else
+ t1=abs(row[thread][i]-row[thread][j]);
+ dmat[thread][i][j]+=t1*t1;
+ }
+ }
+ fill(row[thread].begin(), row[thread].end(), 0);
+ }
+ }
+ else if(distance==EV){
+ fill(row[thread].begin(), row[thread].end(), 0);
+ for(uint i=0; i<cnt;){
+ unsigned long long key_t=words[i].key;
+ row[thread][words[i].seq]++;
+ while(++i<cnt &&key_t==words[i].key){
+ row[thread][words[i].seq]++;
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ if(revComp)
+ dmat[thread][i][j]+=min(row[thread][i],row[thread][j]+row[thread][j+seqNum]);
+ else
+ dmat[thread][i][j]+=min(row[thread][i],row[thread][j]);
+ }
+ }
+ fill(row[thread].begin(), row[thread].end(), 0);
+ }
+ }
+ else if(distance==JS){
+ fill(row[thread].begin(), row[thread].end(), 0);
+ for(uint i=0; i<cnt;){
+ unsigned long long key_t=words[i].key;
+ row[omp_get_thread_num()][words[i].seq]++;
+ while(++i<cnt &&key_t==words[i].key){
+ row[thread][words[i].seq]++;
+ }
+ for(int i=0; i<seqNum; i++){
+ if(revComp){
+ row[thread][i+seqNum]=(row[thread][i]+row[thread][i+seqNum])/((sequences[i].seqWordEnd-sequences[i].seqStart)*2);
+ row[thread][i]/=sequences[i].seqWordEnd-sequences[i].seqStart;
+ }
+ else
+ row[thread][i]/=sequences[i].seqWordEnd-sequences[i].seqStart;
+ }
+ for(int i=0;i<seqNum;i++){
+ double row_i=row[thread][i];
+ for(int j=0; j<i; j++){
+ double m;
+ double row_j;
+ if(revComp)
+ row_j=row[thread][j+seqNum];
+ else
+ row_j=row[thread][j];
+ m=0.5*(row_i+row_j);
+ if(m==0)
+ continue;
+ if(row_i!=0 &&row_j!=0){
+ dmat[thread][i][j]+= log2(row_i/m ) * row_i;
+ dmat[thread][j][i]+= log2(row_j/m) * row_j;
+ }
+ else if(row_i!=0)
+ dmat[thread][i][j]+= log2(row_i/m) * row_i;
+ else if(row_j!=0)
+ dmat[thread][j][i]+= log2(row_j/m) * row_j;
+ }
+ }
+ fill(row[thread].begin(), row[thread].end(), 0);
+ }
+ }
+ }
+ matchPos.clear();
+ dCarePos.clear();
+ if(distance==EU || distance==JS){
+ for(int p=1; p<threads; p++){
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<seqNum; j++){
+ dmat[0][i][j]+=dmat[p][i][j];
+ dmat[p][i][j]=0;
+ }
+ }
+ }
+ if(distance==EU){
+ #ifdef _OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmat[0][i][j]=sqrt(dmat[0][i][j]);
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ dmatFinal[i][j]+=dmat[0][i][j];
+ dmat[0][i][j]=0;
+ }
+ }
+ }
+ if(distance==JS){
+ #ifdef _OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ dmat[0][i][j]=(0.5*dmat[0][i][j]+0.5*dmat[0][j][i]);
+ }
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ dmatFinal[i][j]+=dmat[0][i][j];
+ dmat[0][i][j]=0;
+ dmat[0][j][i]=0;
+ }
+ }
+ }
+ }
+ }
+ if(patternSet.size()>1&&distance!=EV){
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmatFinal[i][j]/=patternSet.size();
+ }
+ }
+ else if(distance==EV){
+ for(int p=1; p<threads; p++){
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<seqNum; j++){
+ dmat[0][i][j]+=dmat[p][i][j];
+ }
+ }
+ }
+ #ifdef _OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ double min=std::min(sequences[i].length,sequences[j].length)-ell;
+ double max=std::max(sequences[i].length,sequences[j].length)-ell;
+ double q=0;
+ for(int k=0; k<4;k++){
+ if(revComp)
+ q+=((sequences[i].frequencies[k]+sequences[i].frequenciesRevComp[k])*0.5)*((sequences[j].frequencies[k]+sequences[j].frequenciesRevComp[k])*0.5);
+ else
+ q+=sequences[i].frequencies[k]*sequences[j].frequencies[k];
+ }
+ double valueUnderRoot;
+ if(revComp)
+ valueUnderRoot=dmat[0][i][j]/(patternSet.size()*min)-2*max*pow(q,weight);
+ else
+ valueUnderRoot=dmat[0][i][j]/(patternSet.size()*min)-max*pow(q,weight);
+ if(valueUnderRoot>=0){
+ double p = pow(valueUnderRoot, (1.0 / weight));
+ dmat[0][i][j]=-0.75 * log((4.0/3.0)*p-(1.0/3.0));
+ }
+ else{
+ dmat[0][i][j]=1.2;
+ }
+ }
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmatFinal[i][j]+=dmat[0][i][j];
+ }
+ }
+ writeDmat<uint>(dmatFinal,sequences, output, seqData);
+}
+
+template<typename uint>
+void spacedProt(vector<string>& patternSet, unsigned char* seqData, uint n, int distance, int threads, int weight, int dontCare, string output){
+ #ifdef _OPENMP
+ omp_set_dynamic(0);
+ int ell=dontCare+weight-1;
+ omp_set_num_threads(threads);
+ #endif
+ unsigned char alphabet[256]={};
+ alphabet['A']=0;alphabet['B']=1;alphabet['C']=2;alphabet['D']=3;alphabet['E']=4;alphabet['F']=5;alphabet['G']=6;alphabet['H']=7;alphabet['I']=8;alphabet['J']=9;
+ alphabet['K']=10;alphabet['L']=11;alphabet['M']=12;alphabet['N']=13;alphabet['O']=14;alphabet['P']=15;alphabet['Q']=16;alphabet['R']=17;alphabet['S']=18;alphabet['T']=19;
+ alphabet['U']=20;alphabet['V']=21;alphabet['W']=22;alphabet['X']=23;alphabet['Y']=24;alphabet['Z']=25;
+ vector<sequence<uint> > sequences;
+ int seqNum;
+ uint totalSeqLength=0;
+ for(uint i=0,j=0; i<n;i++){
+ if(seqData[i]==';'){ //remove comments
+ while(seqData[++i]!='\n');
+ }else if(seqData[i]=='>'){//read header
+ sequence<uint> tmpSeq={};
+ seqData[j++]=seqData[i];
+ tmpSeq.headerStart=j;
+ totalSeqLength++;
+ if(sequences.size()>0){
+ sequences.back().seqEnd=j-1;
+ sequences.back().seqWordEnd=j-(weight+dontCare);
+ }
+ while(seqData[++i]!='\n'){
+ seqData[j++]=seqData[i];
+ totalSeqLength++;
+ }
+ tmpSeq.headerEnd=j-1;
+ tmpSeq.seqStart=j;
+ sequences.push_back(tmpSeq);
+ }else if(isalpha(seqData[i])){
+ char c=toupper(seqData[i]);
+ seqData[j++]=c;
+ sequences.back().length++;
+ totalSeqLength++;
+ }
+ }
+ sequences.back().seqEnd=totalSeqLength;
+ sequences.back().seqWordEnd=totalSeqLength-(weight+dontCare)+1;
+ seqNum=sequences.size();
+ vector<uint> prefixSum(PROT_BUCKETS,0);
+ vector<unsigned char> matchPos;
+ vector<unsigned char> dCarePos;
+ vector<vector <vector <double> > > dmat(threads, vector< vector<double> >(seqNum, vector<double>(seqNum)));
+ vector<vector <double> > dmatFinal(seqNum, vector<double>(seqNum,0));
+ vector<vector<double> > row(threads, vector<double>(seqNum));
+ for(int p = 0; p < patternSet.size(); p++ ){
+ for(int k = 0; k < patternSet[p].length(); k++){
+ if (patternSet[p][k] == '1')
+ matchPos.push_back(k);
+ else
+ dCarePos.push_back(k);
+ }
+ fill(prefixSum.begin(), prefixSum.end(), 0);
+ for(int i=0;i<seqNum;i++){
+ for(uint j=sequences[i].seqStart; j<sequences[i].seqWordEnd;j++){
+ ++prefixSum[alphabet[seqData[j]]];
+ }
+ }
+ #pragma omp parallel for schedule(runtime)
+ for(unsigned int k=0; k<PROT_BUCKETS; k++){
+ vector<word<uint> > words(prefixSum[k]);
+ double cnt=0;
+ unsigned long long w;
+ unsigned char bits;
+ for(int i=0;i<seqNum;i++){
+ for(uint j=sequences[i].seqStart; j<sequences[i].seqWordEnd;j++){
+ if(alphabet[seqData[j]]==k){
+ w=0;
+ bits=weight*5-5;
+ unsigned char c;
+ for(int o=0; o<matchPos.size();o++){
+ uint pos=matchPos[o];
+ c= seqData[j+pos];
+ w|=(unsigned long long) alphabet[c] << bits;
+ bits-=5;
+ }
+ words[cnt].key=w;
+ words[cnt++].seq=i;
+ }
+ }
+ }
+ std::sort(words.begin(),words.begin()+cnt);
+ int thread=0;
+ #ifdef _OPENMP
+ thread=omp_get_thread_num();
+ #endif
+ if(distance==EU){
+ fill(row[thread].begin(), row[thread].end(), 0);
+ for(uint i=0; i<cnt;){
+ unsigned long long key_t=words[i].key;
+ row[thread][words[i].seq]++;
+ while(++i<cnt &&key_t==words[i].key){
+ row[thread][words[i].seq]++;
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ double t1;
+ t1=abs(row[thread][i]-row[thread][j]);
+ dmat[thread][i][j]+=t1*t1;
+ }
+ }
+ fill(row[thread].begin(), row[thread].end(), 0);
+ }
+ }
+ else if(distance==JS){
+ fill(row[thread].begin(), row[thread].end(), 0);
+ for(uint i=0; i<cnt;){
+ unsigned long long key_t=words[i].key;
+ row[omp_get_thread_num()][words[i].seq]++;
+ while(++i<cnt &&key_t==words[i].key){
+ row[thread][words[i].seq]++;
+ }
+ for(int i=0; i<seqNum; i++){
+ row[thread][i]/=sequences[i].seqWordEnd-sequences[i].seqStart;
+ }
+ for(int i=0;i<seqNum;i++){
+ double row_i=row[thread][i];
+ for(int j=0; j<i; j++){
+ double m;
+ double row_j;
+ row_j=row[thread][j];
+ m=0.5*(row_i+row_j);
+ if(m==0)
+ continue;
+ if(row_i!=0 &&row_j!=0){
+ dmat[thread][i][j]+= log2(row_i/m ) * row_i;
+ dmat[thread][j][i]+= log2(row_j/m) * row_j;
+ }
+ else if(row_i!=0)
+ dmat[thread][i][j]+= log2(row_i/m) * row_i;
+ else if(row_j!=0)
+ dmat[thread][j][i]+= log2(row_j/m) * row_j;
+ }
+ }
+ fill(row[thread].begin(), row[thread].end(), 0);
+ }
+ }
+ }
+ matchPos.clear();
+ dCarePos.clear();
+ for(int p=1; p<threads; p++){
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<seqNum; j++){
+ dmat[0][i][j]+=dmat[p][i][j];
+ }
+ }
+ }
+ if(distance==EU){
+ #ifdef _OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmat[0][i][j]=sqrt(dmat[0][i][j]);
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmatFinal[i][j]+=dmat[0][i][j];
+ }
+ }
+ if(distance==JS){
+ #ifdef _OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++){
+ dmat[0][i][j]=(0.5*dmat[0][i][j]+0.5*dmat[0][j][i]);
+ }
+ }
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmatFinal[i][j]+=dmat[0][i][j];
+ }
+ }
+ }
+ if(patternSet.size()>1&&distance!=EV){
+ for(int i=0;i<seqNum;i++){
+ for(int j=0; j<i; j++)
+ dmatFinal[i][j]/=patternSet.size();
+ }
+ }
+ writeDmat<uint>(dmatFinal,sequences, output, seqData);
+}
+
+#endif
+
diff --git a/spaced b/spaced
new file mode 100755
index 0000000..55eaf90
Binary files /dev/null and b/spaced differ
diff --git a/spaced.cc b/spaced.cc
new file mode 100644
index 0000000..0d804e3
--- /dev/null
+++ b/spaced.cc
@@ -0,0 +1,248 @@
+/**
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details at
+ * http://www.gnu.org/copyleft/gpl.html
+ *
+ */
+
+#include "sort.h"
+#include "variance.h"
+
+bool revComp=true;
+int number = 5, dontcare = 15, weight = 14, threads = 15, distanceType = EV;
+string patternFile="", output="DMat";
+
+void printHelp(){
+ string help = "Alignment-free sequence comparison tool using spaced word frequencies"
+ "\nUsage: ./spaced [options] <sequence file>"
+ "\nOptions:"
+ "\n\t -h: print this help and exit"
+ "\n\t -o <file>: output file name (default: DMat)"
+ "\n\t -k <integer>: pattern weight (default 14)"
+ "\n\t -l <integer>: pattern don't care positions (default 15)"
+ "\n\t -n <integer>: number of patterns (default: 5)"
+ "\n\t -f <file>: use patterns in <file> instead of random patterns"
+ "\n\t -t <integer>: numer of threads (default: 25 threads)"
+ "\n\t -r: don't consider the reverse complement"
+ "\n\t -d EU | JS | EV: change distance type to Euclidean, Jensen-Shannon, evolutionary distance (default: EV)"
+ "\n";
+ cout << help << endl;
+}
+
+void parsePattern(string filePath, vector<string>& patternSet){
+ ifstream file(filePath.c_str());
+ string str;
+ while (getline(file, str)){
+ patternSet.push_back(str);
+ }
+ file.close();
+}
+
+void parseParameters(int argc, char** argv){
+ // too few?
+ if (argc < 2) {
+ printHelp(); // print help
+ exit(0); // and exit
+ }
+ // only parameter -h?
+ if (argc == 2)
+ if (strcmp("-h", argv[1]) == 0) {
+ printHelp();
+ exit(0);
+ }
+ for (int i = 1; i < argc - 1; i++) { // first argument is run dir, last should be sequence file
+ // help & exit
+ if (strcmp("-h", argv[i]) == 0) {
+ printHelp();
+ exit(0);
+ }
+ else if (strcmp("-r", argv[i]) == 0) {
+ revComp=false;
+ }
+ // patterns that should be used
+ else if (strcmp("-f", argv[i]) == 0) {
+ patternFile = argv[++i];
+ if (i == argc - 1) {
+ cerr << "Flag -f needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // output matrix file
+ else if (strcmp("-o", argv[i]) == 0) {
+ output = argv[++i];
+ if (i == argc - 1) {
+ cerr << "Flag -o needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // auto-generate pattern weight
+ else if (strcmp("-k", argv[i]) == 0) {
+ weight = atoi(argv[++i]);
+ if (i == argc - 1) {
+ cerr << "Flag -k needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // auto-generate pattern length
+ else if (strcmp("-l", argv[i]) == 0) {
+ dontcare = atoi(argv[++i]);
+ if (i == argc - 1) {
+ cerr << "Flag -l needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // number of combined matrices
+ else if (strcmp("-n", argv[i]) == 0) {
+ number = atoi(argv[++i]);
+ if (i == argc - 1) {
+ cerr << "Flag -n needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // number of threads
+ else if (strcmp("-t", argv[i]) == 0) {
+ threads = atoi(argv[++i]);
+ if (i == argc - 1) {
+ cerr << "Flag -t needed an argument!" << endl;
+ exit(-1);
+ }
+ }
+ // distance type
+ else if (strcmp("-d", argv[i]) == 0) {
+ i++;
+ if (i == argc - 1) {
+ cerr << "Flag -d needed an argument!" << endl;
+ exit(-1);
+ }
+ if (strcmp("eu", argv[i]) == 0 || strcmp("EU", argv[i]) == 0) {
+ distanceType = EU;
+ } else if (strcmp("js", argv[i]) == 0 || strcmp("JS", argv[i]) == 0) {
+ distanceType = JS;
+ } else if (strcmp("ev", argv[i]) == 0 || strcmp("EV", argv[i]) == 0) {
+ distanceType = EV;
+ }
+ }
+ // unknown flag
+ else {
+ printHelp();
+ cerr << "Unknown flag: " << argv[i] << endl;
+ exit(-1);
+ }
+ }
+}
+
+unsigned char* readData(const char * const filename, unsigned long long& n) {
+
+ FILE *file;
+ if (!(file = fopen(filename, "r"))) {
+ printf("Unable to open file %s \n", filename);
+ return NULL;
+ }
+ fseek (file , 0 , SEEK_END);
+ n = ftell (file);
+ unsigned char *result = new unsigned char[n*2];
+ rewind(file);
+ if (n > fread(result, sizeof(unsigned char), n, file)) {
+ printf("Error reading file %s \n", filename);
+ fclose(file);
+ delete[] result;
+ return NULL;
+ }
+ fclose(file);
+ return result;
+}
+
+int main(int argc, char *argv[]){
+ parseParameters(argc, argv);
+ char *inputFile = argv[argc - 1];
+ unsigned long long n = 0;
+ unsigned char *str = readData(inputFile, n);
+ if (str == NULL)
+ return 0;
+ vector<string> patternSet;
+ if(patternFile.length()==0){
+ variance* variance_obj;
+ int* laenge = new int[2];
+ laenge[0]=laenge[1]=dontcare+weight;
+ variance_obj = new variance(number,laenge, weight);
+ variance_obj->Silent();
+ variance_obj->Improve(500); //5000 mal versuchen ein Pattern zu verbessern
+ patternSet = variance_obj->GetPattern();
+ }
+ else
+ parsePattern(patternFile, patternSet);
+ double seq=0;
+ int length=0;
+ double dna=0;
+ char c;
+ for(unsigned long long i=0; i<n;i++){
+ c=toupper(str[i]);
+ if(c=='>'){
+ while(str[++i]!='\n');
+ seq++;
+ }
+ if(isalpha(c)){
+ switch(c){
+ case 'A':
+ dna++;
+ break;
+ case 'C':
+ dna++;
+ break;
+ case 'G':
+ dna++;
+ break;
+ case 'T':
+ dna++;
+ break;
+ case 'N':
+ dna++;
+ break;
+ }
+ length++;
+ }
+ }
+ if(dna/length>0.9){
+ if(seq>0x7FFF){
+ cout << "Too many sequences" << endl;
+ exit(0);
+ }
+ if(weight>32){
+ cout << "Maximum weight for DNA sequences is 32" << endl;
+ exit(0);
+ }
+ cout << seq << " DNA sequences read of total length "<< length << endl;
+ if(n>=0xffffffff)
+ spacedDNA<unsigned long long>(patternSet, str,n, distanceType, threads, weight, dontcare, revComp, output);
+ else
+ spacedDNA<unsigned int>(patternSet, str,n, distanceType, threads,weight, dontcare, revComp, output);
+ }
+ else{
+ if(seq>0xFFFF){
+ cout << "Too many sequences" << endl;
+ exit(0);
+ }
+ if(weight>12){
+ cout << "Maximum weight for protein sequences is 12" << endl;
+ exit(0);
+ }
+ cout << seq << " protein sequences read of total length "<< length << endl;
+ if(distanceType==EV){
+ cout << "Evolutionary distance not available for protein sequences"<< endl;
+ exit(0);
+ }
+ if(n>=0xffffffff)
+ spacedProt<unsigned long long>(patternSet, str,n, distanceType, threads, weight, dontcare, output);
+ else
+ spacedProt<unsigned int>(patternSet, str,n, distanceType, threads, weight, dontcare, output);
+ }
+}
+
+
diff --git a/variance.cpp b/variance.cpp
new file mode 100644
index 0000000..9a3d887
--- /dev/null
+++ b/variance.cpp
@@ -0,0 +1,788 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * variance object file
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#include "variance.h"
+
+
+/*---Variables---------------------------------------------------------------*/
+
+
+
+/*===Main-Part===============================================================*/
+/*---Constructor-------------------------------------------------------------*/
+/**
+ * Default constructor, sets the default vaulues, pattern will be generated automatically.
+ */
+variance::variance() {
+ int *length = new int[2];
+ length[0] = 14;
+ length[1] = 14;
+ Ctor(NULL, 10, length, 8, 10000, 10000, 10000, 0.75, 0.25, 64);
+}
+
+
+/**
+ * File constructor, sets values only from files; resets automatically if there are problems.
+ */
+variance::variance(char* pattern_file) {
+ int *length = new int[2];
+ length[0] = 14;
+ length[1] = 14;
+ Ctor(pattern_file, 10, length, 8, 10000, 10000, 10000, 0.75, 0.25, 64);
+}
+
+
+/**
+ * Short constructor, sets some default vaulues, just pattern dimension is set.
+ *
+ * @param size The amount of patterns; pattern number.
+ *
+ * @param length The pattern length for each pattern of the pattern set.
+ *
+ * @param weigth The weight (match positions; '1') for each pattern of the pattern set.
+ */
+variance::variance(int size, int *length, int weight) {
+ Ctor(NULL, size, length, weight, 10000, 10000, 10000, 0.75, 0.25, 64);
+}
+
+
+/**
+ * Long constructor, sets the values; resets automatically, if there are problems.
+ *
+ * @param pattern_file File, that may contains submitted pattern
+ *
+ * @param align_file File, that may contains an alignment file to estimate p, q, l_hom, li and lj
+ *
+ * @param size The amount of patterns; pattern number.
+ *
+ * @param length The pattern length for each pattern of the pattern set.
+ *
+ * @param weigth The weight (match positions; '1') for each pattern of the pattern set.
+ *
+ * @param l_hom In theory, the amount of homologous positions of two sequences in an multiple alignment.
+ *
+ * @param l1 In theory, the first(represents each sequence i) sequence length of two observed sequences.
+ *
+ * @param l2 In theory, the second(represents each sequence j) sequence length of two observed sequences.
+ *
+ * @param p The match probability ( = #matches / #l_hom)
+ *
+ * @param q The background probability for each nucleotide A,C,G,T
+ */
+variance::variance(char* pattern_file, int size, int *length, int weight, int l_hom, int l1, int l2, double p, double q, int H) {
+ Ctor(pattern_file, size, length, weight, l_hom, l1, l2, p, q, H);
+}
+
+
+/**
+ * Default destructor, deletes all vectors and matrices in the object.
+ */
+variance::~variance() {
+ Clear();
+}
+
+void variance::Ctor(char* pattern_file, int size, int *length, int weight, int l_hom, int l1, int l2, double p, double q, int H){
+ this->l_hom = l_hom;
+ this->l1 = l1;
+ this->l2 = l2;
+ this->p = p;
+ this->q = q;
+ this->H = H;
+ this->length = length;
+ variance_val = 0;
+ loop_opt = 1;
+ quiet = false;
+ silent = false;
+ oc = false;
+ update = false;
+ outvar = "variance: ";
+ InitVariance(pattern_file, size, length, weight);
+}
+
+/*---Init--------------------------------------------------------------------*/
+/**
+ * Complete Initializing of the Variance, which also creates the patternsets
+ * calculates the first variance/oc and checks if each parameter is in its
+ * own correct size and dimension.
+ */
+void variance::InitVariance(char* pattern_file, int size, int *length, int weight){
+ Clear();
+ pattern_set = new patternset(pattern_file, size, length, weight);
+ this->size = pattern_set->GetSize();
+ this->weight = pattern_set->GetWeight();
+ this->length = pattern_set->GetLength();
+ length_mean = pattern_set->GetLengthMean();
+ improve = pattern_set->GetImprove();
+ update = pattern_set->GetUpdate();
+ if (q > p || q >= 1 || p > 1 || q < 0 || p < 0) { /*Incorrect values for variance have to be corrected*/
+ SecureMessage("pq");
+ p = 0.9;
+ q = 0.25;
+ update = true;
+ }
+
+ if (l_hom <= 0 || l1 <= 0 || l2 <= 0) {
+ SecureMessage("length");
+ l_hom = 10000;
+ l1 = 10000;
+ l2 = 10000;
+ update = true;
+ }
+ InitVarMatrix();
+ InitVar();
+ if(update){
+ Quiet();
+ }
+}
+
+/**
+ * Part Initialzinig of the Variance, it will not reset complete.
+ * Only for ReInitializin a new pattternset with most values in
+ * common.
+ * Public Access, so other programs can 'reload' it
+ */
+void variance::ReInitVariance(){
+ for(int i = 0; i < size; i++){
+ var_sum[i].clear();
+ }
+ var_sum.clear();
+ pattern_order.clear();
+ pattern_set->ReInitPattern();
+ length_mean = pattern_set->GetLengthMean();
+ InitVarMatrix();
+ InitVar();
+}
+
+
+/**
+ * Part of the destructor; could be used in other cases for
+ * reinitializing the variance/oc new.
+ */
+void variance::Clear(){
+ for(int i = 0; i < size; i++){
+ var_sum[i].clear();
+ }
+ var_sum.clear();
+ pattern_order.clear();
+ delete pattern_set;
+}
+
+/**
+ * Up to different patternsizes a specific matrix will be created.
+ * Necessary for reinitliazing.
+ */
+void variance::InitVarMatrix() {
+ std::vector<double> tmp;
+ for (int i = 0; i < size; i++) {
+ tmp.push_back(0.0);
+ }
+ for (int i = 0; i < size; i++) {
+ var_sum.push_back(tmp);
+ }
+}
+
+/**
+ * First Initializing of the variance/oc where each value of the matrix has
+ * to be created, instead of updated variance, where it is not necessary to
+ * recalculate each value/entry. Only of those, which alter due to pattern
+ * modification.
+ */
+void variance::InitVar() {
+ double var_hom, var_bac, hom, back, tmp;// ,hom_val, back_val;
+ int shift, lower_length, upper_length;
+ uint64_t oc_hom;
+
+ oc_hom = 1;
+ //hom_val = 0;
+ //back_val = 0;
+
+ for (int i = 0; i < size; i++) { /*i and j represents Pi and Pj of the set of pattern*/
+ lower_length = (int) log2(pattern_set->GetPattern(i));
+ for (int j = i; j < size; j++) {
+ upper_length = (int) log2(pattern_set->GetPattern(j)) + 1;
+
+ if(i == j && !oc){
+ upper_length = 1;
+ }
+
+ var_hom = 0.0;
+ var_bac = 0.0;
+ for (int s = -1 * lower_length; s < upper_length; s++) { /*As in the formula, the shift goes from max shift left to max shift right*/
+ shift = ShiftPos(i, j, s); /*At least one position has to overlap*/
+ if (!oc) {
+ var_hom += (pow(p, shift) - pow(p, 2 * weight)); /*summation of the homologue first part*/
+ var_bac += (pow(q, shift) - pow(q, 2 * weight));
+ }
+ else{
+ var_hom += (double) (oc_hom << (shift));
+ }
+ }
+ if(!oc){
+ hom = (l_hom - length_mean + 1);
+ hom *= var_hom;
+
+ back = (l1 - length_mean + 1);
+ back *= (l2 - length_mean);
+ back *= var_bac;
+ }
+ else{
+ hom = var_hom;
+ back = 0.0;
+ }
+ // hom_val += hom;
+ // back_val += back;
+ var_sum[i][j] = hom + back; /*For each pair Pi and Pj this is the direct contribute to the complete variance*/
+ var_sum[j][i] = var_sum[i][j];
+ }
+ }
+ variance_val = Variance();
+
+ pattern_order.clear();
+
+ for(int i = 0; i < size; i++){
+ tmp = 0.0;
+ for(int j = 0; j < size; j++){
+ tmp += var_sum[i][j]; /*Sets the contribute of each pattern in the beginning*/
+ }
+ pattern_order.push_back(tmp);
+ }
+ ResetPatternOrder();
+
+ //std::cout << "homologue contribute: " << hom_val/variance_val << std::endl;
+ //std::cout << "background contribute: " << back_val/variance_val << std::endl;
+}
+
+
+/*---Variance-----------------------------------------------------------------*/
+/**
+ * IF neccessary, it calculates the current variance/oc from the matrix.
+ *
+ * @return returns current variance/oc
+ */
+double variance::Variance() {
+ double var;
+
+ var = 0.0;
+
+ for (int i = 0; i < size; i++) {
+ for(int j = i; j < size; j++){
+ var += var_sum[i][j];
+ }
+ }
+ variance_val = var;
+ return var;
+}
+
+/**
+ * Instead of recalculating each value, only the matrix entries which
+ * belong to the modified pattern have to be recalculated.
+ *
+ * @param pat Patternindex of the modified pattern
+ */
+void variance::UpdateVariance(int pat) {
+ double var_hom, var_bac, hom, back;
+ int shift, lower_length, upper_length;
+ uint64_t oc_hom;
+
+ oc_hom = 1;
+
+ for (int j = 0; j < size; j++) {
+ var_hom = 0.0;
+ var_bac = 0.0;
+
+ lower_length = (int) log2(pattern_set->GetPattern(pat));
+ upper_length = (int) log2(pattern_set->GetPattern(j)) + 1;
+
+ if(pat == j && !oc){
+ upper_length = 1;
+ }
+
+ for (int s = -1 * lower_length; s < upper_length; s++) { /*As in the formula, the shift goes from max shift left to max shift right*/
+ shift = ShiftPos(pat, j, s); /*At least one position has to overlap*/
+ if (!oc) {
+ var_hom += (pow(p, shift) - pow(p, 2 * weight)); /*summation of the homologue first part*/
+ var_bac += (pow(q, shift) - pow(q, 2 * weight));
+ }
+ else{
+ var_hom += (double) (oc_hom << (shift));
+ }
+ }
+
+ variance_val -= var_sum[pat][j];
+ pattern_order[pat] -= var_sum[pat][j];
+ pattern_order[j] -= var_sum[pat][j];
+
+ if(!oc){
+ hom = (l_hom - length_mean + 1);
+ hom *= var_hom;
+ back = (l1 - length_mean + 1);
+ back *= (l2 - length_mean);
+ back *= var_bac;
+ }
+ else{
+ hom = var_hom;
+ back = 0.0;
+ }
+
+ var_sum[pat][j] = hom + back; /*For each pair Pi and Pj this is the direct share of the complete variance...*/
+ var_sum[j][pat] = var_sum[pat][j];
+
+ variance_val += var_sum[pat][j];
+ pattern_order[pat] += var_sum[pat][j];
+ pattern_order[j] += var_sum[pat][j];
+ }
+}
+
+/**
+ * Shifts pattern and counts the number of common match positions
+ *
+ * @param p1 Position of the first used pattern of the pattern set
+ *
+ * @param p2 Position of the second used pattern of the pattern set
+ * NOTE: possible is p1 = p2
+ *
+ * @param s The shift of the second pattern, s < 0 := shift left pattern 2, s > := shift right pattern 2
+ *
+ * @return Calculates and returns current variance
+ */
+int variance::ShiftPos(int number, int number2, int s) {
+ int counter, maxi, maxa, maxb;
+ uint64_t pata, patb, c;
+
+ pata = pattern_set->GetPattern(number);
+ patb = pattern_set->GetPattern(number2);
+
+ counter = 0;
+ c = 1; /*LookUp Counter*/
+
+ if (s < 0) { /*Negative shift means, the upper pattern is shifted right*/
+ s = 0 - s;
+
+ pata = pata >> s;
+ }
+ else{ /*otherweise the lower pattern is shifted right*/
+ patb = patb >> s;
+ }
+
+ maxa = (int) log2(pata)+1;
+ maxb = (int) log2(patb)+1;
+ maxi = std::min(maxa,maxb);
+ for (int i = 0; i < maxi; i++) {
+ if(((pata & patb) & c << i) != 0){
+ counter++;
+ }
+ }
+
+ if(!oc){
+ counter = 2*weight - counter; /*Allows us to calculate th n(p,p',s) value in a way...*/
+ } /*...with not so much space. We do not need to look for...*/
+ return counter; /*...each position containing at least one match*/
+}
+
+/**
+ * Returns the position of the worst pattern, estimated by the maximum variancepart
+ * for each pattern pair
+ *
+ * @param n The pattern index of the chosen pattern.
+ *
+ * @return returns position worst matrix by max_value
+ */
+int variance::WorstPattern(int n){
+ return pattern_order_sort[n]->GetPos();
+}
+
+void variance::ResetPatternOrder(){
+ extkey* tmp;
+ int end;
+
+ if(pattern_order_sort.size() != 0){
+ pattern_order_sort.clear();
+ }
+
+ if(length[0]==weight){
+ end = size-1;
+ }
+ else{
+ end = size;
+ }
+
+ for(int i = 0; i < end; i++){
+ tmp = new extkey(i,pattern_order[i]); /*Using an extended key to store each contribute and the ...*/
+ pattern_order_sort.push_back(tmp); /*pattern number.*/
+ }
+
+ std::sort(pattern_order_sort.rbegin(), pattern_order_sort.rend());
+}
+
+/**
+ * The improvement method
+ *
+ * @param limit The number of patterns which have to be selected and mayby modified.
+ */
+void variance::Improve(int limit) {
+ std::vector<uint64_t> vec;
+ uint64_t pat_save;
+ double var_save, tmp_var, best_by;
+ int worst_pat, better_pattern, counter;
+ //std::string str_out = "Output_variance";
+ //std::ofstream fout;
+ //fout.open(str_out);
+
+ best_by = variance_val;
+
+ if (!silent) {
+ std::cout << "First " << outvar << "\t" << GetVariance() << std::endl;
+ std::cout << "First norm_" << outvar << "\t" << GetNormVariance() << std::endl << std::endl;
+ }
+ //fout << 0 << "\t" << GetNormVariance() << std::endl;
+
+ if(loop_opt > 1){
+ pattern_set->Silent();
+ for(int i = 0; i < size; i++){
+ vec.push_back(pattern_set->GetPattern(i));
+ }
+ }
+
+
+ for(int d = 0; d < loop_opt; d++){
+
+ counter = 0;
+ better_pattern = 0;
+ worst_pat = 0;
+ if(d > 0){
+ if(!silent && !quiet){
+ std::cout << "\n\n=================\nAnd do it again!" << std::endl;
+ }
+ ReInitVariance();
+ }
+ if (improve) {
+ for (int i = 1; i <= limit; i++) {
+ var_save = variance_val;
+ worst_pat = GetWorstPat(counter%((int)pattern_order_sort.size()));
+ //std::cout << counter << "\t" << worst_pat << std::endl;
+ pat_save = pattern_set->GetPattern(worst_pat);
+ pattern_set->ChangeBits(worst_pat);
+
+ UpdateVariance(worst_pat);
+ tmp_var = variance_val;
+
+ if (tmp_var < var_save && pattern_set->UniqPattern(worst_pat)) {
+ better_pattern++;
+ var_save = variance_val;
+ if (quiet && !silent) {
+ std::cout << "\r*** BETTER PATTERN " << better_pattern << " ***";
+ std::cout.flush();
+ }
+ else if (!quiet && !silent) {
+ std::cout << "*** BETTER PATTERN " << better_pattern << " ***" << std::endl;
+ std::cout << "Step " << i << " / " << limit << std::endl << "Patternset: \n";
+ pattern_set->Print();
+ std::cout << outvar << GetVariance() << std::endl;
+ std::cout << "norm_" << outvar << GetNormVariance() << std::endl << std::endl;
+ }
+ else{
+ //Do-Nothing
+ }
+ counter = 0;
+ ResetPatternOrder();
+ //fout << i << "\t" << GetNormVariance() << std::endl;
+ }
+ else{
+ pattern_set->SetPattern(worst_pat, pat_save);
+ UpdateVariance(worst_pat);
+ counter++;
+ }
+ }
+ if(best_by > variance_val && loop_opt > 2){
+ best_by = variance_val;
+ for(int k = 0; k < size; k++){
+ vec[k] = pattern_set->GetPattern(k);
+ }
+ }
+ if (!silent && loop_opt < 2) {
+ std::cout << "\n\nBest pattern:\n\n";
+ pattern_set->Print();
+ std::cout << "\nBest " << outvar << GetVariance() << std::endl;
+ std::cout << "Best norm_"<< outvar << GetNormVariance() << std::endl;
+ }
+ //fout << limit << "\t" << GetNormVariance() << std::endl;
+ }
+ }
+ //if(update){
+ //SecureMessage("update");
+ //}
+ if(loop_opt > 1){
+ for(int i = 0; i < size; i++){
+ pattern_set->SetPattern(i, vec[i]);
+ }
+ variance_val = best_by;
+ if(!silent){
+ std::cout << "\n\nBest pattern:\n\n";
+ pattern_set->Print();
+ std::cout << "\nBest " << outvar << GetVariance() << std::endl;
+ std::cout << "Best norm_"<< outvar << GetNormVariance() << std::endl;
+ }
+ }
+ vec.clear();
+ //fout.close();
+}
+
+
+/**
+ * Turns the current calculation from variance to overlap complexity.
+ * The matrix has to be recalculated complete.
+ */
+void variance::ImproveOC() {
+ this->oc = true;
+ outvar = "oc: ";
+ InitVar();
+}
+
+/*---stuff-------------------------------------------------------------------*/
+/**
+ * Method calculate the number of all pattern combinations
+ * Actually the gauss summation (n*(n+1)/2)
+ *
+ * @return the number of maximum pattern combinations
+ */
+double variance::Gauss() {
+ return 0.5*size*size + 0.5*size;
+}
+
+/**
+ * Changes to quiet output, only the number of better patterns will be printed
+ * automatically, and the best patternset by the variance/oc minimization.
+ * Errors arr going to be printed!
+ */
+void variance::Quiet() {
+ this->quiet = true;
+}
+
+/**
+ * Changes to silent output, nothing will be printed automatically!
+ * Errors are going to be printed!
+ */
+void variance::Silent() {
+ this->quiet = true;
+ this->silent = true;
+ pattern_set->Silent();
+}
+
+/**
+ * Resets the pattern lengths from evenly distributed lengths between maximum
+ * and minimum to random created pattern lengths.
+ * One quarter of the set will have maximum length, the rest will have the
+ * length between maximum and minimum randomly distributed.
+ */
+void variance::RandPatLength(){
+ pattern_set->RandPatLength();
+ InitVar();
+}
+
+/**
+ * Sets the number of times the improvement method should be executed again.
+ * It will save the best variance/oc and its patternset.
+ *
+ * @param n The number of iterations for the improvement method
+ */
+void variance::LoopOpt(int n){
+ if(n == 0){
+ n = 1;
+ }
+ loop_opt = n;
+
+}
+
+/**
+ * Prints complete the current patternset
+ */
+void variance::PrintPattern(){
+ pattern_set->Print();
+}
+
+/**
+ * A Method to collect all errormessages. Just easier for programmer to change
+ * the text or extend.
+ *
+ * @param errmsg
+ * Due to a few possible errormessages, this is the option, which has to be printed.
+ *
+ * @param pos
+ * The position of the incorrect patterns.
+ *
+ */
+void variance::SecureMessage(std::string errmsg) {
+ /*if (errmsg == "noimprove") {
+ std::cerr << "Using your pattern conditions it is not sensible to improve your pattern, sorry!" << std::endl;
+ std::cerr << "Deactivating improve mode\n" << std::endl;
+ return;
+ }
+ if (errmsg == "update"){
+ std::cerr << "\n\n--> IMPORTANT <--\nDue to some configuration errors, your submitted parameters have been updated!\n" << std::endl;
+ return;
+ }
+ if (errmsg == "pq") {
+ std::cerr << "Error while parsing your p and/or q value: \t0 < q <= p <= 1!" << std::endl;
+ std::cerr << "Return to default values:\tp = 0.9 \tq=0.25\n" << std::endl;
+ return;
+ }
+ if (errmsg == "length") {
+ std::cerr << "Error while parsing your sequence length S: \t 0 < S!" << std::endl;
+ std::cerr << "Return to default value: \tS = 10000\n" << std::endl;
+ return;
+ }
+ if (errmsg == "wrongindex") {
+ std::cerr << "ERROR! Pattern does not exist... do nothing\n" << std::endl;
+ return;
+ }
+ if (errmsg == "speed"){
+ std::cerr << "\n\n\nUPDATE! SpEED lenghts optimization is used, creating a new patternset!!!\n\n\n\n" << std::endl;
+ return;
+ }*/
+}
+
+/*---Set-&-GetFunc-----------------------------------------------------------*/
+
+std::vector<std::string> variance::GetPattern(){
+ return pattern_set->GetStringPattern();
+}
+
+std::string variance::GetPattern(int number){
+ return pattern_set->GetString(number);
+}
+
+/**
+ * Returns the current variance
+ *
+ * @return returns variance
+ */
+double variance::GetVariance() {
+ return variance_val;
+}
+
+/**
+ * Returns the current variance, normalized
+ *
+ * @return returns norm_variance
+ */
+double variance::GetNormVariance() {
+ return variance_val / Gauss();
+}
+
+
+/**
+ * Returns the position of the worst pattern, estimated by the maximum variancepart
+ * for each pattern pair
+ *
+ * @return returns position worst matrix by max_value
+ */
+int variance::GetWorstPat(int number) {
+ if(number >= size){
+ return 0;
+ }
+ return WorstPattern(number);
+}
+
+/**
+ * Returns the weight of each pattern, the match positions
+ *
+ * @return returns weight
+ */
+int variance::GetWeight() {
+ return pattern_set->GetWeight();
+}
+
+/**
+ * Returns the amount of patters; number of patterns
+ *
+ * @return returns size
+ */
+int variance::GetSize() {
+ return size;
+}
+
+/**
+ * Returns the length of each pattern
+ *
+ * @return returns length
+ */
+int* variance::GetLength() {
+ length = pattern_set->GetLength();
+ return length;
+}
+
+/**
+ * Returns the match probability
+ *
+ * @return returns p value
+ */
+double variance::GetP() {
+ return p;
+}
+
+/**
+ * Returns the background probabillity, summation over all nucleotids
+ *
+ * @return returns q value
+ */
+double variance::GetQ() {
+ return q;
+}
+
+/**
+ * Returns the length of the homologous sequence pair
+ *
+ * @return returns homologous positions
+ */
+int variance::GetLHom() {
+ return l_hom;
+}
+
+/**
+ * Returns the length of the first observed sequence
+ *
+ * @return returns length sequence 1
+ */
+int variance::GetL1() {
+ return l1;
+}
+
+/**
+ * Returns the length of the second observed sequence
+ *
+ * @return returns length sequence 2
+ */
+int variance::GetL2() {
+ return l2;
+}
+
+/**
+ * Returns the current format, variance or overlap complexity
+ *
+ * @return the format string, either 'variance' or 'oc'
+ */
+std::string variance::GetFormat() {
+ return outvar;
+}
+
+/**
+ * Returns the length of a random homolgue region on a dataset
+ *
+ * @return The length of random homolgue region
+ */
+int variance::GetH(){
+ return H;
+}
diff --git a/variance.h b/variance.h
new file mode 100644
index 0000000..458b100
--- /dev/null
+++ b/variance.h
@@ -0,0 +1,112 @@
+/**
+ * This programm calculates the variance/OC and/or the sensitivity of a set of pattern with the same weight.
+ * It is possible to improve your patternset and read patterns from a file.
+ *
+ * variance object header
+ *
+ * For theory please have a look at:
+ *
+ * B. Morgenstern, B. Zhu, S. Horwege, C.-A Leimeister (2015)
+ * Estimating evolutionary distances between genomic sequences from spaced-word matches
+ * Algorithms for Molecular Biology 10, 5. (http://www.almob.org/content/10/1/5/abstract)
+ *
+ *
+ * @author: Lars Hahn - 26.10.2015, Georg-August-Universitaet Goettingen
+ * @version: 1.0.2 11/2015
+ */
+#ifndef VARIANCE_H_
+#define VARIANCE_H_
+
+
+#include <iostream>
+#include <fstream>
+#include <random>
+#include <vector>
+#include <string>
+#include <algorithm>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <inttypes.h>
+#include "patternset.h"
+#include "extkey.h"
+
+class variance {
+public:
+ variance();
+ variance(char* pattern_file);
+ variance(int size, int *length, int weight);
+ variance(char* pattern_file, int size, int *length, int weight, int l_hom, int l1, int l2, double p, double q, int H);
+ ~variance();
+
+ void ReInitVariance();
+ void RandPatLength();
+
+ void Improve(int limit);
+ void ImproveOC();
+
+ std::vector<std::string> GetPattern();
+ std::string GetPattern(int number);
+ int GetWorstPat(int number);
+ int GetWeight();
+ int GetSize();
+ int* GetLength();
+ void PrintPattern();
+
+ std::string GetFormat();
+ double GetVariance();
+ double GetNormVariance();
+ double GetP();
+ double GetQ();
+ int GetLHom();
+ int GetL1();
+ int GetL2();
+ int GetH();
+
+ void LoopOpt(int n);
+ void Quiet();
+ void Silent();
+
+protected:
+ void Ctor(char* pattern_file, int size, int *length, int weight, int l_hom, int l1, int l2, double p, double q, int H);
+ void InitVariance(char* pattern_file, int size, int *length, int weight);
+ void InitVarMatrix();
+ void InitVar();
+ void SetPatOrder();
+ void Clear();
+
+ double Variance();
+ void UpdateVariance(int pat);
+
+ int ShiftPos(int number, int number2, int s);
+ int WorstPattern(int number);
+ void ResetPatternOrder();
+
+ double Gauss();
+ void SecureMessage(std::string errmsg);
+
+private:
+ std::vector<std::vector<double> > var_sum;
+ patternset* pattern_set;
+ std::vector<double> pattern_order;
+ std::vector<extkey*> pattern_order_sort;
+ std::string outvar;
+ double variance_val;
+ double p;
+ double q;
+ int size;
+ int weight;
+ int length_mean;
+ int *length;
+ int l_hom;
+ int l1;
+ int l2;
+ int H;
+ int loop_opt;
+ bool oc;
+ bool improve;
+ bool quiet;
+ bool silent;
+ bool update;
+};
+#endif
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/spaced.git
More information about the debian-med-commit
mailing list