[med-svn] [bagpipe] 04/06: New upstream version 2012.02.15
Andreas Tille
tille at debian.org
Sat Dec 2 08:23:42 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository bagpipe.
commit d0e28521417525a1e7b28df096378f814675f6f9
Author: Andreas Tille <tille at debian.org>
Date: Sat Dec 2 09:22:25 2017 +0100
New upstream version 2012.02.15
---
GNUcopyleft.txt | 675 ++++++++++++++++++++++
bagpipe/bagpipe.pl | 657 ++++++++++++++++++++++
bagpipe/lib/Algorithm.pm | 467 ++++++++++++++++
bagpipe/lib/Array.pm | 1111 +++++++++++++++++++++++++++++++++++++
bagpipe/lib/Array2d.pm | 214 +++++++
bagpipe/lib/Assert.pm | 139 +++++
bagpipe/lib/CCloneable.pm | 28 +
bagpipe/lib/CConfigFile.pm | 172 ++++++
bagpipe/lib/CDataFrame.pm | 920 ++++++++++++++++++++++++++++++
bagpipe/lib/CRing.pm | 130 +++++
bagpipe/lib/Debug.pm | 68 +++
bagpipe/lib/FileIoHelper.pm | 461 +++++++++++++++
bagpipe/lib/Hash.pm | 168 ++++++
bagpipe/lib/NamedArgs.pm | 113 ++++
bagpipe/lib/Number.pm | 34 ++
bagpipe/lib/StdDefs.pm | 28 +
bagpipe/lib/String.pm | 90 +++
bagpipe/lib/VarType.pm | 67 +++
bagpipe/lib/fitgev.R | 63 +++
bagpipe/lib/nullscangenome.R | 55 ++
bagpipe/lib/plotscans.R | 119 ++++
bagpipe/lib/posboot.R | 87 +++
bagpipe/lib/scangenome.R | 50 ++
bagpipe/lib/scanloci.R | 57 ++
bagpipe/lib/write_nomenu_config.R | 25 +
debian/bagpipe.1 | 150 -----
debian/bagpipe.install | 3 -
debian/bagpipe.manpages | 1 -
debian/changelog | 6 -
debian/compat | 1 -
debian/control | 39 --
debian/copyright | 27 -
debian/rules | 16 -
debian/script/bagpipe | 4 -
debian/source/format | 1 -
debian/upstream/metadata | 22 -
debian/watch | 4 -
37 files changed, 5998 insertions(+), 274 deletions(-)
diff --git a/GNUcopyleft.txt b/GNUcopyleft.txt
new file mode 100644
index 0000000..1b6fc8f
--- /dev/null
+++ b/GNUcopyleft.txt
@@ -0,0 +1,675 @@
+
+ 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/bagpipe/bagpipe.pl b/bagpipe/bagpipe.pl
new file mode 100755
index 0000000..81d2520
--- /dev/null
+++ b/bagpipe/bagpipe.pl
@@ -0,0 +1,657 @@
+#!/usr/bin/perl
+
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+# TODO: replace window.radius with peak.separation everywhere!!!
+BEGIN
+{
+ unless ($ENV{'BAGPIPE_LIBS'})
+ {
+ die "Error: shell variable \$BAGPIPE_LIBS must be defined.\n";
+ }
+}
+# General packages
+use strict;
+use English;
+use Data::Dumper;
+use Getopt::Long;
+use warnings;
+# My packages
+use lib $ENV{'BAGPIPE_LIBS'};
+use Assert;
+use StdDefs;
+use CConfigFile;
+use CDataFrame;
+use File::Copy;
+use File::Path;
+use FileIoHelper qw(OpenFilesForReading OpenFilesForWriting
+ FileNotEmpty FileEmptyOrAbsent FileExists
+ GetFileBasename GetFileSuffix
+ GetDirContents );
+
+# Global variables accessed by functions
+my $DRY_RUN = false;
+my $CLEAN_UP = true;
+my $DEVEL = false;
+my $SHOW_WARNINGS = true;
+my $CAUGHT_R_INPUT_ERROR_PATTERN = "Bagphenotype .* Error";
+
+# Documentation ----------------------------------------
+
+my $VERSION_STRING = <<EOV;
+--------------------------------------------------------------------
+Bagpipe - a pipeline for genetic association using HAPPY genome data
+Version: 2012_02_14
+Author: William Valdar <william.valdar\@unc.edu>.
+Cite: Valdar et al (2009) Genetics 182(4):1263-77.
+Please note: this is software is distributed under a GNU copyleft
+license, is not formally validated, and comes with no guarantees
+whatsoever.
+--------------------------------------------------------------------
+EOV
+
+my $HELP_USAGE_STRING =<<EOU;
+Type 'bagpipe -h' for summary usage
+ 'bagpipe --help' for detailed usage
+ 'bagpipe --examples' for examples
+EOU
+
+my $BASIC_USAGE_STRING = <<EOU;
+bagpipe configfile [-scan] [-perm] [-perm_collate]
+ [-nullsim] [-nullsim_collate] [-plot] [other options]
+EOU
+
+my $SUMMARY_USAGE_STRING = <<EOU;
+Usage:
+ bagpipe configfile [-action(s)] [-option(s)=argument]
+
+actions: scan, perm, perm_collate, nullsim, nullsim_collate,
+ plot
+
+options: locus_groups, outdir, scandir
+ plot_height, plot_width, plot_scale, plot_type,
+ plot_score, plot_tailprob
+ memory, overwrite
+EOU
+
+my $DETAILED_USAGE_STRING = << 'EOU';
+
+Usage:
+ bagpipe configfile [options ...]
+
+Requirements:
+Config file: a file of configuration parameters
+ HAPPY genome cache: a directory containing cached HAPPY matrices
+ Phenotype file: a tab delimited table with phenotype and covariate
+ values for each animal. The table must include
+ 1) a column SUBJECT.NAME giving the id of the animal in
+ the corresponding HAPPY data;
+ 2) a column with the name of the phenotype as listed in the
+ Phenotype column of the model menu file;
+ 3) columns corresponding to any covariates used in the
+ model menu formulae.
+Key:
+n - an integer
+d - a decimal
+s - a string
+b - 1 for true and 0 for false.
+xL - a comma separated list of type x
+locus.group - either the name of a chromosome (eg, 1 or X) or a file of marker names.
+
+Actions:
+ scan - Genome scans
+ perm - Permutations
+ perm_collate - Calculate permutation thresholds
+ nullsim - Null simulations (not yet implemented)
+ nullsim_collate - Calculate null simulation thresholds (not yet implemented)
+ plot - Plot QTL scans to pdf file.
+
+ scanloci=sL - [Experimental] Scan individual loci, saving all fit information
+ posboot=sL - [Experimental] Perform k positional bootstraps between marker
+ intervals m1 and m2, where s=m1,m2,k
+
+Settings:
+ outdir - Directory for output. Default is ./
+ scandir - Directory for output of scan data files within outdir (./)
+ locus_groups=sL - Run analyses on certain locus groups only [Default: do all in configfile]
+ save.at.loci=sL - [Experimental] Save model fit when scanning at specified loci.
+ tailprobs=dL - Calculate nullscan tail probabilities (ie, significance thresholds).
+ Eg, 0.05,0.01. [Default nL=0.001,0.01,0.05,0.1,0.2]
+
+Plotting:
+ plot_box - Draw box around each plot
+ plot_scale=s - Scale used for plot summary s={cM,bp,Mb}
+ plot_score=s - Specify locus score to plot s={LOD,modelcmp}
+ plot_type=s - s = {chromosomes,genome}
+ plot_tailprob=sL - Include threshold lines as calculated by the *_collate options.
+ Each string s should have format <nullscantype>_<d>, where
+ nullscantype={perm,nullsim} and <d> is the tailprob looked up in the
+ corresponding thresh file. Eg, s=perm_0.05 plots the 0.05 threshold
+ calculated by permutation.
+ plot_height=d - Notional height of plotsummary output in inches
+ plot_width=d - Notional width of plotsummary output in inches
+
+Performance:
+ cleanup - Delete Rout files that do not report errors or warnings
+ dryrun - Print commands but do not execute them
+ memory=d - Store marker matrices in memory upto n Mb [default n=0]
+ nomenuconfig - Write the complete config file for an analysis
+ ie, incorporating all the information from the menu file
+ but minus the menu file
+ overwrite=b - Overwrite existing analysis files [Default overwrite=0]
+
+Help:
+ examples - Give examples
+ example_config - Write an example config file to config.example
+
+EOU
+
+# Functions --------------------------------------------
+
+sub WriteBlankConfigFile($)
+{
+ my $ost = shift;
+
+ $ost->print(<<END_OF_CONFIG);
+genome.cache.dir ../SaveMatricesFromHAPPY
+phenotype.file ./PhenotypesPlusCovariates.txt
+analysis.id MyPhenotypeModelName
+scan.formula.null anxiety.score ~ sex
+scan.formula.test anxiety.score ~ sex + interval.additive(THE.LOCUS)
+scan.function linear
+num.nullscans 200
+nullscan.seed 1
+locus.groups 1,2,X,FileOfOtherMarkers.txt
+END_OF_CONFIG
+}
+
+sub ErrorUsage($)
+{
+ my $errorMessage=shift;
+ print STDERR "$errorMessage\n$HELP_USAGE_STRING\n";
+}
+
+sub BasicUsage()
+{
+ print STDERR "$VERSION_STRING\n$BASIC_USAGE_STRING\n$HELP_USAGE_STRING\n";
+}
+
+sub SummaryUsage()
+{
+ print STDERR "$VERSION_STRING\n$SUMMARY_USAGE_STRING\n$HELP_USAGE_STRING\n";
+}
+
+sub LongUsage()
+{
+ print STDERR "$VERSION_STRING\n$DETAILED_USAGE_STRING\n";
+}
+
+sub PrintExamples()
+{
+ print << 'EOE';
+Bagpipe Examples:
+
+bagpipe configfile --scan
+ Run genome scans on the locus.groups listed in configfile
+
+bagpipe config --scan --locus_groups=1,2
+ Scan only chr 1 and 2.
+
+bagpipe configfile --scan --dryrun
+ Write commands to screen but do not execute.
+
+bagpipe config --perm_collate --tailprobs=0.025
+ Perform permutations and calculate only the upper p=0.025 quantile from an existing set of permutations
+
+bagpipe configfile --scan --perm --perm_collate --plot --plot_score=modelcmp --plot_tailprob=perm_0.05,perm_0.01
+ Scan all locus groups, do permutation scans, calculate thresholds a
+
+EOE
+}
+
+sub SystemR
+{
+ my ($cmd, $rout, $errorMessage, $force) = @ARG;
+ $cmd = "Rscript --vanilla $cmd "
+ .($DEVEL ? "--devel " : "")
+ ."> $rout 2>&1";
+ if ($DRY_RUN and not $force)
+ {
+ print "$cmd\n";
+ return;
+ }
+ print "Executing: $cmd\n";
+ system($cmd);
+ if (my $msg = CheckRout($rout, $CLEAN_UP))
+ {
+ if ($msg =~ m/$CAUGHT_R_INPUT_ERROR_PATTERN/)
+ {
+ die "Bagphenotype Caught Error: $errorMessage. More info:\n$msg\n";
+ }
+ else
+ {
+ die "Runtime Error from R script: $errorMessage. See $rout for details.\nLast lines from $rout:\n$msg\n";
+ }
+ }
+}
+
+sub CheckRout($$)
+{
+ my($file, $cleanUp) = @ARG;
+ my $ist = OpenFilesForReading($file);
+ my $hasWarnings = false;
+ while (my $line = $ist->getline)
+ {
+ if ($line =~ s/.*($CAUGHT_R_INPUT_ERROR_PATTERN.*)/$1/)
+ {
+ return $line;
+ }
+ if ($line =~ m/Error/ or $line =~ m/Segmentation\sfault/)
+ {
+ return " ". join(" ",$ist->getlines);
+ }
+ if ($line =~ m/Warning/)
+ {
+ $hasWarnings = true;
+ }
+ }
+ $ist->close;
+ if ($cleanUp)
+ {
+ unless ($hasWarnings and $SHOW_WARNINGS)
+ {
+ unlink $file;
+ }
+ }
+ return 0;
+}
+
+sub GetOptionFromConfig
+{
+ my($config, $arg, $aDeprecated) = @ARG;
+
+ if ($config->has($arg))
+ {
+ return $config->get($arg);
+ }
+ if (defined $aDeprecated)
+ {
+ for my $darg (@$aDeprecated)
+ {
+ if ($config->has($darg))
+ {
+ warn "Deprecated option `'$darg`'. Use `'$arg`' instead\n";
+ return $config->get($darg);
+ }
+ }
+ }
+ die "Could not find option $arg";
+}
+
+
+# Main program ---------------------------------
+
+my $R_SOURCE_DIR = $ENV{'BAGPIPE_LIBS'};
+my $OVERWRITE = 0;
+my $memory = 0;
+my $cpus = 1;
+my $outDir = "./";
+my $scanOutDir = "./SCANS/";
+my $nullOutDir = "./NULLSCANS/";
+my $peakOutDir = "./";
+my $plotScale = "cM";
+my $plotScore = "modelcmp";
+my $plotType = "chromosomes";
+my $plotWidth = 12;
+my $plotHeight = 6;
+my %options = ();
+my $tailprobs = "0.001,0.01,0.05,0.1,0.2";
+
+if (0==@ARGV)
+{
+ BasicUsage();
+ exit(1);
+}
+
+my $correctArgList = GetOptions(\%options,
+ 'h',
+ 'help',
+ 'examples',
+ 'example_config',
+ # Directories for input and output
+ 'rsourcedir=s' => \$R_SOURCE_DIR,
+ 'outdir=s' => \$outDir,
+ 'scandir=s' => \$scanOutDir,
+ 'nulldir=s' => \$nullOutDir,
+ 'peakdir=s' => \$peakOutDir,
+ # analysis options
+ 'locus_groups=s',
+ 'overwrite=i' => \$OVERWRITE,
+ # Performance options
+ 'cleanup=i' => \$CLEAN_UP,
+ 'dryrun',
+ 'memory=i' => \$memory,
+ 'cpus=i' => \$cpus,
+ 'warnings=i' => \$SHOW_WARNINGS,
+ 'tailprobs=s' => \$tailprobs,
+ # Actions
+ 'scan',
+ 'scanloci=s',
+ 'perm',
+ 'perm_collate',
+ 'nullsim',
+ 'nullsim_collate',
+ 'nullsimperm',
+ 'nullsimperm_collate',
+ 'posboot_between=s',
+ 'plot',
+ # Plotting cosmetics
+ 'plot_box',
+ 'plot_height=f' => \$plotHeight,
+ 'plot_scale=s' => \$plotScale,
+ 'plot_score=s' => \$plotScore,
+ 'plot_tailprobs=s',
+ 'plot_title=s',
+ 'plot_type=s' => \$plotType,
+ 'plot_width=f' => \$plotWidth,
+ 'qtls=s',
+ # Private options
+ 'devel' => \$DEVEL,
+ );
+
+
+if ($options{'h'})
+{
+ SummaryUsage();
+ exit(0);
+}
+if ($options{'help'})
+{
+ LongUsage();
+ exit(0);
+}
+if ($options{'examples'})
+{
+ PrintExamples();
+ exit(0);
+}
+if ($options{'example_config'})
+{
+ my $configFile = "./config.example";
+ WriteBlankConfigFile(OpenFilesForWriting($configFile));
+ print "Example config file written to $configFile\n";
+ exit(0);
+}
+if (not $correctArgList or 1!=@ARGV)
+{
+ ErrorUsage(1!=@ARGV ? "Must specify configfile" : $ERRNO );
+ exit(1);
+}
+$DRY_RUN = true if $options{'dryrun'};
+
+my $scanDir = "$outDir/$scanOutDir/";
+my $hNullscanProperties = {
+ 'perm' => {
+ 'dir' => "$outDir/$nullOutDir",
+ 'type'=>"permute.scan",
+ 'suffix' => "permscan",
+ },
+ 'nullsim' => {
+ 'dir' => "$outDir/$nullOutDir",
+ 'type'=>"nullsim.scan",
+ 'suffix' => "nullsimscan",
+ },
+ 'nullsimperm' => {
+ 'dir' => "$outDir/$nullOutDir",
+ 'type'=>"nullsimpermute.scan",
+ 'suffix' => "nullsimpermscan",
+ }
+};
+
+#-------------------------------
+# Prescan and update config file
+#-------------------------------
+
+my $configFile = $ARGV[0];
+my $config = new CConfigFile($configFile);
+my $id = $config->get("analysis.id");
+mkpath($outDir, 1, 0755);
+
+# Locus Groups -------------------------------
+
+my $aAvailableLocusGroups = [split m/,/, $config->get("locus.groups")];
+my $aLocusGroups = $aAvailableLocusGroups;
+if (exists $options{'locus_groups'})
+{
+ $aLocusGroups = [split m/,/, $options{'locus_groups'}];
+ unless (Array::IsSubSet($aLocusGroups, $aAvailableLocusGroups))
+ {
+ die "Not all requested locus groups are listed in config file\n";
+ }
+}
+
+# Single locus scans of genome --------------------------
+if (exists $options{'scan'})
+{
+
+ mkpath($scanDir, true, 0755);
+ for my $locusGroup (@$aLocusGroups)
+ {
+ my $scanFile = "$id.group_$locusGroup.scan";
+ my $scanFilePathName = "$scanDir/$scanFile";
+
+ next if FileNotEmpty($scanFilePathName) and not $OVERWRITE;
+
+ SystemR("$R_SOURCE_DIR/scangenome.R --args"
+ ." --configfile=$configFile"
+ ." --locus.group=$locusGroup"
+ ." --outdir=$scanDir"
+ ." --outfile=$scanFile "
+ ." --id=$id",
+ "$outDir/$scanFile.Rout",
+ "genome scan of locus group $locusGroup");
+
+ if (FileEmptyOrAbsent($scanFilePathName) and not $DRY_RUN)
+ {
+ die "Genome scan failed for $locusGroup\n";
+ }
+ }
+}
+
+# Single locus scans of individual loci ------------------
+if (exists $options{'scanloci'})
+{
+ my $aSloci= [split m/,/, $options{'scanloci'}];
+ my $aSfiles= [map{$id."_".$ARG.".save"} @$aSloci];
+
+ SystemR("$R_SOURCE_DIR/scanloci.R --args"
+ ." --configfile=$configFile"
+ ." --outdir=$outDir"
+ ." --loci=".join(',', @$aSloci)
+ ." --loci.files=".join(',', @$aSfiles)
+ ." --id=$id",
+ "$outDir/scanloci.$PID.Rout",
+ "scanning individual loci with detailed output");
+ # todo check for output files
+}
+
+# Null scans -----------------------------------
+
+my $aNullScans = [];
+push(@$aNullScans, 'perm') if exists $options{'perm'};
+push(@$aNullScans, 'nullsim') if exists $options{'nullsim'};
+push(@$aNullScans, 'nullsimperm') if exists $options{'nullsimperm'};
+for my $nullscan (@$aNullScans)
+{
+ my $dir = $hNullscanProperties->{$nullscan}{'dir'};
+ mkpath($dir, true, 0755);
+ my $suffix = $hNullscanProperties->{$nullscan}{'suffix'};
+ my $aScanFiles = [];
+ for my $locusGroup (@$aLocusGroups)
+ {
+ my $scanFile = "$id.group_$locusGroup.$suffix";
+ push @$aScanFiles, $scanFile;
+
+ next if FileNotEmpty("$dir/$scanFile") and not $OVERWRITE;
+
+ SystemR("$R_SOURCE_DIR/nullscangenome.R --args"
+ ." --configfile=$configFile"
+ ." --cpus=$cpus"
+ ." --id=$id"
+ ." --scan.type=".$hNullscanProperties->{$nullscan}{'type'}
+ ." --outdir=$dir"
+ ." --outfile=$scanFile"
+ ." --locus.group=$locusGroup"
+ ." --memory=$memory",
+ "$outDir/$scanFile.Rout",
+ "Error in $nullscan scan of locus group $locusGroup");
+
+ if (FileEmptyOrAbsent("$dir/$scanFile") and not $DRY_RUN)
+ {
+ die "$nullscan failed for $locusGroup\n";
+ }
+ }
+}
+
+my $aNullScanCollations = [];
+push(@$aNullScanCollations, 'perm') if exists $options{'perm_collate'};
+push(@$aNullScanCollations, 'nullsim') if exists $options{'nullsim_collate'};
+push(@$aNullScanCollations, 'nullsimperm') if exists $options{'nullsimperm_collate'};
+for my $nullscan (@$aNullScanCollations)
+{
+ my $dir = $hNullscanProperties->{$nullscan}{'dir'};
+ my $suffix = $hNullscanProperties->{$nullscan}{'suffix'};
+ my $aScanFiles = [];
+ for my $locusGroup (@$aAvailableLocusGroups)
+ {
+ my $scanFile = "$id.group_$locusGroup.$suffix";
+ if (FileEmptyOrAbsent("$dir/$scanFile"))
+ {
+ die "Cannot collate $nullscan output: missing results for locus group $locusGroup\n";
+ }
+ push @$aScanFiles, $scanFile;
+ }
+ my $gevFile = "$outDir/$id.gev.$suffix";
+ my $threshFile = "$outDir/$id.thresh.$suffix";
+ SystemR("$R_SOURCE_DIR/fitgev.R --args"
+ ." --configfile=$configFile"
+ ." --infiles=".join(",", @$aScanFiles)
+ ." --indir=$dir"
+ ." --threshfile=$threshFile"
+ ." --gevfile=$gevFile"
+ ." --tailprobs=$tailprobs",
+ "$threshFile.Rout",
+ "Error in GEV fit of permutation scan");
+}
+
+# Plotting -----------------------------------
+if (exists $options{'plot'})
+{
+ my $aThresholds = [];
+ if (exists $options{'plot_tailprobs'})
+ {
+ for my $string (split(m/,/, $options{'plot_tailprobs'}))
+ {
+ my($nullscan, $tailprob) = split(m/_/, $string);
+ unless (exists $hNullscanProperties->{$nullscan})
+ {
+ die "Could not recognize nullscan type $nullscan in plot_tailprob argument\n";
+ }
+ my $threshFile = "$outDir/$id.thresh.".$hNullscanProperties->{$nullscan}{'suffix'};
+ my $df = new CDataFrame(
+ 'header'=> true,
+ 'stream'=> OpenFilesForReading($threshFile),
+ 'format'=> 'delim');
+ my $p = $df->getCol("upper.tail.prob");
+ my $quant = undef;
+ for (my $i=0; $i<$df->getNumRows; $i++)
+ {
+ if ($tailprob==$p->[$i])
+ {
+ $quant = $df->getCell($i,$plotScore);
+ last;
+ }
+ }
+ unless (defined $quant)
+ {
+ die "Could not find threshold for $plotScore tail prob $tailprob in $threshFile\n";
+ }
+ push @$aThresholds, $quant;
+ }
+ }
+
+ my $outFile = "$outDir/$id.scans.pdf";
+ my $aScanFiles = [];
+ for my $locusGroup (@$aAvailableLocusGroups)
+ {
+ my $scanFile = "$id.group_$locusGroup.scan";
+ if (FileEmptyOrAbsent("$scanDir/$scanFile"))
+ {
+ warn "Skipping missing data for locus group $locusGroup\n";
+ next;
+ }
+ push @$aScanFiles, $scanFile;
+ }
+ my $scanFileList = join ',', @$aScanFiles;
+ SystemR("$R_SOURCE_DIR/plotscans.R --args"
+ ." --configfile=$configFile"
+ ." --plot.height=$plotHeight"
+ ." --plot.width=$plotWidth"
+ ." --plot.type=$plotType"
+ .($options{'plot_box'} ? " --plot.box" : "")
+ ." --indir=$scanDir/"
+ ." --infiles=$scanFileList"
+ .($options{'plot_title'} ? " --title=".$options{'plot_title'} : "")
+ ." --plot.scale=$plotScale"
+ ." --locus.score.type=$plotScore"
+ ." --outfile=$outFile"
+ .(@$aThresholds ? " --thresholds=".join(",", @$aThresholds) : " ")
+ ." --flag.missing",
+ "$outDir/$id.scans.pdf.Rout",
+ "Error in plot");
+}
+
+# Positional bootstrap -------------------------------
+
+if (exists $options{'posboot_between'})
+{
+ my($locus1, $locus2, $numBoots) = split m/,/, $options{'posboot_between'};
+ my $tag="$id.$locus1-$locus2";
+ my $posbootResults="$outDir/$tag.posboot";
+ my $posbootSummary="$outDir/$tag.posbootsum";
+ my $posbootQuickSummary="$outDir/$tag.posbootqsum";
+ my $posbootPlot="$outDir/$tag.posboot.pdf";
+
+ SystemR("$R_SOURCE_DIR/posboot.R --args"
+ ." --configfile=$configFile"
+ ." --memory=$memory"
+ ." --num.posboots=$numBoots"
+ ." --posboot.between=$locus1,$locus2"
+ ." --posboot.results=$posbootResults"
+ ." --posboot.summary=$posbootSummary"
+ ." --posboot.quicksummary=$posbootQuickSummary"
+ ." --posboot.plot=$posbootPlot"
+ ." --plot.scale=$plotScale"
+ ." --id=$id",
+ "$outDir/$posbootResults.Rout",
+ "bootstrapping $numBoots times between $locus1 and $locus2");
+}
+
+
+
+
+
diff --git a/bagpipe/lib/Algorithm.pm b/bagpipe/lib/Algorithm.pm
new file mode 100644
index 0000000..92780f9
--- /dev/null
+++ b/bagpipe/lib/Algorithm.pm
@@ -0,0 +1,467 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+#==============================================================================
+# MODULE: Algorithm
+#======================================
+# DESCRIPTION:
+#
+
+#_ Module declaration
+
+package Algorithm;
+
+#_ Include libraries
+
+use English;
+use StdDefs;
+use Assert;
+use Array;
+use Exporter;
+use VarType qw(IsArrayRef IsHashRef);
+use base qw(Exporter);
+use vars qw( @EXPORT_OK );
+
+ at EXPORT_OK = qw(
+ Accumulate
+ CopyIf
+ Count
+ CountIf
+ CountString
+ FindIf
+ ForEach
+ Generate
+ MaxElement
+ MinElement
+ RandomShuffle
+ Transform
+ );
+
+#------------------------------------------------------------------------------
+# INTERFACE
+#----------
+
+#_ Public constants
+
+use constant NOT_FOUND => null;
+
+#_ Public functions
+
+sub Accumulate;
+# Function:
+# Accumulate the results of an operation on an array.
+# ARGUMENTS:
+# EITHER
+# 1. <\@> input
+# 2. <real> initial value
+# OR
+# 1. <\@> input
+# 2. <?> initial value
+# 3. <\sub> binary operation that "adds" two elements
+# RETURN:
+# 1. <real> or <?> accumulated value
+
+sub CopyIf($$$);
+# Function:
+# Copy elements that match a predicate from one array into another array
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <\@> output
+# 3. <\sub> predicate
+# RETURN:
+# 1. <\@> output
+
+sub Count($$);
+# Function:
+# Count occurences of a value in an array. Uses the '==' operator to assess
+# equivalence.
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <SCALAR> value
+# RETURN:
+# 1. <int> count
+
+sub CountIf($$);
+# Function:
+# Count matches of a predicate in an array.
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <\sub> predicate
+# RETURN:
+# 1. <int> count
+
+sub CountString($$);
+# Function:
+# Count occurences of a value in an array. Uses the 'eq' operator to assess
+# equivalence.
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <SCALAR> value
+# RETURN:
+# 1. <int> count
+
+sub FindIf($$);
+# Function:
+# Find first match of a predicate in an array.
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <\sub> predicate
+# RETURN:
+# 1. <?> element if found, <CONSTANT> NOT_FOUND otherwise.
+
+sub ForEach($$);
+# Function:
+# Do operation for every element in an array.
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <\sub> operation
+# RETURN:
+# 1. <\sub> operation
+
+sub Generate($$$);
+
+sub Transform($$$);
+# Function:
+# Apply an operation to every element in an array
+# ARGUMENTS:
+# 1. <\@> input
+# 2. <\@> output (may be the same as input)
+# 3. <\sub> operation
+# RETURN:
+# 1. <\@> output
+# SEE ALSO:
+
+
+
+#------------------------------------------------------------------------------
+# IMPLEMENTATION
+#---------------
+
+sub Accumulate
+{
+ if (2 == @ARG)
+ {
+ return _Accumulate_2Args(@ARG);
+ }
+ elsif (3 == @ARG)
+ {
+ return _Accumulate_3Args(@ARG);
+ }
+ WrongNumArgsError();
+}
+
+sub _Accumulate_2Args($$)
+{
+ my ($array, $init) = @ARG;
+ for my $elem (@$array)
+ {
+ $init += $elem;
+ }
+ return $init;
+}
+
+sub _Accumulate_3Args($$$)
+{
+ my($array, $init, $binOp) = @ARG;
+ for my $elem (@$array)
+ {
+ $init += &$binOp($init, $elem);
+ }
+ return $init;
+}
+
+sub CopyIf($$$)
+{
+ my($aIn, $aOut, $unOp) = @ARG;
+
+ for my $elem (@$aIn)
+ {
+ if (&$unOp($elem))
+ {
+ push @$aOut, $elem;
+ }
+ }
+ return $aOut;
+}
+
+sub Count($$)
+{
+ my($aIn, $value) = @ARG;
+
+ my $count = 0;
+ for my $elem (@$aIn)
+ {
+ $count++ if $elem == $value;
+ }
+ return $count;
+}
+
+sub CountIf($$)
+{
+ my($aIn, $op) = @ARG;
+
+ my $count = 0;
+ for my $elem (@$aIn)
+ {
+ if (&$op($elem))
+ {
+ $count++;
+ }
+ }
+ return $count;
+}
+
+sub CountString($$)
+{
+ my($aIn, $value) = @ARG;
+
+ my $count = 0;
+ for my $elem (@$aIn)
+ {
+ $count++ if $elem eq $value;
+ }
+ return $count;
+}
+
+sub FindIf($$)
+{
+ my($array, $pred) = @ARG;
+ for my $elem (@$array)
+ {
+ if (&$pred($elem))
+ {
+ return $elem;
+ }
+ }
+ return NOT_FOUND;
+}
+
+sub ForEach($$)
+{
+ my($array, $op) = @ARG;
+ for my $elem (@$array)
+ {
+ &$op($elem);
+ }
+ return $op;
+}
+
+sub Generate($$$)
+{
+ my($array, $nElem, $func) = @ARG;
+
+ return Array::Generate($array, $nElem, $func);
+}
+
+sub MaxElement
+{
+ 1 == @ARG || 2 == @ARG or WrongNumArgsError();
+
+ my $aIn = $ARG[0];
+ my $binPred = (2 == @ARG) ? $ARG[1]
+ : sub{$ARG[0] <=> $ARG[1]};
+ my $max = $aIn->[0];
+ for (my $i=1; $i<@$aIn; $i++)
+ {
+ if (-1 == &$binPred($max, $aIn->[$i]))
+ {
+ $max = $aIn->[$i];
+ }
+ }
+ return $max;
+}
+
+sub MinElement
+{
+ 1 == @ARG || 2 == @ARG or WrongNumArgsError();
+
+ my $aIn = $ARG[0];
+ my $binPred = (2 == @ARG) ? $ARG[1]
+ : sub{$ARG[0] <=> $ARG[1]};
+ my $min = $aIn->[0];
+ for (my $i=1; $i<@$aIn; $i++)
+ {
+ if (1 == &$binPred($min, $aIn->[$i]))
+ {
+ $min = $aIn->[$i];
+ }
+ }
+ return $min;
+}
+
+sub RandomShuffle($)
+{
+ return Array::Shuffle($ARG[0]);
+}
+
+sub Transform($$$)
+{
+ my($inputCont, $outputCont, $rsub) = @ARG;
+
+ if (IsArrayRef($inputCont) && IsArrayRef($outputCont))
+ {
+ return Transform_array($inputCont, $outputCont, $rsub);
+ }
+ elsif (IsHashRef($inputCont) && IsHashRef($outputCont))
+ {
+ return Transform_hash($inputCont, $outputCont, $rsub);
+ }
+ else
+ {
+ Fatal("non-container type passed to function");
+ }
+}
+
+sub Transform_hash($$$)
+{
+ my($hIn, $hOut, $op) = @ARG;
+
+ while (my($key, $val) = each %$hIn)
+ {
+ $hOut->{$key} = &$op($val);
+ }
+ return $hOut;
+}
+
+sub Transform_array($$$)
+{
+ my($aIn, $aOut, $op) = @ARG;
+
+ return Array::Transform($aIn, $aOut, $op);
+}
+
+#====
+# END
+#==============================================================================
+true;
+__END__
+
+=head1 Module
+
+Algorithm
+
+=head1 Description
+
+Loops can be more confusing than they're worth, especially when they perform trivial tasks. This module is loosely based on the C++ algorithm.h module, which aims to cut down the number of explicit loops in code.
+
+=head1 Notes
+
+=over
+
+=item predicate
+
+This is a function that returns nonzero iff its argument satisfy a condition. If it take 1 argument, it is a unary predicate, 2 arguments, a binary predicate, and so on.
+
+=item operation
+
+This is a function that returns a value. If it takes 1 argument, it is unary, 2 args then binary, etc.
+
+=back
+
+=head1 Functions
+
+=over
+
+=item Accumulate
+
+ Function:
+ Accumulate the results of an operation on an array.
+ ARGUMENTS:
+ EITHER
+ 1. <\@> input
+ 2. <real> initial value
+ OR
+ 1. <\@> input
+ 2. <?> initial value
+ 3. <\sub> binary operation that "adds" two elements
+ RETURN:
+ 1. <real> or <?> accumulated value
+
+=item CopyIf
+
+ Function:
+ Copy elements that match a predicate from one array into another array
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <\@> output
+ 3. <\sub> predicate
+ RETURN:
+ 1. <\@> output
+
+=item Count
+
+ Function:
+ Count occurences of a value in an array. Uses the '==' operator to assess
+ equivalence.
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <SCALAR> value
+ RETURN:
+ 1. <int> count
+
+=item CountIf
+
+ Function:
+ Count matches of a predicate in an array.
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <\sub> predicate
+ RETURN:
+ 1. <int> count
+
+=item CountString
+
+ Function:
+ Count occurences of a value in an array. Uses the 'eq' operator to assess
+ equivalence.
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <SCALAR> value
+ RETURN:
+ 1. <int> count
+
+=item FindIf
+
+ Function:
+ Find first match of a predicate in an array.
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <\sub> predicate
+ RETURN:
+ 1. <?> element if found, <CONSTANT> NOT_FOUND otherwise.
+
+=item ForEach
+
+ Function:
+ Do operation for every element in an array.
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <\sub> operation
+ RETURN:
+ 1. <\sub> operation
+
+=item Transform
+
+ Function:
+ Apply an operation to every element in an array
+ ARGUMENTS:
+ 1. <\@> input
+ 2. <\@> output (may be the same as input)
+ 3. <\sub> operation
+ RETURN:
+ 1. <\@> output
+
diff --git a/bagpipe/lib/Array.pm b/bagpipe/lib/Array.pm
new file mode 100644
index 0000000..ca04649
--- /dev/null
+++ b/bagpipe/lib/Array.pm
@@ -0,0 +1,1111 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package Array;
+use strict;
+use English;
+use StdDefs;
+use Assert;
+use Hash;
+use CRing;
+use VarType;
+use Debug;
+
+use constant NOT_FOUND => -1;
+
+my $RECORD_START = "#begin#Array\n";
+my $RECORD_END = "#end#Array\n";
+
+sub Array::Add
+{
+ my @arrayList = @ARG;
+
+ WrongNumArgsError() unless @ARG;
+
+ my $nElem = scalar(@{$arrayList[0]});
+ my $aAdded = Array::New($nElem, 0);
+ for my $array (@arrayList)
+ {
+ if (VarType::IsArrayRef($array))
+ {
+ for (my $i=0; $i<@$aAdded; $i++)
+ {
+ $aAdded->[$i] += $array->[$i];
+ }
+ }
+ else
+ {
+ for my $elem (@$aAdded)
+ {
+ $elem += $array;
+ }
+ }
+ }
+ return $aAdded;
+}
+
+sub Array::All
+{
+ my $array = $ARG[0];
+ my $predFunc = (2==@ARG) ? $ARG[1]
+ : 0;
+ for my $elem (@$array)
+ {
+ if ($predFunc)
+ {
+ if (not $predFunc->($elem))
+ {
+ return false;
+ }
+ }
+ else
+ {
+ return false if not $elem;
+ }
+ }
+ return true;
+}
+
+sub Array::Any
+{
+ my $array = $ARG[0];
+ my $predFunc = (2==@ARG) ? $ARG[1]
+ : 0;
+ for my $elem (@$array)
+ {
+ if ($predFunc)
+ {
+ if ($predFunc->($elem))
+ {
+ return true;
+ }
+ }
+ else
+ {
+ return true if $elem;
+ }
+ }
+ return false;
+}
+
+sub Array::Apply($$)
+{
+ my($array, $func) = @ARG;
+
+ return Array::Transform($array, [], $func);
+}
+
+sub Array::ApplyInSitu($$)
+{
+ my($array, $func) = @ARG;
+
+ for my $elem (@$array)
+ {
+ $elem = $func->($elem);
+ }
+ return $array;
+}
+
+
+sub Array::Back($)
+{
+ my $array = shift;
+ return $array->[scalar(@$array)-1];
+}
+
+sub Array::Common
+{
+ my @arrayList = @ARG;
+
+ my $flatUniqArrays = [];
+ for my $array (@arrayList)
+ {
+ push @$flatUniqArrays, @{ Unique($array) };
+ }
+ my $hElemFreq = ElementFrequencies($flatUniqArrays);
+
+ my $aCommon = [];
+ while( my($elem, $elemFreq) = each %$hElemFreq )
+ {
+ if( scalar(@arrayList) == $elemFreq )
+ {
+ push @$aCommon, $elem;
+ }
+ }
+ return $aCommon;
+}
+
+sub Array::ConditionToMask($$)
+{
+ my($array, $predicate) = @ARG;
+
+ my $mask = [];
+ for (my $i=0; $i<@$array; $i++)
+ {
+ if ($predicate->($array->[$i]))
+ {
+ push @$mask, true;
+ }
+ else
+ {
+ push @$mask, false;
+ }
+ }
+ return $mask;
+}
+
+sub Array::ContainsString($$)
+{
+ my($aIn, $value) = @ARG;
+ for my $elem (@$aIn)
+ {
+ return true if $elem eq $value;
+ }
+ return false;
+}
+
+
+sub Copy($)
+{
+ my $array = shift;
+ return [@$array];
+}
+
+sub CopyRange($$$)
+{
+ my($array, $start, $end) = @ARG;
+
+ my $copy = [];
+ for (my $i=$start; $i<=$end; $i++)
+ {
+ push @$copy, $array->[$i];
+ }
+ return $copy;
+}
+
+sub Array::CountUnique
+{
+ return scalar(@{Array::Unique(@ARG)});
+}
+
+sub Diff($$)
+{
+ my( $array1, $array2 ) = @ARG;
+
+ my $aArray1Not2 = [];
+ my $aArray2Not1 = [];
+ my $aIntersect = [];
+
+ my %hArray1UniqueElem = ();
+ my %hArray2UniqueElem = ();
+ my $hUnionElemFreq;
+
+ my $elem;
+ my $elemFreq;
+
+ #---------------------------------
+ #- get unique elements of each set
+ #---------------------------------
+
+ @hArray1UniqueElem{ @$array1 } = ();
+ @hArray2UniqueElem{ @$array2 } = ();
+
+ #--------------------------------------------------
+ #- get the union of set1(unique) u set2(unique)
+ #- make a frequency table of elements in this union
+ #--------------------------------------------------
+
+ $hUnionElemFreq = ElementFrequencies([ keys %hArray1UniqueElem, keys %hArray2UniqueElem ]);
+
+ while( ($elem, $elemFreq) = each %$hUnionElemFreq )
+ {
+ #----------------------------------------------------
+ #- non-intersection elements will appear only once...
+ #----------------------------------------------------
+ if( 1 == $elemFreq )
+ {
+ if(exists $hArray1UniqueElem{ $elem } )
+ {
+ push @$aArray1Not2, $elem;
+ }
+ else
+ {
+ push @$aArray2Not1, $elem;
+ }
+ }
+ #--------------------------------------------------
+ #- ...while intersection elements will appear twice
+ #--------------------------------------------------
+ else
+ {
+ push @$aIntersect, $elem;
+ }
+ }
+
+ Prefer( wantarray, "Client should accept array from method" );
+
+ return ( $aArray1Not2, $aArray2Not1, $aIntersect );
+}
+
+# forward decl
+sub Array::Dimensionality($);
+#
+
+sub Array::Dimensionality($)
+{
+ my $array = shift;
+
+ if (VarType::IsArrayRef($array))
+ {
+ return 1 + Array::Dimensionality($array->[0]);
+ }
+ return 0;
+}
+
+sub Array::Divide($$$)
+{
+ my($aIn, $aOut, $divisor) = @ARG;
+
+ if (VarType::IsArrayRef($divisor))
+ {
+ Fatal("Divisor array must be same length as target") if @$aIn!=@$divisor;
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]/$divisor->[$i];
+ }
+ }
+ elsif ($aIn==$aOut)
+ {
+ for my $elem (@$aOut) {$elem /= $divisor;}
+ }
+ else
+ {
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]/$divisor;
+ }
+ }
+ return $aOut;
+}
+
+sub ElementFrequencies($)
+{
+ my ($array) = @ARG;
+ my $hFreq = {};
+
+ for my $elem ( @$array )
+ {
+ if( not exists($hFreq->{$elem}) )
+ {
+ $hFreq->{$elem} = 1;
+ }
+ else
+ {
+ $hFreq->{$elem}++;
+ }
+ }
+ return $hFreq;
+}
+
+sub ElementRelativeFreqs($)
+{
+ my $array = shift;
+
+ my $hFreqs = ElementFrequencies($array);
+
+ for my $key (keys %$hFreqs)
+ {
+ $hFreqs->{$key} = $hFreqs->{$key} / scalar(@$array);
+ }
+ return $hFreqs;
+}
+
+sub Array::EqualContents($$)
+{
+ my($arrayA, $arrayB) = @ARG;
+
+ my $hFreqA = Array::ElementFrequencies($arrayA);
+ my $hFreqB = Array::ElementFrequencies($arrayB);
+
+ return Hash::Equals($hFreqA, $hFreqB);
+}
+
+sub Array::EqualUnique($$)
+{
+ my($arrayA, $arrayB) = @ARG;
+
+ my $uniqA = Array::Unique($arrayA);
+ my %seenB = ();
+ @seenB{@$arrayB} = ();
+
+ return false unless @$uniqA == scalar keys %seenB;
+
+ for my $key (@$uniqA)
+ {
+ return false unless exists $seenB{$key};
+ }
+ return true;
+}
+
+sub Equals($$)
+{
+ my($arrayA, $arrayB) = @ARG;
+
+ # compare sizes
+ return false if scalar(@$arrayA) != scalar(@$arrayB);
+
+ # compare elements
+ for (my $i=0; $i<@$arrayA; $i++)
+ {
+ return false if $arrayA->[$i] ne $arrayB->[$i];
+ }
+ return true;
+}
+
+sub Exclude($$)
+{
+ my($array, $not) = @ARG;
+ my $hNot = Hash::New(VarType::AsArray($not), true);
+
+ my $out = [];
+ for my $elem (@$array)
+ {
+ unless (exists $hNot->{$elem})
+ {
+ push @$out, $elem;
+ }
+ }
+ return $out;
+}
+
+sub Array::FindFirst($$)
+{
+ my($array, $target) = @ARG;
+
+ for (my $i=0; $i<@$array; $i++)
+ {
+ return $i if $array->[$i] eq $target;
+ }
+ return NOT_FOUND;
+}
+
+sub Array::FindFirstIf($$)
+{
+ my($array, $func) = @ARG;
+
+ for (my $i=0; $i<@$array; $i++)
+ {
+ return $i if $func->($array->[$i]);
+ }
+ return NOT_FOUND;
+}
+
+sub Generate($$$)
+{
+ my($array, $nElem, $func) = @ARG;
+
+ for (my $i=0; $i<$nElem; $i++)
+ {
+ $array->[$i] = &$func;
+ }
+ return $array;
+}
+
+
+sub GroupElementsBy($$)
+{
+ my($array, $groupingFunc) = @ARG;
+
+ my $hm = {}; # hash of arrays
+ for my $elem (@$array)
+ {
+ my $key = $groupingFunc->($elem);
+ unless (exists $hm->{$key})
+ {
+ $hm->{$key} = [];
+ }
+ push @{$hm->{$key}}, $elem;
+ }
+ return [values %$hm];
+}
+
+sub HashByIndex($)
+{
+ my $array = shift;
+
+ warn "deprecated method\n";
+
+ return HashElement2Index($array);
+}
+
+sub HashElement2Index($)
+{
+ my $array = shift;
+ my $hash = {};
+
+ my $i = 0;
+ for my $elem (@$array)
+ {
+ $hash->{$elem} = $i++;
+ }
+ return $hash;
+}
+
+sub HashIndex2Element($)
+{
+ my $array = shift;
+ my $hash = {};
+ for (my $i=0; $i<@$array; $i++)
+ {
+ $hash->{$i} = $array->[$i];
+ }
+ return $hash;
+}
+
+sub HasIntersect
+{
+ my $aCommon = Array::Common(@ARG);
+ return 0 != @$aCommon
+}
+
+sub IndicesToMask($$)
+{
+ my($nElem, $aIndices) = @ARG;
+
+ my $mask = Array::New($nElem, false);
+ for my $index (@$aIndices)
+ {
+ $mask->[$index] = true;
+ }
+ return $mask;
+}
+
+sub Array::InvertMask($)
+{
+ my $mask = shift;
+ return Array::Transform($mask, [], sub{$ARG[0]?false:true});
+}
+
+sub IsSubSet($$)
+{
+ my($aSmall, $aBig) = @ARG;
+
+ my $hBig = Array::ToHash($aBig, sub{$ARG[0]}, sub{true});
+
+ for my $elem (@$aSmall)
+ {
+ if (not exists $hBig->{$elem})
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
+
+sub Join
+{
+ my(@aArrays) = @ARG;
+
+ my $aJoined = [];
+ for my $array (@aArrays)
+ {
+ push @$aJoined, @$array;
+ }
+ return $aJoined;
+}
+
+sub Array::LastIndex($)
+{
+ return scalar(@{$ARG[0]})-1;
+}
+
+sub Length($)
+{
+ return scalar(@{$ARG[0]});
+}
+
+sub Load($)
+{
+ my $ist = shift;
+
+ # skip to record start
+ return undef if $ist->eof;
+ while ($ist->getline ne $RECORD_START)
+ {
+ return undef if $ist->eof;
+ }
+
+ my $array = [];
+ while ((my $line = $ist->getline) ne $RECORD_END)
+ {
+ if ($ist->eof)
+ {
+ Fatal("Premature EOF encountered when reading object from file");
+ }
+ $line =~ m/^(\S+)\s+(\S.*)$/;
+ $array->[$1] = $2;
+ }
+ return $array;
+}
+
+sub Array::Mask($$)
+{
+ my($array, $mask) = @ARG;
+
+ if (@$array != @$mask)
+ {
+ Fatal("mask must be as long as array");
+ }
+
+ my $newArray = [];
+ for (my $i=0; $i<@$array; $i++)
+ {
+ if ($mask->[$i])
+ {
+ push @$newArray, $array->[$i];
+ }
+ }
+ return $newArray;
+}
+
+sub Array::Match($$)
+{
+ my($aAll, $aTargets) = @ARG;
+
+ my $hTarget2FoundIndex = Hash::New($aTargets, NOT_FOUND);
+
+ my $i=0;
+ for my $elem (@$aAll)
+ {
+ if (exists $hTarget2FoundIndex->{$elem}
+ and
+ $hTarget2FoundIndex->{$elem} eq NOT_FOUND)
+ {
+ $hTarget2FoundIndex->{$elem} = $i;
+ }
+ $i++
+ }
+ my $aIndices = [];
+ for my $target (@$aTargets)
+ {
+ push @$aIndices, $hTarget2FoundIndex->{$target};
+ }
+ return $aIndices;
+}
+
+sub Max($)
+{
+ my ($array) = @ARG;
+
+ my $elem;
+ my $max;
+ my $firstIteration = true;
+
+ for $elem (@$array)
+ {
+ if($firstIteration)
+ {
+ $max = $elem;
+ $firstIteration = false;
+ }
+ elsif($elem > $max)
+ {
+ $max = $elem;
+ }
+ }
+ return $max;
+}
+
+sub Array::Mean($)
+{
+ my $array = shift;
+
+ return Array::Sum($array) / scalar(@$array);
+}
+
+sub Min($)
+{
+ my ($array) = @ARG;
+
+ my $elem;
+ my $min;
+ my $firstIteration = true;
+
+ for $elem (@$array)
+ {
+ if($firstIteration)
+ {
+ $min = $elem;
+ $firstIteration = false;
+ }
+ elsif($elem < $min)
+ {
+ $min = $elem;
+ }
+ }
+ return $min;
+}
+
+sub ModalFreq($)
+{
+ my $array = shift;
+
+ return Max([values %{ElementFrequencies($array)}]);
+}
+
+sub Modes($)
+{
+ my $array = shift;
+
+ my $hFreq = Array::ElementFrequencies($array);
+ my $modeFreq = Array::Max([values %$hFreq]);
+ my @aModes = ();
+ while (my($elem, $freq) = each %$hFreq)
+ {
+ push @aModes, $elem if $modeFreq == $freq;
+ }
+ return @aModes;
+}
+
+sub Array::Multiply($$$)
+{
+ my($aIn, $aOut, $factor) = @ARG;
+
+ if (VarType::IsArrayRef($factor))
+ {
+ Fatal("Divisor array must be same length as target") if @$aIn!=@$factor;
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]*$factor->[$i];
+ }
+ }
+ elsif ($aIn==$aOut)
+ {
+ for my $elem (@$aOut) { $elem *= $factor };
+ }
+ else
+ {
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]*$factor;
+ }
+ }
+ return $aOut;
+}
+
+sub New
+#_ either NewArray( size ) or
+#_ NewArray( size, value )
+#_ if size == 0, will confess with WrongNumArgsError
+{
+ my $size = $ARG[0] || WrongNumArgsError();
+ my $array;
+
+ if(2 == @ARG)
+ {
+ $array = [($ARG[1]) x $size ];
+ }
+ else
+ {
+ $array = [];
+ $#{$array} = $size - 1;
+ }
+ return $array;
+}
+
+sub NotEmpty($)
+{
+ return scalar(@{$ARG[0]});
+}
+
+sub NumElements
+{
+ return scalar(@{$ARG[0]});
+}
+
+sub NumUnique
+{
+ my @arg = @ARG;
+ Warn("Deprecated method 'NumUnique'");
+ return CountUnique(@arg);
+}
+
+sub Array::NextPermutation($)
+# based on the C++ STL's next.permutation()
+{
+ my $array = shift;
+ if (2>@$array) {return false;}
+
+ my $last = Array::Length($array) - 1;
+ my $i = $last;
+ for (;;)
+ {
+ my $ii = $i--;
+ if ($array->[$i] lt $array->[$ii])
+ {
+ my $j = $last+1;
+ while (not ($array->[$i] lt $array->[--$j])){}
+ Array::SwapElements($array, $i, $j);
+ Array::ReverseRange($array, $ii, $last);
+ return true;
+ }
+ if (0==$i)
+ {
+ Array::ReverseRange($array, 0, $last);
+ return false;
+ }
+ }
+}
+
+sub Paste
+{
+ my $aIns = [];
+ my $aOptions = [];
+
+ for my $arg (@ARG)
+ {
+ if (VarType::IsArrayRef($arg))
+ {
+ push @$aIns, new CRing($arg);
+ }
+ else
+ {
+ push @$aOptions, $arg;
+ }
+ }
+ my %hArg = @$aOptions;
+
+ Assert::CheckArgHash(\%hArg,
+ 'optional' => [qw(sep collapse)],
+ );
+
+ my $sep = exists $hArg{'sep'}
+ ? $hArg{'sep'}
+ : " ";
+
+ my $numItems = Array::Max(
+ Array::Transform($aIns, [], sub{$ARG[0]->getSize})
+ );
+
+ my $out = [];
+ for (my $i=0; $i<$numItems; $i++)
+ {
+ my $string = "";
+ for (my $r=0; $r<@$aIns; $r++)
+ {
+ $string .= $sep unless 0==$r;
+ $string .= $aIns->[$r]->nextElement;
+ }
+ push @$out, $string;
+ }
+
+ if (exists $hArg{'collapse'})
+ {
+ return join($hArg{'collapse'}, @$out);
+ }
+ return $out;
+}
+
+sub Array::Plus($$$)
+{
+ my($aIn, $aOut, $plus) = @ARG;
+
+ if (VarType::IsArrayRef($plus))
+ {
+ Fatal("Divisor array must be same length as target") if @$aIn!=@$plus;
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]+$plus->[$i];
+ }
+ }
+ elsif ($aIn==$aOut)
+ {
+ for my $elem (@$aOut) {$elem += $plus;}
+ }
+ else
+ {
+ for (my $i=0; $i<@$aIn; ++$i)
+ {
+ $aOut->[$i] = $aIn->[$i]+$plus;
+ }
+ }
+ return $aOut;
+}
+
+sub RandomElement($)
+{
+ my $array = shift;
+ return $array->[int rand scalar @$array];
+}
+
+sub Range($)
+{
+ my ($array) = @ARG;
+
+ my $max;
+ my $min;
+ my $firstIteration = true;
+
+ for my $elem (@$array)
+ {
+ if($firstIteration)
+ {
+ $max = $elem;
+ $min = $elem;
+ $firstIteration = false;
+ }
+ elsif($elem < $min)
+ {
+ $min = $elem;
+ }
+ elsif($elem > $max)
+ {
+ $max = $elem;
+ }
+ }
+ my $range = abs( $max - $min );
+
+ return wantarray ?
+ ($range, $min, $max):
+ $range;
+}
+
+sub Array::Rep
+# modeled after the S-plus function rep
+# Will ignore 'length' if 'times' is given
+{
+ my($in, %hArg) = @ARG;
+
+ $in = VarType::AsArray($in);
+
+ Assert::CheckArgHash(\%hArg,
+ 'optional' => [qw(times length each)],
+ );
+ unless (exists $hArg{'times'} or exists $hArg{'length'})
+ {
+ Error::Fatal("Method must have one of times or length as args");
+ }
+ my $each = exists $hArg{'each'}
+ ? $hArg{'each'}
+ : 1;
+ my $length = exists $hArg{'times'}
+ ? scalar(@$in)*$hArg{'times'}*$each
+ : $hArg{'length'};
+ my $out = [];
+ while(true)
+ {
+ for my $elem (@$in)
+ {
+ for (my $i=0; $i<$each; $i++)
+ {
+ push @$out, $elem;
+ return $out if $length == @$out;
+ }
+ }
+ }
+}
+
+sub Resize($$)
+{
+ my($array, $newSize) = @ARG;
+ $#{$array} = $newSize - 1;
+}
+
+sub Array::ReverseRange($$$)
+{
+ my($array, $first, $last) = @ARG;
+ @{$array}[$first..$last] = reverse @{$array}[$first..$last];
+ return $array;
+}
+
+
+sub Shuffle($)
+{
+ my $array = shift;
+
+ my $nElem = @$array;
+
+ for (my $i=0; $i<$nElem; $i++)
+ {
+ my $j = int rand($nElem);
+
+ next if $i==$j;
+
+ # swap i and j
+ my $tmp = $array->[$i];
+ $array->[$i] = $array->[$j];
+ $array->[$j] = $tmp;
+ }
+ return $array;
+}
+
+
+sub Save($$)
+{
+ my($ost, $array) = @ARG;
+
+ $ost->print($RECORD_START);
+ my $i = 0;
+ for my $elem (@$array)
+ {
+ $ost->print("$i\t$elem\n");
+ $i++;
+ }
+ $ost->print($RECORD_END);
+}
+
+sub Array::Size($)
+{
+ return scalar(@{$ARG[0]});
+}
+
+sub Array::Subset($$)
+{
+ my($aIn, $wantedIndices) = @ARG;
+
+ my $aOut = [];
+ for my $i (@$wantedIndices)
+ {
+ push @$aOut, $aIn->[$i];
+ }
+ return $aOut;
+}
+
+sub Array::SuccessiveDifference($)
+{
+ my $array = shift;
+ my $diff = [];
+ for (my $i=0; $i < @$array-1 ; $i++)
+ {
+ $diff->[$i] = $array->[$i+1] - $array->[$i];
+ }
+ return $diff;
+}
+
+
+sub Array::Sum($)
+{
+ my $array = shift;
+
+ my $sum = 0;
+ for my $elem (@$array)
+ {
+ $sum += $elem;
+ }
+ return $sum;
+}
+
+sub Array::SwapElements($$$)
+{
+ my($array, $i, $j) = @ARG;
+ my $tmp = $array->[$i];
+ $array->[$i] = $array->[$j];
+ $array->[$j] = $tmp;
+ return $array;
+}
+
+sub Array::Transform($$$)
+{
+ my($aIn, $aOut, $op) = @ARG;
+
+ for (my $i=0; $i<@$aIn; $i++)
+ {
+ $aOut->[$i] = &$op($aIn->[$i]);
+ }
+ return $aOut;
+}
+
+sub ToHash
+{
+ 2 == @ARG || 3 == @ARG or WrongNumArgsError();
+
+ my($array, $keyFunc) = @ARG;
+ my $valueFunc = (3 == @ARG) ? $ARG[2] : sub{$ARG[0]};
+
+ my $hash = {};
+ for my $elem (@$array)
+ {
+ $hash->{ &$keyFunc($elem) } = &$valueFunc($elem);
+ }
+ return $hash;
+}
+
+sub WhichMax($)
+{
+ my $array = shift;
+
+ my $maxValue = $array->[0];
+ my $maxIndex = 0;
+ for (my $i=1; $i<@$array; $i++)
+ {
+ if ($maxValue < $array->[$i])
+ {
+ $maxValue = $array->[$i];
+ $maxIndex = $i;
+ }
+ }
+ return $maxIndex;
+}
+
+sub Array::UniqueOrdered($)
+{
+ my $array = shift;
+ unless (@$array)
+ {
+ Fatal("array must be of non-zero size");
+ }
+
+ my $prev = $array->[0];
+ my $uniqArray = [$prev];
+ for (my $i=1; $i<@$array; ++$i)
+ {
+ if ($array->[$i] ne $prev)
+ {
+ push @$uniqArray, ($prev = $array->[$i]);
+ }
+ }
+ return $uniqArray;
+}
+
+
+sub Unique
+{
+ if (1==@ARG){return Unique_1arg(@ARG)}
+ if (2==@ARG){return Unique_2arg(@ARG)}
+ WrongNumArgsError();
+}
+
+sub Unique_1arg
+{
+ my $array = shift;
+
+ my %seen;
+
+ @seen{@$array} = ();
+
+ return [keys %seen];
+}
+
+sub Unique_2arg
+{
+ my($array, $cmpFunc) = @ARG;
+
+ my $hUniq = Array::ToHash($array, $cmpFunc);
+ return [values %$hUniq];
+}
+
+true;
diff --git a/bagpipe/lib/Array2d.pm b/bagpipe/lib/Array2d.pm
new file mode 100644
index 0000000..147216f
--- /dev/null
+++ b/bagpipe/lib/Array2d.pm
@@ -0,0 +1,214 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package Array2d;
+use strict;
+use English;
+use StdDefs;
+use Assert;
+use Array;
+use NamedArgs;
+
+my $RECORD_START = "#begin#Array2d\n";
+my $RECORD_END = "#end#Array2d\n";
+
+sub Flatten($)
+{
+ my $array2d = shift;
+ my $flatArray = [];
+
+ for my $array (@$array2d)
+ {
+ push @$flatArray, @$array;
+ }
+ return $flatArray;
+}
+
+sub Load($)
+{
+ my $ist = shift;
+
+ # skip to record start
+ return undef if $ist->eof;
+ while ($ist->getline ne $RECORD_START)
+ {
+ return undef if $ist->eof;
+ }
+
+ my $array = [];
+ while ((my $line = $ist->getline) ne $RECORD_END)
+ {
+ if ($ist->eof)
+ {
+ Fatal("Premature EOF encountered when reading object from file");
+ }
+ $line =~ m/^(\S+)\s+(\S+)\s+(\S.*)$/;
+ $array->[$1][$2] = $3;
+ }
+ return $array;
+}
+
+sub Array2d::Matrix
+# nrow is the 1st dimension
+# ncol is the 2nd dimension
+{
+ my($data, $hArgs) = NamedArgs::Parse(\@ARG, 1,
+ {'optional'=> {
+ 'dim1' => undef,
+ 'dim2' => undef,
+ 'bydim1' => false}
+ });
+
+ my $dim1 = $hArgs->{'dim1'};
+ my $dim2 = $hArgs->{'dim2'};
+
+ if (VarType::IsArrayRef($data))
+ {
+ if (not defined $dim1)
+ {
+ if (not defined $dim2)
+ {
+ $dim1 = @$data;
+ $dim2 = 1;
+ }
+ else
+ {
+ $dim1 = @$data/$dim2;
+ }
+ }
+ elsif (not defined $dim2)
+ {
+ $dim2 = @$data/$dim1
+ }
+
+ my $array2d = [];
+ my $k = 0;
+ if ($hArgs->{"bydim1"})
+ {
+ for (my $j = 0; $j<$dim2; $j++)
+ {
+ for (my $i = 0; $i<$dim1; $i++)
+ {
+ $array2d->[$i][$j] = $data->[$k++];
+ return $array2d if $k==@$data;
+ }
+ }
+ }
+ else
+ {
+ for (my $i = 0; $i<$dim1; $i++)
+ {
+ for (my $j = 0; $j<$dim2; $j++)
+ {
+ $array2d->[$i][$j] = $data->[$k++];
+ return $array2d if $k==@$data;
+ }
+ }
+
+ }
+ }
+ Fatal("not yet implemented\n");
+}
+
+sub New
+{
+ my $dim1Size = $ARG[0] || WrongNumArgsError();
+ my $dim2Size = (2 <= @ARG) ? $ARG[1] : 0;
+ my $value = (3 <= @ARG) ? $ARG[2] : 0;
+
+ my $array2d = [];
+ $#{$array2d} = $dim1Size -1;
+
+ for (my $i=0; $i<$dim1Size; $i++)
+ {
+ if ($dim2Size)
+ {
+ $array2d->[$i] = Array::New($dim2Size, $value);
+ }
+ else
+ {
+ $array2d->[$i] = [];
+ }
+ }
+ return $array2d;
+}
+
+sub NumElements($)
+{
+ my $array2d = shift;
+
+ my $sum = 0;
+ for my $array (@$array2d)
+ {
+ $sum += scalar(@$array);
+ }
+ return $sum;
+}
+
+sub Reset($$$$)
+{
+ my($array2d, $value) = @ARG;
+
+ for my $array (@$array2d)
+ {
+ for my $elem (@$array)
+ {
+ $elem = $value;
+ }
+ }
+}
+
+sub Save($$)
+{
+ my($ost, $array2d) = @ARG;
+
+ $ost->print($RECORD_START);
+ my $i = 0;
+ for my $array (@$array2d)
+ {
+ my $j = 0;
+ for my $elem (@$array)
+ {
+ $ost->print("$i\t$j\t$elem\n");
+ $j++;
+ }
+ $i++;
+ }
+ $ost->print($RECORD_END);
+}
+
+sub SwapAxes ($)
+{
+ return Transpose($ARG[0]);
+}
+
+sub Transpose($)
+{
+ my $array2d = shift;
+
+ my $swapped = [];
+ for (my $i=0; $i<@$array2d; $i++)
+ {
+ for (my $j=0; $j<@{$array2d->[$i]}; $j++)
+ {
+ $swapped->[$j][$i] = $array2d->[$i][$j];
+ }
+ }
+ return $swapped;
+}
+
+
+true;
diff --git a/bagpipe/lib/Assert.pm b/bagpipe/lib/Assert.pm
new file mode 100644
index 0000000..54e322f
--- /dev/null
+++ b/bagpipe/lib/Assert.pm
@@ -0,0 +1,139 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package Assert;
+use base qw(Exporter);
+use Exporter;
+use English;
+use Carp qw( confess cluck );
+use StdDefs;
+use strict;
+use vars qw(@EXPORT);
+
+ at EXPORT = qw( Assert Prefer WrongNumArgsError Fatal Warn Moan);
+
+sub Assert($$)
+{
+ my ($boolean, $message) = @ARG;
+
+ unless ($boolean)
+ {
+ confess "Null Assertion: $message\n";
+ }
+}
+
+
+sub Assert::CheckArgHash
+# usage:
+# CheckArgHash(\%hArg,
+# 'optional' => [qw(times length each)],
+# );
+#
+# CheckArgHash(\%hArg,
+# 'required' => [qw(
+# -parent
+# )],
+# 'optional' => [qw(
+# -initialfile
+# -dialogtitle
+# )],
+{
+ my($hSubArg, %hArg) = @ARG;
+
+ if (exists $hArg{'required'})
+ {
+ for my $arg (@{$hArg{'required'}})
+ {
+ unless (exists $hSubArg->{$arg})
+ {
+ Fatal("Missing required argument \'$arg\' to method");
+ }
+ }
+ }
+ if (exists $hArg{'optional'})
+ {
+ my @allowedList = @{$hArg{'optional'}};
+ push @allowedList, @{$hArg{'required'}}
+ if exists $hArg{'required'};
+ my %allowed = ();
+ @allowed{@allowedList} = ();
+ for my $arg (keys %$hSubArg)
+ {
+ unless (exists $allowed{$arg})
+ {
+ Fatal("Invalid argument \'$arg\' to method\n");
+ }
+ }
+ }
+ if (exists $hArg{'oneof'})
+ {
+ my $count = 0;
+ for my $arg (keys %$hSubArg)
+ {
+ $count++ if exists $hSubArg->{$arg};
+ }
+ if (1 != $count)
+ {
+ Fatal("Method must take exactly one of "
+ .join(',', $hArg{'oneof'})
+ );
+ }
+ }
+}
+
+sub Prefer
+{
+ my ($boolean, $message) = @ARG;
+
+ if( $boolean == false )
+ {
+ cluck "Warning: $message\n";
+ }
+}
+
+sub Fatal
+{
+ my $message = $ARG[0] || null;
+
+ confess "Exception: $message\n";
+}
+
+sub Warn
+{
+ my $message = $ARG[0] || null;
+
+ cluck "Warning: $message\n";
+}
+
+sub Moan
+{
+ my $message = $ARG[0] || null;
+
+ print STDERR "$message\n";
+}
+
+sub DeprecatedMethodWarning
+{
+ my $message = $ARG[0] || null;
+ cluck "Deprecated method: $message\n";
+}
+
+sub WrongNumArgsError ()
+{
+ confess "Error: Wrong number of arguments to method\n";
+}
+
+true;
diff --git a/bagpipe/lib/CCloneable.pm b/bagpipe/lib/CCloneable.pm
new file mode 100644
index 0000000..98f3ec5
--- /dev/null
+++ b/bagpipe/lib/CCloneable.pm
@@ -0,0 +1,28 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+package CCloneable;
+use English;
+use Hash;
+
+sub clone($)
+{
+ my $this = shift;
+
+ return bless Hash::Copy($this), ref($this);
+}
+1;
diff --git a/bagpipe/lib/CConfigFile.pm b/bagpipe/lib/CConfigFile.pm
new file mode 100644
index 0000000..d3da755
--- /dev/null
+++ b/bagpipe/lib/CConfigFile.pm
@@ -0,0 +1,172 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+package CConfigFile;
+use base qw(CCloneable);
+use English;
+use StdDefs;
+use Assert;
+use FileIoHelper qw(OpenFilesForReading OpenFilesForWriting);
+
+my $KEY_FILE_NAME = 'CConfigFile::KEY_FILE_NAME';
+
+sub new($$)
+{
+ my($class, $file) = @ARG;
+
+ my $table = {};
+
+ # file table from file
+ _ReadConfig(OpenFilesForReading($file), $table);
+
+ $table->{$KEY_FILE_NAME} = $file;
+
+ # return as object
+ return bless $table, $class;
+}
+
+sub add($$$)
+{
+ my($this, $key, $value) = @ARG;
+ if ($this->has($key))
+ {
+ Fatal("can only add() new keys");
+ }
+ $this->put($key, $value);
+}
+
+sub get($$)
+{
+ my($this, $key) = @ARG;
+ if ($this->has($key))
+ {
+ return $this->{$key};
+ }
+ else
+ {
+ Fatal("config key \'$key\' not found");
+ }
+}
+
+sub getKeys($)
+{
+ my $this = shift;
+ my @keys = ();
+ for my $key (keys %$this)
+ {
+ next if $key eq $KEY_FILE_NAME;
+ push @keys, $key;
+ }
+ return @keys;
+}
+
+sub getName($)
+{
+ my $this = shift;
+ return $this->{$KEY_FILE_NAME};
+}
+
+sub has($$)
+{
+ my($this, $key) = @ARG;
+ return exists $this->{$key};
+}
+
+sub put($$$)
+{
+ my($this, $key, $value) = @ARG;
+ my $oldValue = $this->{$key};
+ $this->{$key} = $value;
+ return $oldValue;
+}
+
+sub reload($)
+{
+ my $this = shift;
+ _ReadConfig(OpenFilesForReading($this->getName), $this);
+ return $this;
+}
+
+sub remove($$)
+{
+ my($this, $key) = @ARG;
+ if ($this->has($key))
+ {
+ return delete $this->{$key};
+ }
+ else
+ {
+ Fatal("cannot delete, key \'$key\' doesn't exist");
+ }
+}
+
+sub set($$$)
+{
+ my($this, $key, $value) = @ARG;
+ unless ($this->has($key))
+ {
+ Fatal("can only set() existing values");
+ }
+ return $this->put($key, $value);
+}
+
+sub setName($$)
+{
+ my($this, $file) = @ARG;
+ $this->{$KEY_FILE_NAME} = $file;
+ return $this;
+}
+
+sub save($)
+{
+ my $this = shift;
+ _SaveConfig($this, $this->getName);
+}
+
+sub _ReadConfig($$)
+{
+ my($ist, $table) = @ARG;
+ # read config file
+ while (my $line = $ist->getline)
+ {
+ if ($line =~ m/^(\S+)\s+(\S.*?)\s*$/)
+ {
+ $table->{$1} = $2;
+ }
+ }
+}
+
+sub _SaveConfig($$)
+{
+ my($this, $file) = @ARG;
+ my $ost = OpenFilesForWriting($file);
+ $this->write($ost);
+}
+
+sub write($$)
+{
+ my($this, $ost) = @ARG;
+ for my $key (sort keys %$this)
+ {
+ next if $key eq $KEY_FILE_NAME;
+ $ost->print("$key\t".$this->{$key}."\n");
+ }
+}
+
+true;
+
+
diff --git a/bagpipe/lib/CDataFrame.pm b/bagpipe/lib/CDataFrame.pm
new file mode 100644
index 0000000..423541e
--- /dev/null
+++ b/bagpipe/lib/CDataFrame.pm
@@ -0,0 +1,920 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+package CDataFrame;
+use English;
+use VarType;
+use StdDefs;
+use Array;
+use Array2d;
+use Assert;
+use NamedArgs;
+use String;
+use VarType qw(AsArray);
+use Algorithm qw(ForEach);
+use Debug qw(DumpObject);
+
+my $COLUMN_ORDER_FIELD = '__COLUMN_ORDER';
+
+sub CDataFrame::new
+{
+ my($class, $hArgs) = NamedArgs::Parse(\@ARG, 1,
+ {'optional'=>
+ {
+ 'format' => 'normal',
+ 'header' => false,
+ 'hash' => null,
+ 'matrix' => null,
+ 'pad_string' => null,
+ 'stream' => null,
+ 'strip.white' => true,
+ 'hrmatrix' => null, # deprecated
+ #S/R read.table: 'fill' => true,
+ # 'blank.lines.skip'=> true,
+ # 'dec' => ".",
+
+ }
+ });
+ my $this = bless {
+ '__COLUMN_ORDER'=>{},
+ }, $class;
+
+
+ #---------------------------------------
+ # deprecated part
+ if ($hArgs->{'hrmatrix'})
+ {
+ $hArgs->{'matrix'} = $hArgs->{'hrmatrix'};
+ $hArgs->{'header'} = true;
+ Warn("deprecated argument \'hrmatrix\'");
+ }
+ #----------------------------------------
+
+ if ($hArgs->{'matrix'})
+ {
+ $this->_readMatrix($hArgs);
+ }
+ elsif ($hArgs->{'stream'})
+ {
+ $this->_readStream($hArgs);
+ }
+ elsif ($hArgs->{'hash'})
+ {
+ $this->_readHash($hArgs);
+ }
+ return $this;
+}
+
+sub CDataFrame::_readStream($$)
+{
+ my($this, $hArgs) = @ARG;
+
+ if ($hArgs->{'format'} =~ m/normal/i)
+ {
+ $hArgs->{'sep'} = '\s';
+ $hArgs->{'quote'} = "\"";
+ $this->_parseStream($hArgs);
+ }
+ elsif ($hArgs->{'format'} =~ m/delim/i)
+ {
+ $hArgs->{'sep'} = "\t";
+ $hArgs->{'quote'} = "\"";
+ $this->_parseStream($hArgs);
+ }
+ elsif ($hArgs->{'format'} =~ m/splus/i)
+ {
+ $this->_readStreamSplus($hArgs);
+ }
+ elsif ($hArgs->{'format'} =~ m/csv/i)
+ {
+ $hArgs->{'sep'} = ",";
+ $hArgs->{'quote'} = "\"";
+ $this->_parseStream($hArgs);
+ }
+ else
+ {
+ Fatal("Unsupported format \'".$hArgs->{'format'}."\'");
+ }
+}
+
+sub CDataFrame::_parseStream($$)
+{
+ my($this, $hArgs) = @ARG;
+
+ my $ist = $hArgs->{'stream'};
+ my $quote = quotemeta $hArgs->{'quote'};
+ my $sep = $hArgs->{'sep'};
+
+ my $aRows = [];
+ while (my $line = $ist->getline)
+ {
+ next unless $line =~ m/\S/;
+ $line = String::Trim($line);
+
+ # avoid dealing with quotes if possible
+ if ($line !~ m/$quote/)
+ {
+ #warn("no quote\n");
+ push @$aRows, [split m/$sep/o, $line];
+ next;
+ }
+
+ # deal with quotes
+ my $row = [];
+ unless (Number::IsEven(String::CountChar($line, $quote)))
+ {
+ Fatal("Odd number of quote characters [$quote] in [$line]");
+ }
+
+ while ($line =~ m/\S/)
+ {
+ if ($line =~ s/^([^$quote]+)$sep*//)
+ {
+ push @$row, split m/$sep/o, $1;
+ }
+ if ($line =~ s/^(${quote}[^$quote]+${quote})$sep?//)
+ {
+ push @$row, $1;
+ }
+ }
+ push @$aRows, $row;
+
+ }
+
+ # strip white if desired
+ if ($hArgs->{'strip.white'})
+ {
+ for my $row (@$aRows)
+ {
+ $row = String::TrimArray($row);
+ }
+ }
+
+ $hArgs->{'matrix'} = $aRows;
+
+ $this->_readMatrix($hArgs);
+}
+
+sub CDataFrame::_readStreamSplus($$)
+{
+ my($this, $hArgs) = @ARG;
+
+ my $ist = $hArgs->{'stream'};
+ my $aRows = [];
+
+ if ($hArgs->{'header'})
+ {
+ my $line = null;
+ while (($line = $ist->getline) !~ m/\S/){}
+ push @$aRows, _ParseSplusLine($line);
+ }
+ while (my $line = $ist->getline)
+ {
+ next unless $line =~ m/\S/;
+ my $lineArray = _ParseSplusLine($line);
+ shift @$lineArray;
+ push @$aRows, $lineArray;
+ }
+ $hArgs->{'matrix'} = $aRows;
+
+ $this->_readMatrix($hArgs);
+}
+
+sub CDataFrame::_readMatrix($$)
+{
+ my($this, $hArgs) = @ARG;
+
+ my $aRows = $hArgs->{'matrix'};
+
+ # get max number of columns
+ my $maxRowLength = 0;
+ my $maxRowLengthIndex = 0;
+ for (my $i=0; $i<@$aRows; ++$i)
+ {
+ my $row = $aRows->[$i];
+ if ($maxRowLength < @$row)
+ {
+ $maxRowLength = @$row;
+ $maxRowLengthIndex = $i;
+ }
+ }
+
+ # get or make col names
+ my $aColNames = null;
+ if ($hArgs->{'header'})
+ {
+ $aColNames = shift @$aRows;
+ if (@$aColNames < $maxRowLength)
+ {
+ Fatal("Must have at least as many column headers as data columns:"
+ ." headers = ".@$aColNames.", max row length = "
+ .$maxRowLength. " at row ".($maxRowLengthIndex+1));
+ }
+ }
+ else
+ {
+ $aColNames = Array::Paste(["X"], [1..scalar(@{$aRows->[0]})],
+ "sep"=>null);
+ }
+ $this->setColOrder($aColNames);
+
+ # pad data rows if necessary
+ for my $row (@$aRows)
+ {
+ if ($maxRowLength > @$row)
+ {
+ my $pad = Array::New($maxRowLength-@$row,
+ $hArgs->{'pad_string'});
+ push @$row, @$pad;
+ }
+ }
+
+ # make df
+ for my $colName (@$aColNames)
+ {
+ $this->{$colName} = [];
+ }
+ for my $row (@$aRows)
+ {
+ for (my $i=0; $i<@$row; $i++)
+ {
+ push @{$this->{$aColNames->[$i]}}, $row->[$i];
+ }
+ }
+}
+
+sub CDataFrame::_readHash($$)
+{
+ my($this, $hArgs) = @ARG;
+
+ my $hash = $hArgs->{'hash'};
+
+ while (my($key, $value) = each %$hash)
+ {
+ $this->addCol($key, $value);
+ }
+}
+
+sub CDataFrame::addCol($$$)
+{
+ my($this, $name, $value, $hArgs) = NamedArgs::Parse(\@ARG, 3,
+ {
+ 'optional' =>
+ {
+ 'expand' => false,
+ }
+ }
+ );
+
+ unless ($this->hasCol($name))
+ {
+ my $rank = undef;
+ if (scalar keys %{ $this->{$COLUMN_ORDER_FIELD} })
+ {
+ $rank = 1 + Array::Max([values %{$this->{$COLUMN_ORDER_FIELD}}]);
+ }
+ else
+ {
+ $rank = 1;
+ }
+ $this->{$COLUMN_ORDER_FIELD}{$name} = $rank;
+ }
+
+ # expand if necessary
+ if ($hArgs->{'expand'})
+ {
+ return $this->addExpandingCol($name, $value);
+ }
+
+ my $array = VarType::AsArray($value);
+ if ($this->getNumCols and $this->getNumRows != @$array)
+ {
+ Fatal("cannot add array of size ".scalar(@$array)
+ ." to CDataFrame of size ".$this->getNumRows);
+ }
+ $this->{$name} = $array;
+ return $this;
+}
+
+sub CDataFrame::addExpandingCol($$$)
+{
+ my($this, $name, $array) = @ARG;
+ $array = VarType::AsArray($array);
+
+ if ($this->getNumRows < @$array)
+ {
+ Fatal("array is too large to expand");
+ }
+ my $aExpanded = Array::Rep(
+ $array,
+ 'length' => $this->getNumRows,
+ 'each' => 1
+ );
+ $this->{$name} = $aExpanded;
+ return $this;
+}
+
+sub CDataFrame::addRow($$)
+{
+ my($this, $hash) = @ARG;
+ unless ($this->getNumCols == scalar keys %$hash)
+ {
+ Fatal("added rows must have same column names as dataframe:\n"
+ ."is => ".join("|", sort keys %$hash)."\n"
+ ."ought => ".join("|", sort keys %$this)."\n");
+ }
+ while (my($key, $value) = each %$hash)
+ {
+ unless (exists $this->{$key})
+ {
+ Fatal("could not find column named \'$key\' in dataframe");
+ }
+ push @{$this->{$key}}, $value;
+ }
+ return $this;
+}
+
+sub CDataFrame::cbind($$)
+{
+ my($this, $df, $hArgs) = NamedArgs::Parse(\@ARG, 2,
+ {
+ 'optional' =>
+ {
+ 'expand' => false,
+ }
+ });
+ for my $col (@{$df->getColOrder})
+ {
+ $this->addCol($col, $df->getCol($col), 'expand'=>$hArgs->{'expand'});
+ }
+
+
+}
+
+sub CDataFrame::clone($)
+{
+ my $this = shift;
+ my $copy = {};
+ while (my($key, $value) = each %$this)
+ {
+ if ($key eq $COLUMN_ORDER_FIELD)
+ {
+ $copy->{$key} = Hash::Copy($value);
+ }
+ else
+ {
+ $copy->{$key} = [@$value];
+ }
+ }
+ return bless $copy, "CDataFrame";
+}
+
+sub CDataFrame::crop
+{
+ my($this, %hArg) = @ARG;
+
+ # rm null column
+ delete $this->{''} if exists $this->{''};
+
+ # search for trailing blanks
+ my $maxLastIndex = 0;
+ for my $name (keys %$this)
+ {
+ my $array = $this->{$name};
+ my $lastIndex = Array::RevFindIf($array,
+ sub
+ {
+ my $str = shift;
+ return false if not defined $str;
+ if ($hArg{'rmnulls'} and 0==length($str))
+ {
+ return false;
+ }
+ return true;
+ }
+ );
+ if (Array::NOT_FOUND eq $lastIndex)
+ {
+ $lastIndex = scalar(@$array) - 1;
+ }
+ if ($lastIndex > $maxLastIndex)
+ {
+ $maxLastIndex = $lastIndex;
+ }
+ }
+
+ # crop trailing blanks
+ while (my($name, $array) = each %$this)
+ {
+ Array::Crop($array, 0, $maxLastIndex);
+ }
+
+ return $this;
+}
+
+
+sub CDataFrame::filterRows
+{
+ my ($this, %hFilter) = @ARG;
+
+ my $nrow = $this->getNumRows;
+ my $mask = Array::Rep(true, 'length'=>$nrow);
+
+ while (my($name, $value) = each %hFilter)
+ {
+ my $array = $this->getCol($name);
+ for (my $i=0; $i<$nrow; ++$i)
+ {
+ $mask->[$i] = false if $array->[$i] ne $value;
+ }
+ }
+ return $this->mask($mask);
+}
+
+sub CDataFrame::getCell($$$)
+{
+ my($this, $rowIndex, $colName) = @ARG;
+ unless ($this->hasCell($rowIndex, $colName))
+ {
+ Fatal("No such cell ($rowIndex,\'$colName\')");
+ }
+ return $this->{$colName}[$rowIndex];
+}
+
+sub CDataFrame::getCol($$)
+{
+ my($this, $name) = @ARG;
+ unless (exists $this->{$name})
+ {
+ Fatal("Column \'$name\' does not exist");
+ }
+ return $this->{$name};
+}
+
+sub CDataFrame::getColNames($)
+{
+ my $this = shift;
+ my $aNames = [];
+ for my $name (sort keys %$this)
+ {
+ next if $name eq $COLUMN_ORDER_FIELD;
+ push @$aNames, $name;
+ }
+ return $aNames;
+}
+
+sub CDataFrame::getColOrder($)
+{
+ my $this = shift;
+ my $hKnown2Rank = $this->{$COLUMN_ORDER_FIELD};
+
+ my $aKnown = [];
+ my $aUnknown = [];
+ for my $name (@{ $this->getColNames })
+ {
+ if (exists $hKnown2Rank->{$name})
+ {
+ push @$aKnown, $name;
+ }
+ else
+ {
+ push @$aUnknown, $name;
+ }
+ }
+
+ my $aNames = [
+ sort
+ { $hKnown2Rank->{$a} <=> $hKnown2Rank->{$b} }
+ @$aKnown
+ ];
+
+ push @$aNames, sort @$aUnknown;
+ return $aNames;
+}
+
+sub CDataFrame::getNumCols($)
+{
+ my $this = shift;
+ return (scalar keys %$this) -1;
+}
+
+sub CDataFrame::getNumRows($)
+{
+ my $this = shift;
+ my $aCols = $this->getColNames;
+ if (@$aCols)
+ {
+ return scalar(@{ $this->{$aCols->[0]} });
+ }
+ return 0;
+}
+
+sub CDataFrame::getRow
+{
+ my $this = shift @ARG;
+ my $i = shift @ARG;
+ my $cols = (@ARG) ? shift @ARG : $this->getColOrder;
+
+ unless ($this->hasRow($i))
+ {
+ Fatal("Cannot get non-existant row $i\n");
+ }
+
+ my $row = [];
+ for my $col (@$cols)
+ {
+ push @$row, $this->{$col}[$i];
+ }
+ return $row;
+}
+
+sub CDataFrame::getRows($$)
+{
+ my($this, $aIndices) = @ARG;
+
+ my $subDf = new CDataFrame;
+ for my $field (@{$this->getColNames})
+ {
+ $subDf->addCol($field,
+ Array::Subset($this->getCol($field), $aIndices));
+ }
+ $subDf->setColOrder($this->getColOrder);
+
+ return $subDf;
+}
+
+sub CDataFrame::hasCell($$$)
+{
+ my($this, $rowIndex, $colName) = @ARG;
+ return ($this->hasRow($rowIndex) and $this->hasCol($colName))
+ ? true : false;
+}
+
+sub CDataFrame::hasCol($$)
+{
+ my($this, $name) = @ARG;
+ return exists $this->{$name};
+}
+
+sub CDataFrame::hasRow($$)
+{
+ my($this, $rowIndex) = @ARG;
+ return 0 <= $rowIndex and $rowIndex < $this->getNumRows;
+}
+
+sub CDataFrame::isEmpty($)
+{
+ my $this = shift;
+ return 0==$this->getNumCols;
+}
+
+
+sub CDataFrame::mask($$)
+{
+ my($this, $mask) = @ARG;
+
+ unless (scalar(@$mask) == $this->getNumRows)
+ {
+ Fatal("mask must have as many elements as data frame has rows");
+ }
+
+ my $maskedDf = new CDataFrame;
+ for my $field (@{$this->getColNames})
+ {
+ $maskedDf->addCol($field,
+ Array::Mask($this->getCol($field), $mask));
+ }
+ $maskedDf->setColOrder($this->getColOrder);
+ return $maskedDf;
+}
+
+sub CDataFrame::rbind($$)
+{
+ my($a, $b, $hArgs) = NamedArgs::Parse(\@ARG, 2,
+ {
+ 'optional' =>
+ {
+ 'expand' => false,
+ }
+ }
+ );
+
+ # deal with trivial cases
+ return $a->clone if $b->isEmpty;
+ return $b->clone if $a->isEmpty;
+
+ # see if column names match...
+ unless (Array::Equals($a->getColNames, $b->getColNames))
+ {
+ return $a->_rbindExpand($b) if $hArgs->{'expand'};
+ Fatal("CDataFrames must have the same column names");
+ }
+
+ # bind tables together
+ my $c = new CDataFrame;
+ for my $name (@{$a->getColNames})
+ {
+ $c->addCol($name, [@{$a->getCol($name)}, @{$b->getCol($name)}]);
+ }
+ return $c;
+}
+
+sub CDataFrame::_rbindExpand($$)
+{
+ my($a, $b) = @ARG;
+
+ my($notB, $notA, $comm) = Array::Diff($a->getColNames, $b->getColNames);
+
+ ForEach($notB, sub{$b->addExpandingCol($ARG[0], null)});
+ ForEach($notA, sub{$a->addExpandingCol($ARG[0], null)});
+
+ return $a->rbind($b);
+}
+
+sub CDataFrame::rqbind($$)
+{
+ my($this, $df) = @ARG;
+
+ return $this if $df->isEmpty;
+
+ unless ($this->isEmpty
+ or
+ Array::Equals($this->getColNames, $df->getColNames))
+ {
+ Fatal("CDataFrames must have the same column names:\n[",
+ join(@{$this->getColOrder}, ","), "]\n[",
+ join(@{$df->getColOrder}, ","), "]");
+ }
+
+ for my $name (@{$df->getColNames})
+ {
+ push @{$this->{$name}}, @{$df->{$name}};
+ }
+ return $this;
+}
+
+sub CDataFrame::removeCol($$)
+{
+ my($this, $colName) = @ARG;
+
+ my $array = $this->getCol($colName);
+
+ delete $this->{$colName};
+ delete $this->{$COLUMN_ORDER_FIELD}{$colName};
+
+ return $array;
+}
+
+sub CDataFrame::removeCols($$)
+{
+ my($this, $aColNames) = @ARG;
+ for my $colName (@$aColNames)
+ {
+ $this->removeCol($colName);
+ }
+}
+
+sub CDataFrame::renameCol($$$)
+{
+ my($this, $oldName, $newName) = @ARG;
+
+ unless ($this->hasCol($oldName))
+ {
+ Fatal("Cannot rename non-existant column \"$oldName\"");
+ }
+
+ $this->{$newName} = delete $this->{$oldName};
+ $this->{$COLUMN_ORDER_FIELD}{$newName} =
+ delete $this->{$COLUMN_ORDER_FIELD}{$oldName};
+ return $this;
+}
+
+sub CDataFrame::renameCols($$)
+{
+ my($this, $hOld2New) = @ARG;
+
+ while (my($old, $new)=each %$this)
+ {
+ $this->renameCol($old, $new);
+ }
+ return $this;
+}
+
+sub CDataFrame::setCol($$$)
+{
+ my($this, $name, $value) = @ARG;
+ return $this->addCol($name, $value);
+}
+
+sub CDataFrame::setCell($$$$)
+{
+ my($this, $rowIndex, $colName, $value) = @ARG;
+ unless ($this->hasCell($rowIndex, $colName))
+ {
+ Fatal("No such cell ($rowIndex,\'$colName\')");
+ }
+ $this->{$colName}[$rowIndex] = $value;
+}
+
+sub CDataFrame::setColOrder($$)
+{
+ my($this, $aColNames) = @ARG;
+ $this->{$COLUMN_ORDER_FIELD} = {};
+ my $count = 0;
+ for my $name (@$aColNames)
+ {
+ $this->{$COLUMN_ORDER_FIELD}{$name} = ++$count;
+ }
+}
+
+sub CDataFrame::transpose($$)
+{
+ my($this, $pivot) = @ARG;
+
+ unless ($this->hasCol($pivot))
+ {
+ Fatal("Pivot column \"$pivot\" does not exist");
+ }
+
+ my $df = new CDataFrame;
+ my $oldCols = Array::Exclude($this->getColOrder, $pivot);
+ $df->addCol($pivot, $oldCols);
+ for (my $i=0; $i<$this->getNumRows; ++$i)
+ {
+ my $newCol = $this->getCell($i, $pivot);
+ my $data = $this->getRow($i, $oldCols);
+ $df->addCol($newCol, $data);
+ }
+ return $df;
+}
+
+sub CDataFrame::write
+{
+ my($this, $hArgs) = NamedArgs::Parse(\@ARG, 1,
+ {
+ 'required' => [qw(stream)],
+ 'optional' =>
+ {
+ 'format' => 'normal',
+ 'header' => true,
+ 'cols' => null,
+ }
+ }
+ );
+
+ my $matrix = $this->asMatrix(
+ 'header'=>$hArgs->{'header'},
+ 'transpose'=>false,
+ 'cols' => $hArgs->{'cols'},
+ );
+
+ my $ost = $hArgs->{'stream'};
+ if ('normal' eq $hArgs->{'format'} or 'delim' eq $hArgs->{'format'})
+ {
+ for my $row (@$matrix)
+ {
+ $ost->print(join("\t", @$row), "\n");
+ }
+ }
+ elsif ($hArgs->{'format'} =~ m/csv/i)
+ {
+ for my $row (@$matrix)
+ {
+ $ost->print(join(",", @$row), "\n");
+ }
+ }
+ elsif ($hArgs->{'format'} =~ m/splus/i)
+ {
+ $ost->print(_MakeSplusLine($matrix->[0]), "\n");
+ for (my $i=1; $i<@$matrix; $i++)
+ {
+ unshift @{$matrix->[$i]}, $i;
+ $ost->print(_MakeSplusLine($matrix->[$i]), "\n")
+ }
+ }
+ else
+ {
+ Fatal("Unsupported write format \'".$hArgs->{'format'}."\'");
+ }
+}
+
+sub CDataFrame::asMatrix
+{
+ my($this, $hArgs) = NamedArgs::Parse(\@ARG, 1,
+ {
+ 'optional' =>
+ {
+ 'header' => true,
+ 'transpose' => false,
+ 'cols' => null,
+ }
+ });
+ my $matrixT = [];
+ my $aColNames = $hArgs->{'cols'}
+ ? $hArgs->{'cols'}
+ : $this->getColOrder;
+
+ for my $colName (@$aColNames)
+ {
+ my $rowT = Array::Copy($this->getCol($colName));
+ if ($hArgs->{'header'})
+ {
+ unshift @$rowT, $colName;
+ }
+ push @$matrixT, $rowT;
+ }
+
+ if ($hArgs->{'transpose'})
+ {
+ return $matrixT;
+ }
+ else
+ {
+ return Array2d::Transpose($matrixT);
+ }
+}
+
+sub CDataFrame::subset
+{
+ my($this, $hArgs) = NamedArgs::Parse(\@ARG, 1,
+ {
+ 'optional' =>
+ {
+ 'rows' => null,
+ 'cols' => null,
+ }
+ }
+ );
+ my $aWantedRows = (null ne $hArgs->{'rows'})
+ ? AsArray($hArgs->{'rows'})
+ : [0..($this->getNumRows-1)];
+ my $aWantedCols = (null ne $hArgs->{'cols'})
+ ? AsArray($hArgs->{'cols'})
+ : $this->getColNames;
+
+ my $df = new CDataFrame;
+ for my $wantedCol (@$aWantedCols)
+ {
+ my $data = Array::Subset($this->getCol($wantedCol),$aWantedRows);
+ $df->addCol($wantedCol, $data);
+ }
+ return $df;
+}
+
+sub _ParseSplusLine($)
+{
+ my $line = String::Trim(shift);
+
+ my $array = [];
+ while (true)
+ {
+ if ($line =~ s/^\"\"\s*//)
+ {
+ push @$array, null;
+ }
+ elsif ($line =~ s/^(\"?)([^\"]+)(\1?)\s*//)
+ {
+ my $leader = $1;
+ my $middle = $2;
+ my $trailer = $3;
+
+ if ($leader eq '"')
+ {
+ push @$array, $middle;
+ }
+ else
+ # for items not enclosed by ""s,
+ # split on white space
+ {
+ push @$array, split(m/\s+/, $middle);
+ }
+ }
+ last unless $line =~ m/\S/;
+ }
+ return $array;
+}
+
+sub _MakeSplusLine($)
+{
+ my $array = shift;
+ $array = Array::Apply($array, sub{"\"".$ARG[0]."\""});
+ return join(" ", @$array);
+}
+
+
+
+1;
diff --git a/bagpipe/lib/CRing.pm b/bagpipe/lib/CRing.pm
new file mode 100644
index 0000000..def23de
--- /dev/null
+++ b/bagpipe/lib/CRing.pm
@@ -0,0 +1,130 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+#--------------------------------------------------------------------------
+package CRingNode;
+use English;
+use StdDefs;
+
+sub new($$)
+{
+ my($class, $value) = @ARG;
+ my $this = {
+ 'leftLink' => null,
+ 'rightLink' => null,
+ 'value' => $value,
+ };
+ return bless $this, $class;
+}
+
+sub getLeft($) {$ARG[0]->{'leftLink'}}
+sub getRight($) {$ARG[0]->{'rightLink'}}
+sub getValue($) {$ARG[0]->{'value'}}
+
+sub setLeft($$) {$ARG[0]->{'leftLink'} = $ARG[1]}
+sub setRight($$) {$ARG[0]->{'rightLink'} = $ARG[1]}
+
+#--------------------------------------------------------------------------
+package CRing;
+use English;
+use StdDefs;
+use Array;
+
+sub new
+{
+ my($class) = shift;
+ my $array = @ARG ? shift
+ : [];
+ my $this = {
+ 'length' => scalar(@$array),
+ 'head' => _MakeRing($array),
+ };
+ return bless $this, $class;
+}
+
+sub get($$)
+{
+ my($this, $i) = @ARG;
+ return $this->_getNode($i)->getValue;
+}
+
+sub getSize($) {$ARG[0]{'length'}}
+
+sub nextElement($)
+{
+ my $this = shift;
+ my $currNode = $this->_getHead;
+ $this->rotate(1);
+ return $currNode->getValue;
+}
+
+sub rotate($$)
+{
+ my($this, $nMoves) = @ARG;
+
+ $this->_setHead($this->_getNode($nMoves));
+
+ return $this;
+}
+
+sub _getHead($) {$ARG[0]->{'head'}}
+
+sub _getNode($$)
+{
+ my($this, $nPlaces) = @ARG;
+
+ my $node = $this->_getHead;
+ if ($nPlaces >= 0)
+ {
+ for (my $i=1; $i<=$nPlaces; $i++)
+ {
+ $node = $node->getRight;
+ }
+ }
+ else # $nPlaces is negative
+ {
+ for (my $i=1; $i<=-$nPlaces; $i++)
+ {
+ $node = $node->getLeft;
+ }
+ }
+ return $node;
+}
+
+sub _setHead($) {$ARG[0]->{'head'} = $ARG[1]}
+
+sub _MakeRing($)
+{
+ my($array) = @ARG;
+
+ my $aNodes = Array::Transform($array, [], sub{new CRingNode($ARG[0])});
+ my $nNodes = scalar(@$aNodes);
+
+ for (my $i=0; $i<$nNodes; $i++)
+ {
+ my $leftIndex = (0==$i)
+ ? $nNodes-1
+ : $i-1;
+ my $rightIndex = ($nNodes-1==$i)
+ ? 0
+ : $i+1;
+ $aNodes->[$i]->setLeft($aNodes->[$leftIndex]);
+ $aNodes->[$i]->setRight($aNodes->[$rightIndex]);
+ }
+ return $aNodes->[0];
+}
+1;
diff --git a/bagpipe/lib/Debug.pm b/bagpipe/lib/Debug.pm
new file mode 100644
index 0000000..da2deaf
--- /dev/null
+++ b/bagpipe/lib/Debug.pm
@@ -0,0 +1,68 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+package Debug;
+use English;
+use Exporter;
+use Dumpvalue;
+use FileIoHelper;
+
+use base ('Exporter');
+use vars qw(@EXPORT_OK @EXPORT);
+ at EXPORT_OK = qw( DumpObject DumpHere );
+ at EXPORT = qw(DumpObject DumpHere DumpSep);
+
+my $ost = GetStdOut();
+
+sub DumpHere()
+{
+ my($pkg, $file, $line) = caller();
+ print $ost "Reached line $line [$file]\n";
+}
+
+
+sub DumpObject
+{
+ my $var = shift;
+ my($pkg, $file, $line) = caller();
+ print $ost "DumpObject() dumping from $file:$line\n";
+ DumpValue($var);
+ return 0;
+}
+
+sub DumpSep()
+{
+ print $ost "-" x 80, "\n";
+}
+
+sub DumpValue
+{
+ my $var = shift;
+ my $oldFh = select $ost;
+ (new Dumpvalue)->dumpValue($var);
+ select $oldFh;
+}
+
+sub SetDumpStream($)
+{
+ $ost = shift;
+ my $oldFh = select $ost;
+ $OUTPUT_AUTOFLUSH = 1;
+ select $oldFh;
+}
+
+1;
diff --git a/bagpipe/lib/FileIoHelper.pm b/bagpipe/lib/FileIoHelper.pm
new file mode 100644
index 0000000..0a04a11
--- /dev/null
+++ b/bagpipe/lib/FileIoHelper.pm
@@ -0,0 +1,461 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+use strict;
+#==============================================================================
+# MODULE: FileIoHelper
+#=====================
+# DESCRIPTION:
+#
+
+#_ Module declaration
+
+package FileIoHelper;
+
+#_ Include libraries
+
+use English;
+use StdDefs;
+use FileHandle;
+use DirHandle;
+use English;
+use Carp;
+use Assert;
+
+#_ Export
+
+use Exporter;
+use base ('Exporter');
+use vars qw(@EXPORT @EXPORT_OK);
+
+ at EXPORT = qw(
+ CloseFiles
+ GetStdError
+ GetStdIn
+ GetStdOut
+ OpenFilesForAppending
+ OpenFilesForReading
+ OpenFilesForWriting
+ );
+
+
+ at EXPORT_OK = qw(
+ CloseFiles
+ ColsFromFile
+ ColsFromStream
+ DirExists
+ FileExists
+ FileIsBinary
+ FileIsEmpty
+ FileEmptyOrAbsent
+ FileNotEmpty
+ GetDirContents
+ GetStdError
+ GetStdIn
+ GetStdOut
+ GetTempFileName
+ GetFileBasename
+ GetFilePath
+ GetFileSuffix
+ LinesFromFile
+ OpenFilesForAppending
+ OpenFilesForReading
+ OpenFilesForWriting
+ SplitPathName
+ println
+ );
+
+
+
+#_ Module variables
+
+my $tempFileNumber = 0;
+my $tempFilePrefix = 'tmp';
+my $tempFileSuffix = '.tmp';
+
+#_ Module init
+
+my $fhStdIn = bless \*STDIN, 'FileHandle';
+my $fhStdOut = bless \*STDOUT, 'FileHandle';
+my $fhStdErr = bless \*STDERR, 'FileHandle';
+
+#_ Public functions
+
+sub CloseFiles
+{
+ foreach my $fhToClose (@ARG)
+ {
+ $fhToClose->close()
+ or confess "Exception: Error closing file: $ERRNO\n";
+ }
+ return;
+}
+
+sub ColsFromFile($$)
+{
+ my($fileName, $aWantedCols) = @ARG;
+
+ return ColsFromStream(OpenFilesForReading($fileName), $aWantedCols);
+}
+
+sub ColsFromStream($$)
+{
+ my($ist, $aWantedCols) = @ARG;
+ my @aLists = ();
+
+ while(my $line = $ist->getline)
+ {
+ if($line =~ m/\S/)
+ {
+ my @aCols = split(m/\s+/, $line);
+ for(my $i=0; $i<@$aWantedCols; $i++)
+ {
+ push @{$aLists[$i]}, $aCols[$aWantedCols->[$i]];
+ }
+ }
+ }
+ return wantarray ? @aLists : $aLists[0];
+}
+
+sub DirExists($) {-d $ARG[0]}
+sub FileExists($) {-e $ARG[0]}
+sub FileIsBinary($) {-B $ARG[0]}
+sub FileIsEmpty($) {-z $ARG[0]}
+
+sub FileNotEmpty($)
+{
+ my($file) = @ARG;
+ return -e $file && !-z $file;
+}
+
+sub FileEmptyOrAbsent($)
+{
+ return !FileNotEmpty($ARG[0]);
+}
+
+sub GetStdIn() {$fhStdIn}
+sub GetStdOut() {$fhStdOut}
+sub GetStdError() {$fhStdErr}
+
+sub GetDirContents
+{
+ my($dirName, %hArgs) = @ARG;
+
+ Assert::CheckArgHash(\%hArgs,
+ 'optional'=>[qw(incpath nodots)]
+ );
+
+ my $aFiles = [];
+ my $aDirs = [];
+
+ my $dir = new DirHandle($dirName);
+ unless ($dir)
+ {
+ Fatal("Exception: cannot open directory $dirName: $ERRNO\n");
+ }
+
+ while(my $fileName = $dir->read)
+ {
+ if ($hArgs{'nodots'})
+ {
+ next if $fileName eq "." or $fileName eq "..";
+ }
+ my $pathName = "$dirName/$fileName";
+ my $name = $hArgs{'incpath'} ? $pathName
+ : $fileName;
+ if(-d $pathName)
+ {
+ push @$aDirs, $name;
+ }
+ else
+ {
+ push @$aFiles, $name;
+ }
+ }
+ $dir->close;
+ return {'files' => $aFiles, 'dirs' => $aDirs};
+}
+
+sub GetTempFileName()
+{
+ my $fileName = $tempFilePrefix."_".$PROCESS_ID."_".($tempFileNumber++).$tempFileSuffix;
+ Assert( !FileExists($fileName),"Temporary file $fileName already exists");
+ return $fileName;
+}
+
+sub LinesFromFile($)
+{
+ my $file = shift;
+ my $ist = OpenFilesForReading($file);
+ my $array = [];
+ while (my $line = $ist->getline)
+ {
+ chomp $line;
+ push @$array, $line;
+ }
+ return $array;
+}
+
+
+sub OpenFilesForReading
+{
+ 0 < @ARG or WrongNumArgsError();
+
+ my @aHandles = ();
+ my $fh;
+ my $filename;
+
+ for $filename (@ARG)
+ {
+ $fh = new FileHandle($filename, "r")
+ or confess "Exception: Could not open $filename for reading: $ERRNO\n";
+ push @aHandles, $fh;
+ }
+
+ return wantarray? @aHandles : $aHandles[0];
+}
+
+sub OpenFilesForAppending
+{
+ 0 < @ARG or WrongNumArgsError();
+
+ my @aFh = ();
+ for my $file (@ARG)
+ {
+ my $fh = new FileHandle($file, "a")
+ or confess "Exception: Could not open $file for reading: "
+ .$ERRNO."\n";
+ push @aFh, $fh;
+ }
+ return wantarray ? @aFh : $aFh[0];
+}
+
+sub OpenFilesForWriting
+{
+ 0 < @ARG or WrongNumArgsError();
+
+ my @aHandles = ();
+ my $fh;
+ my $filename;
+
+ for $filename (@ARG)
+ {
+ $fh = new FileHandle($filename, "w")
+ or confess "Exception: Could not open $filename for writing: $ERRNO\n";
+ push @aHandles, $fh;
+ }
+
+ return wantarray? @aHandles : $aHandles[0];
+}
+
+sub SplitPathName($)
+{
+ my $pathname = shift;
+
+ $pathname =~ m/(.*?)([^\/]+)$/;
+ my $path = $1 || "./";
+ my $basename = $2;
+
+ my $suffix = $basename =~ s/\.([^\.]*)$// # if has extension
+ ? $1 # assign ext
+ : null; # otherwise ext = null
+
+ return $basename, $path, $suffix;
+}
+
+sub GetFileSuffix($)
+{
+ return (SplitPathName(shift))[2];
+}
+
+sub GetFileBasename($)
+{
+ return (SplitPathName(shift))[0];
+}
+
+sub GetFilePath($)
+{
+ return (SplitPathName(shift))[1];
+}
+
+sub println
+{
+ if (ref $ARG[0])
+ {
+ my $ost = shift;
+ $ost->print(@ARG, "\n");
+ }
+ else
+ {
+ print @ARG, "\n";
+ }
+}
+
+
+#====
+# END
+#==============================================================================
+true;
+
+__END__
+
+=head1 Module
+
+FileIoHelper
+
+=head1 Description
+
+Miscellaneous functions for files.
+
+=head1 Functions
+
+=over
+
+=item CloseFiles
+
+ Function:
+ Closes all fhs given as arguments
+ WARNING:
+ CRE if there is an error on closing any stream
+ ARGUMENTS:
+ <@> array of streams
+ RETURN:
+ void
+
+=item ColsFromFile
+
+ Function: convenience function
+ Returns arrays representing the specified columns from a file
+ ARGUMENTS:
+ 1. <string> filename
+ 2. <\@ int> list of column numbers (starting 0..)
+ RETURN:
+ if scalar
+ 1. <\@ string> first column specified
+ elsif array
+ <@@ string> array of arrays of specified columns
+
+=item ColsFromStream
+
+ Function:
+ Returns arrays representing the specified columns from a file
+ ARGUMENTS:
+ 1. <ist> stream
+ 2. <\@ int> list of column numbers (starting 0..)
+ RETURN:
+ if scalar
+ 1. <\@ string> first column specified
+ elsif array
+ <@@ string> array of arrays of specified columns
+
+=item FileExists
+
+ Function: wrapper
+ Returns whether file exists. Readable wrapper for "-e $filename"
+ ARGUMENTS:
+ 1. <string> filename
+ RETURN:
+ <boolean> true iff file exists
+
+=item GetStdError
+
+ Function: accessor
+ Returns stream version of STDERR
+ RETURN:
+ 1. <ost> standard error stream
+ SEE ALSO:
+
+=item GetStdIn
+
+ Function: accessor
+ Returns stream version of STDIN
+ RETURN:
+ 1. <ist> standard input stream
+
+=item GetStdOut
+
+ Function: accessor
+ Returns stream version of STDOUT
+ RETURN:
+ 1. <ost> standard output stream
+
+=item GetTempFileName
+
+ Function:
+ Returns unique filename suitable for a temporary file. The name includes
+ the process id number of the program and a count.
+ RETURN:
+ 1. <string> filename
+
+=item OpenFilesForAppending
+
+ Function:
+ Safely opens files for appending, returning ost objects
+ WARNING:
+ CRE if cannot open file
+ ARGUMENTS:
+ <@ strings> array of filenames
+ RETURN:
+ <@ ost> array of streams
+
+=item OpenFilesForReading
+
+ Function:
+ Safely opens files for reading, returning ist objects
+ WARNING:
+ CRE if cannot open file
+ ARGUMENTS:
+ <@ strings> array of filenames
+ RETURN:
+ <@ ist> array of streams
+
+=item OpenFilesForWriting
+
+ Function:
+ Safely opens files for writing, returning ost objects
+ WARNING:
+ CRE if cannot open file
+ ARGUMENTS:
+ <@ strings> array of filenames
+ RETURN:
+ <@ ost> array of streams
+
+=item SplitPathName
+
+ Function:
+ Splits a pathname into the basename, path and suffix. Eg,
+ for "/home/bsm/bilbo/flange.txt" :
+ basename: "flange"
+ path: "/home/bsm/bilbo/"
+ suffix: "txt"
+ Does the job File::Basename::fileparse is meant to do.
+ WARNING:
+ May not work properly for "/public/fred/.." type names.
+ ARGUMENTS:
+ 1. <string> filename
+ RETURN:
+ 1. <string> basename
+ 2. <string> path
+ 3. <string> suffix
+
+=item println
+
+ Function:
+ Same as CORE::print but automatically adds on an endline character.
+ ARGUMENTS:
+ <@ ?> things to print
+ RETURN:
+ void
diff --git a/bagpipe/lib/Hash.pm b/bagpipe/lib/Hash.pm
new file mode 100644
index 0000000..3259cba
--- /dev/null
+++ b/bagpipe/lib/Hash.pm
@@ -0,0 +1,168 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package Hash;
+use strict;
+use English;
+use StdDefs;
+use Assert;
+use Array;
+use VarType qw(IsArrayRef);
+
+
+sub Hash::Apply($$)
+{
+ my($hash, $func) = @ARG;
+ my $out = {};
+ while (my($key, $value)=each %$hash)
+ {
+ $out->{$key} = $func->($value);
+ }
+ return $out;
+}
+
+sub Copy
+{
+ 1 == @ARG || 2 == @ARG or WrongNumArgsError();
+
+ my $hSrc = $ARG[0];
+ my $hDest = 2 == @ARG ? $ARG[1] : {};
+
+ Assert(defined($hSrc), "source not defined");
+ Assert(defined($hDest), "dest not defined");
+
+
+ %$hDest = %$hSrc;
+
+ return $hDest;
+}
+
+sub Equals($$)
+{
+ my($hash1, $hash2) = @ARG;
+
+ # cmp sizes
+ if ((scalar keys %$hash1) != (scalar keys %$hash2))
+ {
+ return false;
+ }
+ # cmp contents
+ while (my($key1, $value1) = each %$hash1)
+ {
+ unless (exists $hash2->{$key1} and $hash2->{$key1} eq $value1)
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
+sub Invert($)
+{
+ my ($hOriginal) = @ARG;
+
+ my $hInverted;
+
+ %{$hInverted} = reverse %$hOriginal;
+
+ return $hInverted;
+}
+
+sub IsEmpty($)
+{
+ return 0 == scalar keys %{$ARG[0]};
+}
+
+sub Join
+{
+ my(@aHashes) = @ARG;
+
+ my $hJoined = {};
+
+ for my $hash (@aHashes)
+ {
+ while (my($key, $value) = each(%$hash))
+ {
+ $hJoined->{$key} = $value;
+ }
+ }
+ return $hJoined;
+}
+
+sub New($$)
+{
+ if (IsArrayRef($ARG[1]))
+ {
+ return New_from2Arrays(@ARG);
+ }
+ else
+ {
+ return New_from1Array(@ARG);
+ }
+}
+
+sub New_from1Array
+{
+ my($aKeys, $value) = @ARG;
+ return New_from2Arrays($aKeys,
+ Array::New(scalar(@$aKeys), $value));
+}
+
+sub New_from2Arrays
+{
+ my ($aKeys, $aValues) = @ARG;
+ my $hash;
+ @$hash{@$aKeys} = @$aValues;
+ return $hash;
+}
+
+sub RandomKey($)
+{
+ my $hash = shift;
+ my @keys = keys %$hash;
+ return Array::RandomElement(\@keys);
+}
+
+sub Size($)
+{
+ return scalar keys %{$ARG[0]};
+}
+
+sub Hash::ToArrays($)
+{
+ my $hash = shift;
+
+ my $aKeys = [];
+ my $aValues = [];
+ while (my($key, $value)=each %$hash)
+ {
+ push @$aKeys, $key;
+ push @$aValues, $value;
+ }
+ return ($aKeys, $aValues);
+}
+
+sub Hash::Transform($$$)
+{
+ my($hIn, $hOut, $func) = @ARG;
+ while (my($key, $value)=each %$hIn)
+ {
+ $hOut->{$key} = $func->($value);
+ }
+ return $hOut;
+}
+
+true;
diff --git a/bagpipe/lib/NamedArgs.pm b/bagpipe/lib/NamedArgs.pm
new file mode 100644
index 0000000..65f898f
--- /dev/null
+++ b/bagpipe/lib/NamedArgs.pm
@@ -0,0 +1,113 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package NamedArgs;
+use base qw(Exporter);
+use Exporter;
+use English;
+use StdDefs;
+use Assert;
+use strict;
+use Number;
+
+sub NamedArgs::Parse($$$)
+{
+ my($aSubArgs, $numPlainArgs, $hSpec) = @ARG;
+
+ #---------------------------
+ # Check number of args given
+ #---------------------------
+
+ my $numSubArgs = scalar(@$aSubArgs);
+
+ # check minimum number of args given
+ my $numReqArgs = $numPlainArgs;
+ if (exists $hSpec->{'required'})
+ {
+ $numReqArgs += 2*scalar(@{$hSpec->{'required'}});
+ }
+ unless ($numSubArgs >= $numReqArgs)
+ {
+ WrongNumArgsError();
+ }
+
+ # check number of args is sensible
+ if (
+ (Number::IsEven($numPlainArgs) and Number::IsOdd($numSubArgs))
+ or
+ (Number::IsOdd($numPlainArgs) and Number::IsEven($numSubArgs))
+ )
+ {
+ WrongNumArgsError();
+ }
+
+ # check maximum number of args
+ if (exists $hSpec->{'optional'})
+ {
+ my $numMaxArgs = $numReqArgs
+ + 2*scalar keys %{$hSpec->{'optional'}};
+ if ($numSubArgs > $numMaxArgs)
+ {
+ WrongNumArgsError();
+ }
+ }
+
+ #------------------
+ # Separate out args
+ #------------------
+
+ # take out unnamed args
+ my @aPlainArgs = splice @$aSubArgs, 0, $numPlainArgs;
+ my %hSubArgs = @$aSubArgs;
+
+ # prepare specification
+ my %hRequired = ();
+ @hRequired{@{$hSpec->{'required'}}} = () if exists $hSpec->{'required'};
+
+ my %hOptional = (exists $hSpec->{'optional'})
+ ? %{$hSpec->{'optional'}}
+ : ();
+
+ # fill args, check for illegals
+ while (my($name, $arg) = each %hSubArgs)
+ {
+ if (exists $hRequired{$name})
+ {
+ $hRequired{$name} = $arg;
+ }
+ elsif (exists $hOptional{$name})
+ {
+ $hOptional{$name} = $arg;
+ }
+ else
+ {
+ Fatal("Invalid argument \'$name\' to method\n");
+ }
+ }
+ # check for unfilled requireds
+ while (my($name, $val) = each %hRequired)
+ {
+ if (not defined($val))
+ {
+ Fatal("Must specify value for required argument \'$name\'");
+ }
+ }
+
+ return @aPlainArgs, {%hRequired, %hOptional};
+}
+1;
+
+
diff --git a/bagpipe/lib/Number.pm b/bagpipe/lib/Number.pm
new file mode 100644
index 0000000..5b00550
--- /dev/null
+++ b/bagpipe/lib/Number.pm
@@ -0,0 +1,34 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+#!./perl
+package Number;
+use strict;
+use English;
+
+sub IsEven($)
+{
+ return IsOdd($ARG[0]) ? 0 : 1;
+}
+
+sub IsOdd($)
+{
+ return ($ARG[0] % 2) ? 1 : 0;
+}
+
+1;
+
+
diff --git a/bagpipe/lib/StdDefs.pm b/bagpipe/lib/StdDefs.pm
new file mode 100644
index 0000000..72bf307
--- /dev/null
+++ b/bagpipe/lib/StdDefs.pm
@@ -0,0 +1,28 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package StdDefs;
+use strict;
+use vars qw( @ISA @EXPORT );
+use Exporter;
+
+ at ISA = qw( Exporter );
+ at EXPORT = qw( true false null );
+
+use constant true => 1;
+use constant false => 0;
+use constant null => '';
+
diff --git a/bagpipe/lib/String.pm b/bagpipe/lib/String.pm
new file mode 100644
index 0000000..dd0f1ee
--- /dev/null
+++ b/bagpipe/lib/String.pm
@@ -0,0 +1,90 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package String;
+use strict;
+use English;
+use StdDefs;
+use Exporter;
+use Array;
+use base qw(Exporter);
+use vars qw(@EXPORT_OK);
+
+ at EXPORT_OK = qw(
+ Trim
+ Chomp
+ Columnify
+ CountChar
+ UniqueChars
+ );
+
+sub Chomp($)
+{
+ my $s = shift;
+ $s =~ s/[\r\n]+$//;
+ return $s;
+}
+
+sub Columnify
+{
+ my($colw, @aData) = @ARG;
+
+ return sprintf( "%-${colw}.${colw}s" x scalar(@aData), @aData );
+}
+
+sub CountChar($$)
+{
+ my($string, $char) = @ARG;
+
+ my $count = 0;
+ $count++ while $string =~ m/$char/go;
+ return $count;
+}
+
+sub Trim($)
+{
+ my $text = shift;
+ $text =~ s/^\s+//;
+ $text =~ s/\s+$//;
+ return $text;
+}
+
+sub TrimArray($)
+{
+ my $array = shift;
+ for my $text (@$array)
+ {
+ $text =~ s/^\s+//;
+ $text =~ s/\s+$//;
+ }
+ return $array;
+}
+
+sub UniqueChars($)
+{
+ my $text = shift;
+ return Array::Unique([split m//, $text]);
+}
+
+sub OrderedStrcat($$$)
+{
+ my($a, $b, $token) = @ARG;
+
+ return $a lt $b ? $a.$token.$b : $b.$token.$a;
+}
+
+
+true;
diff --git a/bagpipe/lib/VarType.pm b/bagpipe/lib/VarType.pm
new file mode 100644
index 0000000..266835e
--- /dev/null
+++ b/bagpipe/lib/VarType.pm
@@ -0,0 +1,67 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+package VarType;
+use strict;
+use English;
+use StdDefs;
+use Exporter;
+use base qw(Exporter);
+use vars qw( @EXPORT_OK );
+
+ at EXPORT_OK = qw(
+ IsArrayRef
+ IsHashRef
+ IsNumber
+ AsArray
+ );
+
+sub IsArrayRef($);
+sub IsInt($);
+sub IsHashRef($);
+sub IsNumber($);
+
+sub AsArray($)
+{
+ my $var = shift;
+
+ return (IsArrayRef($var))
+ ? $var
+ : [$var];
+}
+
+sub IsArrayRef($) { (ref($ARG[0]) =~ m/ARRAY/) ? true : false}
+sub IsInt($) { ($ARG[0] =~ m/^\d+$/) ? true : false}
+sub IsHashRef($) { (ref($ARG[0]) =~ m/HASH/) ? true : false}
+sub IsNumber($) {(defined ParseNumber($ARG[0])) ? true : false}
+
+sub ParseNumber($)
+{
+ use POSIX qw(strtod);
+ my $str = shift;
+ $str =~ s/^\s+//;
+ $str =~ s/\s+$//;
+ $ERRNO = 0;
+ my($num, $unparsed) = strtod($str);
+ if (($str eq '') || ($unparsed != 0) || $ERRNO)
+ {
+ return undef;
+ }
+ else
+ {
+ return $num;
+ }
+}
diff --git a/bagpipe/lib/fitgev.R b/bagpipe/lib/fitgev.R
new file mode 100644
index 0000000..5e2dd96
--- /dev/null
+++ b/bagpipe/lib/fitgev.R
@@ -0,0 +1,63 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+options(warn=1)
+
+inputdir <- cmdline.option("indir", "./")
+null.files <- paste(inputdir, sep="/", cmdline.strings("infiles"))
+tailprobs <- as.numeric(cmdline.strings("tailprobs"))
+
+# read data
+lod.results <- list()
+modelcmp.results <- list()
+print(null.files)
+for (file in null.files)
+{
+ group.data <- read.delim(file)
+ for (i in 1:nrow(group.data))
+ {
+ s <- group.data$scan.number[i]
+ lod.results[[s]] <- c(unlist(lod.results[s]), group.data$best.LOD[i])
+ modelcmp.results[[s]] <- c(unlist(modelcmp.results[s]), group.data$best.modelcmp[i])
+ }
+}
+if (1!=length(unique(sapply(lod.results, length))) | 1!=length(unique(sapply(lod.results, length))))
+{
+ stop("Cannot calculate thresholds when the number of scans vary among locus groups\n")
+}
+overall.modelcmp <- sapply(modelcmp.results, max, na.rm=TRUE)
+overall.lod <- sapply(lod.results, max, na.rm=TRUE)
+# fit GEV
+fit.modelcmp <- try(fit.gev(overall.modelcmp, tailprobs))
+fit.lod <- try(fit.gev(overall.lod, tailprobs))
+if (caught.error(fit.modelcmp) | caught.error(fit.lod))
+{
+ stop("GEV fit failed: could be too few data points (eg, too few nullscans)\n")
+}
+
+gev.out <- cbind(c("LOD", "modelcmp"), c(length(overall.lod), length(overall.modelcmp)), rbind(fit.lod$gev, fit.modelcmp$gev))
+colnames(gev.out)[1:2] <- c("score","num.obs")
+write.delim(gev.out, cmdline.option("gevfile"))
+
+thresh.out <- cbind(fit.lod$thresholds, fit.modelcmp$thresholds$quantile)
+colnames(thresh.out) <- c("upper.tail.prob", "LOD", "modelcmp")
+write.delim(thresh.out, cmdline.option("threshfile"))
diff --git a/bagpipe/lib/nullscangenome.R b/bagpipe/lib/nullscangenome.R
new file mode 100644
index 0000000..6dea5b2
--- /dev/null
+++ b/bagpipe/lib/nullscangenome.R
@@ -0,0 +1,55 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+options(warn=1)
+# TODO: make random seed the same for all chromosomes
+
+## start clock
+(my.start <- Sys.time())
+
+## null or perm scan
+
+config <- bagpipe.read.configfile(cmdline.option("configfile"))
+h <- happy.load.genome(configfile.string(config, "genome.cache.dir"))
+h <- happy.init.reserve(h, memory.limit.Mb=cmdline.numeric("memory"), auto.reserve=TRUE)
+loci <- bagpipe.extract.loci(h, cmdline.option("locus.group"))
+
+do.scan(h,
+ config = config,
+ loci = loci,
+ output.dir = cmdline.option("outdir"),
+ output.file = cmdline.option("outfile"),
+ scan.type = cmdline.option("scan.type",
+ allowed.values=c(
+ "nullsim.scan",
+ "permute.scan",
+ "nullsimpermute.scan")),
+ verbose = TRUE,
+ cpus = cmdline.integer("cpus",1))
+
+## stop clock
+(my.end <- Sys.time())
+my.time <- difftime(my.end, my.start)
+cat("Runtime ", my.time, attr(my.time, "units"), "\n")
+
+
+
diff --git a/bagpipe/lib/plotscans.R b/bagpipe/lib/plotscans.R
new file mode 100644
index 0000000..c782d3e
--- /dev/null
+++ b/bagpipe/lib/plotscans.R
@@ -0,0 +1,119 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+options(warn=1)
+
+config <- bagpipe.read.configfile(cmdline.option("configfile"))
+h <- happy.load.genome(configfile.string(config, "genome.cache.dir"))
+
+inputdir <- cmdline.strings("indir", default="./")
+infiles <- paste(inputdir, sep="/", cmdline.strings("infiles"))
+
+plot.file <- cmdline.option("outfile")
+plot.height <- cmdline.integer("plot.height")
+plot.width <- cmdline.integer("plot.width")
+plot.type <- cmdline.option("plot.type", default="chromosomes")
+plot.box <- cmdline.flag("plot.box")
+thresholds <- as.numeric(cmdline.strings("thresholds", stop.on.fail=FALSE))
+scale <- cmdline.option("plot.scale")
+line.color <- cmdline.option("line.color", default="black")
+main.title <- cmdline.option("title",
+ default= paste(sep="\n",configfile.string(config, "scan.formula.null"), configfile.string(config, "scan.formula.test")))
+
+qtls <- NULL
+if (cmdline.has.option("qtls"))
+{
+ qtls <- read.delim(cmdline.option("qtls"))
+}
+
+gscan.data <- load.gscan(files=infiles, dir=inputdir)
+score.type <- cmdline.option("locus.score.type")
+score.name <- switch(score.type,
+ LOD= "LOD",
+ modelcmp=as.character(unique(na.omit(gscan.data$genome$modelcmp.type))) )
+
+all.loci <- gscan.data$genome$locus
+all.scores <- gscan.data$genome[,score.type]
+max.score <- max(c(all.scores, thresholds), na.rm=TRUE)
+
+flag.missing.y.line <- function(...){}
+if (cmdline.flag("flag.missing"))
+{
+ flag.missing.y.line <- function(x,ylim,...){ mtext("*", side=1, at=mean(x), line=-1)}
+}
+
+pdf(plot.file, height=plot.height, width=plot.width)
+if (plot.type=="chromosomes")
+{
+ for (chr in happy.list.chromosomes(h))
+ {
+ i.chr <- chr==happy.get.chromosome(h, all.loci)
+ if (!any(i.chr)) next
+ loci <- all.loci[i.chr]
+ scores <- all.scores[i.chr]
+ happy.plot.intervals(h, scale=scale, ylim = c(0, 1.1*max.score),
+ loci.lim=loci[c(1,length(loci))],
+ main = main.title,
+ ylab = switch(score.name,
+ logP=expression(-log[10]*P),
+ LOD="LOD score",
+ score.name),
+ xlab=paste("Position on chromosome ", chr, " (",scale,")",sep=""),
+ loci=loci,
+ y = scores,
+ col=line.color,
+ axis.x.extend=0.05,
+ flag.missing.y.line=flag.missing.y.line,
+ draw.edge.support=FALSE)
+ if (!is.null(thresholds))
+ {
+ abline(h=thresholds, col=line.color,
+ lty=rep(2:5,length.out=length(thresholds)))
+ }
+ if (plot.box) box()
+ }
+}
+
+if (plot.type=="genome")
+{
+ happy.plot.intervals(h, scale=scale, ylim = c(0, 1.1*max.score),
+ main = main.title,
+ ylab = switch(score.name,
+ logP=expression(-log[10]*P),
+ LOD="LOD score",
+ score.name),
+ xlab=paste("Chromosome (",scale," scale)",sep=""),
+ loci=all.loci,
+ y = all.scores,
+ col=line.color,
+ flag.missing.y.line=flag.missing.y.line,
+ draw.edge.support=FALSE)
+ if (!is.null(thresholds))
+ {
+ abline(h=thresholds, col=line.color,
+ lty=rep(2:5,length.out=length(thresholds)))
+ }
+ if (plot.box) box()
+}
+dev.off()
+
+
diff --git a/bagpipe/lib/posboot.R b/bagpipe/lib/posboot.R
new file mode 100644
index 0000000..7beae43
--- /dev/null
+++ b/bagpipe/lib/posboot.R
@@ -0,0 +1,87 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+
+options(warn=1)
+cat(commandArgs(), "\n")
+## start clock
+(my.start <- Sys.time())
+
+#--- BODY ---
+
+config <- bagpipe.read.configfile(cmdline.option("configfile"))
+h <- happy.load.genome(configfile.string(config, "genome.cache.dir"))
+num.boots <- cmdline.integer("num.posboots")
+bootfile <- cmdline.option("posboot.results")
+
+locus.range <- cmdline.strings("posboot.between")
+loci=bagpipe.define.posboot.loci(h, locus.range)
+h <- happy.init.reserve(h, memory.limit.Mb=cmdline.numeric("memory"), auto.reserve=TRUE)
+boot.results=bagpipe.posboot.scan(h, loci, num.boots, bootfile, save.every=1)
+boot.summary=make.posboot.summary(h, loci, boot.results)
+write.delim(file=cmdline.string("posboot.summary"), boot.summary)
+
+ci50=rbind(
+ make.posboot.summary.ci(h, boot.summary, prob=0.50, score="LOD"),
+ make.posboot.summary.ci(h, boot.summary, prob=0.50, score="modelcmp"))
+ci80=rbind(
+ make.posboot.summary.ci(h, boot.summary, prob=0.80, score="LOD"),
+ make.posboot.summary.ci(h, boot.summary, prob=0.80, score="modelcmp"))
+ci95=rbind(
+ make.posboot.summary.ci(h, boot.summary, prob=0.95, score="LOD"),
+ make.posboot.summary.ci(h, boot.summary, prob=0.95, score="modelcmp"))
+write.delim(file=cmdline.string("posboot.quicksummary"), rbind(ci50,ci80,ci95))
+
+if (cmdline.has.option("posboot.plot"))
+{
+ plot.scale <- cmdline.string("plot.scale", default="cM")
+ loci.lim <- c(loci[1], tail(loci,1))
+ chr <- happy.get.chromosome(h, loci[1])
+
+ pdf(file=cmdline.string("posboot.plot"), width=6, height=5)
+
+ for (score in c("LOD", "modelcmp"))
+ {
+ y=boot.summary[,paste("times.max.", sep="", score)]
+ max.y=max(y)
+ qy=boot.summary[,paste("cumfrac.max.", sep="", score)]
+ qy=ifelse(qy>0.5, 1-qy, qy)
+
+ i50=as.integer(ci50[ci50$score==score,c("idx.start","idx.end")])
+ i80=as.integer(ci80[ci80$score==score,c("idx.start","idx.end")])
+ i95=as.integer(ci95[ci95$score==score,c("idx.start","idx.end")])
+
+ fill=rep("white", length(y))
+ fill[i95[1]:i95[2]] = "gray80"
+ fill[i80[1]:i80[2]] = "gray50"
+ fill[i50[1]:i50[2]] = "black"
+
+ happy.plot.intervals(h, loci=loci, y=y,
+ loci.lim=loci.lim,
+ ylim=c(0,max.y), fill=fill,
+ ylab="Frequency at peak top",
+ xlab=paste("Position on chromosome ", chr," (",plot.scale,")",sep=""),
+ main=paste(score, "maxima from", sum(y), "positional bootstraps"),
+ scale=plot.scale)
+ }
+ dev.off()
+}
diff --git a/bagpipe/lib/scangenome.R b/bagpipe/lib/scangenome.R
new file mode 100644
index 0000000..830d363
--- /dev/null
+++ b/bagpipe/lib/scangenome.R
@@ -0,0 +1,50 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+
+options(warn=1)
+cat(commandArgs(), "\n")
+## start clock
+(my.start <- Sys.time())
+
+#--- BODY ---
+
+config <- bagpipe.read.configfile(cmdline.option("configfile"))
+h <- happy.load.genome(configfile.string(config, "genome.cache.dir"))
+loci <- bagpipe.extract.loci(h, cmdline.option("locus.group"))
+
+#---main scan code
+
+do.scan(h,
+ config = config,
+ loci = loci,
+ output.dir = cmdline.option("outdir"),
+ output.file = cmdline.option("outfile"),
+ scan.type = "scan",
+ verbose = TRUE)
+
+#--- /BODY ---
+
+## stop clock
+(my.end <- Sys.time())
+my.time <- difftime(my.end, my.start)
+cat(my.time, attr(my.time, "units"), "\n")
diff --git a/bagpipe/lib/scanloci.R b/bagpipe/lib/scanloci.R
new file mode 100644
index 0000000..5d24c6c
--- /dev/null
+++ b/bagpipe/lib/scanloci.R
@@ -0,0 +1,57 @@
+#Copyright 2012 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe.backend)
+if(cmdline.flag("devel"))
+{
+ detach("package:bagpipe.backend")
+ library(bagpipe)
+}
+
+options(warn=1)
+cat(commandArgs(), "\n")
+## start clock
+(my.start <- Sys.time())
+
+#--- BODY ---
+
+config = bagpipe.read.configfile(cmdline.option("configfile"))
+h = happy.load.genome(configfile.string(config, "genome.cache.dir"))
+loci = cmdline.strings("loci")
+loci.files = cmdline.strings("loci.files")
+loci.files = file.path(cmdline.string("outdir"),loci.files)
+
+if (!all(happy.has.markers(h, loci)))
+{
+ bagpipe.input.error("Cannot save locus fit data for unknown loci:",
+ paste(collapse=", ", loci[!happy.has.markers(h, loci)]))
+}
+save.at.loci=as.list(loci.files,all.names=TRUE)
+names(save.at.loci)=loci
+
+do.scan(h,
+ config = config,
+ loci = loci,
+ scan.type = "scan",
+ save.at.loci = save.at.loci,
+ verbose = TRUE)
+
+#--- /BODY ---
+
+## stop clock
+(my.end <- Sys.time())
+my.time <- difftime(my.end, my.start)
+cat(my.time, attr(my.time, "units"), "\n")
diff --git a/bagpipe/lib/write_nomenu_config.R b/bagpipe/lib/write_nomenu_config.R
new file mode 100644
index 0000000..489ef0d
--- /dev/null
+++ b/bagpipe/lib/write_nomenu_config.R
@@ -0,0 +1,25 @@
+#Copyright 2010 William Valdar
+#This file is part of bagpipe.
+
+#bagpipe 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.
+
+#bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
+
+library(bagpipe)
+options(warn=2)
+
+for (file in list.files("~/WorkHome/code/Rbagphenotype/devel/", pattern=".*R$", full.names=TRUE)) source(file)
+
+phenotype <- cmdline.option("phenotype")
+config <- make.phenotype.config(phenotype, cmdline.option("configfile"))
+config$menu.file <- NULL
+write.configfile(config, file = cmdline.option("outfile") )
diff --git a/debian/bagpipe.1 b/debian/bagpipe.1
deleted file mode 100644
index 73219c7..0000000
--- a/debian/bagpipe.1
+++ /dev/null
@@ -1,150 +0,0 @@
-.\" DO NOT MODIFY THIS FILE! It was generated by help2man 1.41.1.
-.TH bagpipe "1" "April 2013" "bagpipe" "User Commands"
-.SH NAME
-bagpipe \- a pipeline for genetic association using HAPPY genome data
-.SS "Usage:"
-.IP
-bagpipe configfile [\-action(s)] [\-option(s)=argument]
-.PP
-actions: scan, perm, perm_collate, nullsim, nullsim_collate, plot
-.PP
-options: locus_groups, outdir, scandir, plot_height, plot_width,
-plot_scale, plot_type, plot_score, plot_tailprob memory, overwrite
-.PP
-Type 'bagpipe \fB\-h\fR' for summary usage
-.PP
-\&'bagpipe \fB\-\-help\fR' for detailed usage
-.PP
-\&'bagpipe \fB\-\-examples\fR' for examples
-.PP
-Requirements:
-Config file: a file of configuration parameters
-.IP
-HAPPY genome cache: a directory containing cached HAPPY matrices
-Phenotype file: a tab delimited table with phenotype and covariate
-.IP
-values for each animal. The table must include
-.IP
-1) a column SUBJECT.NAME giving the id of the animal in
-the corresponding HAPPY data;
-.IP
-2) a column with the name of the phenotype as listed in the
-Phenotype column of the model menu file;
-.IP
-3) columns corresponding to any covariates used in the
-model menu formulae.
-.PP
-Key:
-.IP
-n \- an integer
-.IP
-d \- a decimal
-.IP
-s \- a string
-.IP
-b \- 1 for true and 0 for false.
-.IP
-xL \- a comma separated list of type x
-.IP
-locus.group \- either the name of a chromosome (eg, 1 or X) or a file of marker names.
-.SS "Actions:"
-.TP
-scan
-\- Genome scans
-.TP
-perm
-\- Permutations
-.TP
-perm_collate
-\- Calculate permutation thresholds
-.TP
-nullsim
-\- Null simulations (not yet implemented)
-.TP
-nullsim_collate
-\- Calculate null simulation thresholds (not yet implemented)
-.TP
-plot
-\- Plot QTL scans to pdf file.
-.TP
-scanloci=sL
-\- [Experimental] Scan individual loci, saving all fit information
-.TP
-posboot=sL
-\- [Experimental] Perform k positional bootstraps between marker
-.IP
-intervals m1 and m2, where s=m1,m2,k
-.SS "Settings:"
-.TP
-outdir
-\- Directory for output. Default is ./
-.TP
-scandir
-\- Directory for output of scan data files within outdir (./)
-.TP
-locus_groups=sL
-\- Run analyses on certain locus groups only [Default: do all in configfile]
-.TP
-save.at.loci=sL
-\- [Experimental] Save model fit when scanning at specified loci.
-.TP
-tailprobs=dL
-\- Calculate nullscan tail probabilities (ie, significance thresholds).
-.IP
-Eg, 0.05,0.01. [Default nL=0.001,0.01,0.05,0.1,0.2]
-.SS "Plotting:"
-.TP
-plot_box
-\- Draw box around each plot
-.TP
-plot_scale=s
-\- Scale used for plot summary s={cM,bp,Mb}
-.TP
-plot_score=s
-\- Specify locus score to plot s={LOD,modelcmp}
-.TP
-plot_type=s
-\- s = {chromosomes,genome}
-.IP
-plot_tailprob=sL \- Include threshold lines as calculated by the *_collate options.
-.IP
-Each string s should have format <nullscantype>_<d>, where
-nullscantype={perm,nullsim} and <d> is the tailprob looked up in the
-corresponding thresh file. Eg, s=perm_0.05 plots the 0.05 threshold
-calculated by permutation.
-.TP
-plot_height=d
-\- Notional height of plotsummary output in inches
-.TP
-plot_width=d
-\- Notional width of plotsummary output in inches
-.SS "Performance:"
-.TP
-cleanup
-\- Delete Rout files that do not report errors or warnings
-.TP
-dryrun
-\- Print commands but do not execute them
-.TP
-memory=d
-\- Store marker matrices in memory upto n Mb [default n=0]
-.TP
-nomenuconfig
-\- Write the complete config file for an analysis
-.IP
-ie, incorporating all the information from the menu file
-but minus the menu file
-.TP
-overwrite=b
-\- Overwrite existing analysis files [Default overwrite=0]
-.SS "Help:"
-.TP
-examples
-\- Give examples
-.TP
-example_config
-\- Write an example config file to config.example
-.SS "Web"
-.TP
-More information can be found at:
-http://valdarlab.unc.edu/software/bagpipe/_build/html/bagpipe.html
diff --git a/debian/bagpipe.install b/debian/bagpipe.install
deleted file mode 100644
index 7bb8044..0000000
--- a/debian/bagpipe.install
+++ /dev/null
@@ -1,3 +0,0 @@
-debian/script/bagpipe usr/bin
-bagpipe/bagpipe.pl usr/share/bagpipe
-bagpipe/lib/* usr/share/perl5/Bagpipe
diff --git a/debian/bagpipe.manpages b/debian/bagpipe.manpages
deleted file mode 100644
index 752cbb0..0000000
--- a/debian/bagpipe.manpages
+++ /dev/null
@@ -1 +0,0 @@
-debian/bagpipe.1
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index dc8b2c6..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,6 +0,0 @@
-bagpipe (2012.02.15-1) UNRELEASED; urgency=low
-
- * initial version (Closes: #XXXXXX)
-
- -- Thorsten Alteholz <debian at alteholz.de> Tue, 09 Apr 2013 18:00:07 +0200
-
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index f599e28..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-10
diff --git a/debian/control b/debian/control
deleted file mode 100644
index 1df20f2..0000000
--- a/debian/control
+++ /dev/null
@@ -1,39 +0,0 @@
-Source: bagpipe
-Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
-Uploaders: Thorsten Alteholz <debian at alteholz.de>
-Section: science
-Priority: optional
-Build-Depends: debhelper (>= 10)
-Standards-Version: 3.9.8
-Vcs-Browser: http://anonscm.debian.org/viewvc/debian-med/trunk/packages/bagpipe/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/bagpipe/trunk/
-Homepage: http://valdarlab.unc.edu/software/bagpipe/_build/html/bagpipe.html
-
-Package: bagpipe
-Architecture: all
-Depends: ${shlibs:Depends},
- ${misc:Depends},
- ${perl:Depends},
- r-other-valdar-bagpipe.backend
-Description: genomewide LD mapping
- Bagpipe is a program for performing genomewide linkage disequilibrium
- mapping of quantitative trait loci in populations whose genome structure
- can be accommodated in the HAPPY framework [Mott00]. This includes most
- diploid crosses where the founders of the individuals have known genotypes.
- .
- * Bagpipe is a simplified and streamlined version of Bagphenotype that
- does not currently include resample model averaging (RMA) capabilities.
- * Bagpipe can help fit single locus regression models (with or without
- random effects) to marker intervals whose genetic ancestry is inferred
- using the HAPPY software.
- * Bagpipe cannot help you decide what is a sensible model to fit.
- * Bagpipe does not currently accommodate populations with significant
- population structure, except through the specification of simple random
- intercepts based on unpatterned covariance matrices.
- * Bagpipe is named after the Scottish wind instrument "the bagpipes" and
- after Bagphenotype, which in turn was a PIPEline for BAGging-based
- multiple QTL analysis of phenoTYPEs. Bagphenotype was in turn based
- on software written by Richard Mott and William Valdar to analyze
- heterogeneous stock mice in [Valdar06].
- * Bagpipe is experimental software, is provided free of charge subject to
- copyleft restrictions, and comes with no guarantees whatsoever.
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index 91260ff..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,27 +0,0 @@
-Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: bagpipe
-Source: http://valdarlab.unc.edu/software/bagpipe/_build/html/bagpipe.html
-
-Files: *
-Copyright: 2012 William Valda
-License: GPL3+
- bagpipe 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.
- .
- bagpipe 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 bagpipe. If not, see <http://www.gnu.org/licenses/>.
- .
- On Debian systems, the complete text of the GNU General
- Public License version 2 can be found in "/usr/share/common-licenses/GPL-3".
-
-Files: debian/*
-Copyright: 2013 Thorsten Alteholz <debian at alteholz.de>
-License: GPL3+
- see above
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 6d852ff..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/usr/bin/make -f
-
-# DH_VERBOSE := 1
-
-%:
- dh $@
-
-override_dh_auto_install:
- dh_auto_install
- ls -lR
-# dh_installchangelogs README
-
-get-orig-source:
- mkdir -p ../tarballs
- uscan --verbose --force-download --destdir=../tarballs
-
diff --git a/debian/script/bagpipe b/debian/script/bagpipe
deleted file mode 100755
index 3f5b7c2..0000000
--- a/debian/script/bagpipe
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/bin/bash
-
-export BAGPIPE_LIBS=/usr/share/perl5/Bagpipe
-/usr/share/bagpipe/bagpipe.pl $@
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/upstream/metadata b/debian/upstream/metadata
deleted file mode 100644
index ee49be8..0000000
--- a/debian/upstream/metadata
+++ /dev/null
@@ -1,22 +0,0 @@
-Reference:
- - Author: Richard Mott and Christopher J. Talbot and Maria G. Turri and Allan C. Collins and Jonathan Flint
- Title: "A method for fine mapping quantitative trait loci in outbred animal stocks"
- Journal: Proc Natl Acad Sci U S A.
- Year: 2000
- Volume: 97
- Number: 23
- Pages: 12649-54
- DOI: 10.1073/pnas.230304397
- PMID: 11050180
- URL: http://www.pnas.org/content/97/23/12649
- - Author: William Valdar and Leah C Solberg and Dominique Gauguier and Stephanie Burnett and Paul Klenerman and William O Cookson and Martin S Taylor and J Nicholas P Rawlins and Richard Mott and Jonathan Flint
- Title: "Genome-wide genetic association of complex traits in heterogeneous stock mice"
- Journal: Nature Genetic
- Year: 2006
- Volume: 38
- Number: 8
- Pages: 879-87
- DOI: 10.1038/ng1840
- PMID: 16832355
- URL: http://www.nature.com/ng/journal/v38/n8/full/ng1840.html
- eprint: http://www.nature.com/ng/journal/v38/n8/pdf/ng1840.pdf
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index 76ff926..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,4 +0,0 @@
-version=3
-opts="uversionmangle=s/_/./g" \
- http://valdarlab.unc.edu/software/bagpipe/_build/html/bagpipe.html http://valdarlab.unc.edu/software/bagpipe/install/bagpipe_([\d_]+)\.tar\.gz
-
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/bagpipe.git
More information about the debian-med-commit
mailing list