From 2ee8d0084bd6ab7bb8f66c32d7fd7b1e32e75ae3 Mon Sep 17 00:00:00 2001 From: Carsten Allefeld Date: Sun, 7 Aug 2016 14:47:13 +0200 Subject: [PATCH] version of 2014-2-18 --- LICENSE | 674 ------------------------------------------ README.md | 1 - README.pdf | Bin 0 -> 32601 bytes README.text | 101 +++++++ contrasts.m | 60 ++++ cvManovaSearchlight.m | 119 ++++++++ cvManova_compute.m | 139 +++++++++ cvManova_precompute.m | 40 +++ inestimability.m | 41 +++ loadDataSPM.m | 121 ++++++++ patchPath.m | 44 +++ runSearchlight.m | 93 ++++++ saveMRImage.m | 66 +++++ signPermutations.m | 35 +++ systemFree.m | 26 ++ 15 files changed, 885 insertions(+), 675 deletions(-) delete mode 100644 LICENSE delete mode 100644 README.md create mode 100644 README.pdf create mode 100644 README.text create mode 100644 contrasts.m create mode 100644 cvManovaSearchlight.m create mode 100644 cvManova_compute.m create mode 100644 cvManova_precompute.m create mode 100644 inestimability.m create mode 100644 loadDataSPM.m create mode 100644 patchPath.m create mode 100644 runSearchlight.m create mode 100644 saveMRImage.m create mode 100644 signPermutations.m create mode 100644 systemFree.m diff --git a/LICENSE b/LICENSE deleted file mode 100644 index 9cecc1d..0000000 --- a/LICENSE +++ /dev/null @@ -1,674 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 3, 29 June 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU General Public License is a free, copyleft license for -software and other kinds of works. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -the GNU General Public License is intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. We, the Free Software Foundation, use the -GNU General Public License for most of our software; it applies also to -any other work released this way by its authors. You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - To protect your rights, we need to prevent others from denying you -these rights or asking you to surrender the rights. Therefore, you have -certain responsibilities if you distribute copies of the software, or if -you modify it: responsibilities to respect the freedom of others. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must pass on to the recipients the same -freedoms that you received. You must make sure that they, too, receive -or can get the source code. And you must show them these terms so they -know their rights. - - Developers that use the GNU GPL protect your rights with two steps: -(1) assert copyright on the software, and (2) offer you this License -giving you legal permission to copy, distribute and/or modify it. - - For the developers' and authors' protection, the GPL clearly explains -that there is no warranty for this free software. For both users' and -authors' sake, the GPL requires that modified versions be marked as -changed, so that their problems will not be attributed erroneously to -authors of previous versions. - - Some devices are designed to deny users access to install or run -modified versions of the software inside them, although the manufacturer -can do so. This is fundamentally incompatible with the aim of -protecting users' freedom to change the software. The systematic -pattern of such abuse occurs in the area of products for individuals to -use, which is precisely where it is most unacceptable. Therefore, we -have designed this version of the GPL to prohibit the practice for those -products. If such problems arise substantially in other domains, we -stand ready to extend this provision to those domains in future versions -of the GPL, as needed to protect the freedom of users. - - Finally, every program is threatened constantly by software patents. -States should not allow patents to restrict development and use of -software on general-purpose computers, but in those that do, we wish to -avoid the special danger that patents applied to a free program could -make it effectively proprietary. To prevent this, the GPL assures that -patents cannot be used to render the program non-free. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Use with the GNU Affero General Public License. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU Affero General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the special requirements of the GNU Affero General Public License, -section 13, concerning interaction through a network will apply to the -combination as such. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - {one line to give the program's name and a brief idea of what it does.} - Copyright (C) {year} {name of author} - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . - -Also add information on how to contact you by electronic and paper mail. - - If the program does terminal interaction, make it output a short -notice like this when it starts in an interactive mode: - - {project} Copyright (C) {year} {fullname} - This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, your program's commands -might be different; for a GUI interface, you would use an "about box". - - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU GPL, see -. - - The GNU General Public License does not permit incorporating your program -into proprietary programs. If your program is a subroutine library, you -may consider it more useful to permit linking proprietary applications with -the library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. But first, please read -. diff --git a/README.md b/README.md deleted file mode 100644 index d5700c1..0000000 --- a/README.md +++ /dev/null @@ -1 +0,0 @@ -# cvmanova \ No newline at end of file diff --git a/README.pdf b/README.pdf new file mode 100644 index 0000000000000000000000000000000000000000..74128ce6e5742deb8937d69533603147d747e41e GIT binary patch literal 32601 zcmbTdV~{TGwk+D#Y}>YN+jh^kZDY1=+cswNZQHhO+kL(rac{(0XRp|?@BLFzzp5gh z9GN3W)+jPXQE@sZdNvrc!`rJ@7$#;yMnZcdD;Qp07zS}m8y8b225}oh7gJGFV|x=* z7zSBWJ98HcLMARQE*O4(7-ttJQ$t%AkM$ZYKcYq_Q62dRcHmGZ|&Zc+rQsmB=;3YG86F5gc_j$zwjL23d3Jw)D?m-9s6PJT51E!9C1 zg94|~_{hUZ{G}D0oWs$B5`Wx1ba6N;h+iD1I+{K6hsq)w+C`pykDNeHjQc*h^@$0| z^&mYXWV+olZ4>C$1c0d}9nto!*UA?w3ZZ4g-qbU$YSr9x{IdX$jyK6{FZ~Uk=HBL4 z{_=)ajHKxuP@9Y{V=DMUkc2VQ+ypG{$jxUidY9H{6AgAyMiE8~QVboZXTI&$Cvrzj ziT{ox>o(Sm|0;S;!|r8Xhi|IEn}<18RHP!?qG25-KPKGgbtio!`Rj@FHcWw`$L~@r zAa=+(bGS(MSaBT=mE6FQr6qVD{p01S>lm?W*nNI)F)6wH3cb=qbd|rwR$U;5PjsM1 zx!610fXnqgewBVC529|L`?{gD=}!Z1?q6NIrmt|Jkfl-$eRH{u(iL*I5hOJAC53Sy zA;R^b)QpPL+`~BLXLb%Rf6dI21#`#3Qk(LQg~O0F19CvpK*c z1T86;iW>?gja2LFIC4q_=YA#?cg^OqGto@s;cyazLO{9qm$PH9ab=W~>5vhXdLP_S zP21>-(5W;_ZH4>xCd*8ODFs%K4OaYlWJ-gp20&usI42HdE|3m=3?;a;sRF z%cO|uJ`~K14zUHZ8XHDiAEWJHh><1(ava=7rT=L2Z)VqsbUb)Ys}MQzK+<^lqx~{^ zA@D_PN6zPa)UO2cnC~u|3uJO3lX`=f8bew){)y)pP1<>cu=E_l;Z_SAcl3ZQ+*xuc z8W&a|LvZDF$}PE?M7Q8{9>JD{v@)$J#wE>V_bPk8CXyY-C&M zjb&MfS}-~ZN3kiWBE!I{sR28Y*J$x5Zjk~oQ#4F&^)7ovU3#nAeh`=4mDVLpQX&^? zUtLu7Cc5%WPv%2|$8zAm#|Vt5VfVVTjw_U@0~7x4SlsbmVTt3=KWj04JVl+AORn|Xq)z4AEnm!atF4<&wiCEar#;5e*0UHznb#RbShU6L=m_*vvBFh{IPrn83*bSxCx zz%;>(LGn^O;KVjQRkgIUG%ou{epDoDPkUh~(MZ|Swi;hIPm(MlbUoZ~5SaUgA=X=% zF49{Fbz}l616)lMLbh%kp{LCh<*BG>Puw}cZQFdk`CA?`(!~?YcJPQdTr77JRl30h z?RoQAor3_|mi=JTn_E|mU9Hwt7#8T7i7k-veSMlcX&RvYPAh(=Dq*-{904jpi?}~( zTOy{Mr$%51CH9Q-e6czw!E&oZxg(5rv+1li(vJ0+Mx)dezVz)|uWWhrTZf8tI%eJ* z_k-cSK)w0=PbP=Bv_Hl7M=Tg(PG)m+%^Zv+1C}-|Vqm*g+9_C9Q3hy@1DV^5Tl#sG zJ6)TucM)KPCf>1NizsbROx}op)WQGgQ99 zC2WwB&_ERgFr53zB5j*hqaJIck2j#LBbjIYN4o=V9X?W5X&i8_3kO_7^%W9w}< z9}2n^uJUJQ<8Xv6f*fm;$Ij8Elt_y{=FE)Ei7*ALtyLNGMuP)fbg?TUck2S4v=TP# zu95*3t%=G z;^*U?zcyNQy=4Q$>+F9WwA{pa9*=pN9{)JBnIAFWU3*rl1MA^mVEtJ*_@YTOB0{jb zFKF*_yDBG4)3gR}*RzWGz-VdR4wnjS-kO73#dOxk;ofo8&StgNuxzt1IKDH&2C5MN zXd5k70g=!Cu`+#Wd*rKr%}FSnN0=UW^>J67xDxgu;Ck|Nv9D%makG6%i{P-@+q9YS zbMf;yy5@cA*+P?bom~D^sbKmYl9@K3K}H?3u&~0>qRS&l?z@$3gcX4Xeq0Q!$v!!a z-)E$O&avvtt#r^{?9C&*nEOYoL&Pm8$s42+pR|mrT2^E=!Di<;#*cZeZuXG+ZAl+3 zlV&3(s&+>o$0BG$UlmqrXb5=`uGpHLfJ`CV*Q()js7q^Ti`0hn8FGaGN#d%QW1_1k zmQv-tYGQ&QShge2lHX2uixGQu$yVou1g2{#Hz(#*oz1A-2M6cB`vPsQI9veI^cM}e4}7WqgB{2PH0%eKgPCB)o&g&+ zmy^NIPj7V;Ej}Ml6$R8UxH_W-&8MZl?#~w^dha)L`Y}D_$&Cbtm!IU=o|K+H3O_$@ z6GtmMmqY1v86I!`oj224wbR_E5dbz* zbsq+jP;p1Gjw6nA%}^N@QZ0j<5^xiOzm%C-cO^SMG(B5W!5r5kbk)(AYpBQK64j{Z zk_haw2vz`HGbXX^6(7LY#l+5pXApOttDSTAXke_4Crwy50w?<|6)LhbHDK5sJ#t6| z(Uop0y=B!xfE;?$p4?rS&Ks=QnNFbfl0id-et-oWsaYnee2RY@o9ipNx}&dyk*~2GGLsPd8T^rxG`h~B`Rn1Ez(gv7|sU;sbTb7aFWj-0uK(e z7|Cj@memyFTdA8ziak7YfOTeJ*MblP$N!;wY0kM4kY>v#XmDlVRx=~5tss>OEg73~ zkqM9~k0KDxM(w$Swfe|5pD_CBAcoMfG9uP&83eV*@Cf*DcF~~K1lVs!8Z@n^%C79{ z?%j}5Yf4mE*qA&uv3FLF3LnX8d22HY*U#&tUj{V(l})*>#y3DNLjx7Mb}xj&o<{XR za6bw0QrXl@uNVxS(6^|@1Suq7yolv+_5M7`!L2+O~g&VZ11E^2)W<^mmD$V z%xBdroAtI{wU>4pK;;k52zua)2+WCMC~QsyBXc7x#LHh!9{I->Yoz9#6fWCBca5Gc zs6!yWk(peNNu&tTFI*IG!Rn3JGzQ#^FT2IDl~Univ6R}Q$V0y z@6c@u6as}Ntb%FkPq5|blX&yrZ_T~*AE1%+b!5A*_3Mg_2yUr5o4oGx_*yDXds2K$ zvCg3q)@l;9?cMnF1f#$E_kXXb=Df(X9pVjL;>WBOW#n;9obm$L3Z0G#6PB@EP71C+ zCude4xK6@k$i<{uc+OcI8qKq7lV4R~_!R>5$FDluGAwJ6T^%qm{m4+An5u%p7e^tm zUpegtl%yvGb_+47#cVI5SvU`HK{$<#Yl;baK5ychnx?bhmNX_TBbrnUQ2x5rzVJ@s zH<1iacF}wm`YyvqAcAiaAW;q=GjQ5KHo`Z23_tjQJR@lwOu7Yse<#Sv^Q?dM#eAe- zL%?Y1Bq}n|1iPlk$hJfma!X1fRU0S~c4?fu{`gl1ZynIhok}1VH8JW$Rllyb^0^fh ztipJ*b`f}EUP$ypiGX>VHsY-p40N&ZtJSqb(b=xc9-QXPq82*Akf6^E-7OlpAd-|) zQHxrCb!AY5aSA_?R^eVUtkU5H0b{FNHkKO%VvQU#$NfZ%AjtLFhSEKz`5N|`dN%F! zg-DcM@9y8BKkBRSSL^2xF8T{L{}*Du7=G64I8*xEiu!C| zv6)U&_4U(#)7gv2QR31{Z1(=B8k!jP4~y!4=Y%5oiK}w46+Wq#?rq}b=<3~Y%$&rC zBP=j}@*=-Nk}9Q~OSjTCME|c!Rf(bt)-q;6x++i>+yydrA&Lu~ln&%NHQaNL@=Dz2 z_+?giqBH-0rX=|w|Gr8$SGYSOp>wJ4CLrimyl1J-n+4@kCl<_TAiJyAT*255N+)Y^NkZ}|I4$S zQfky@5)359u&B?f>%329q2L6@cH!EaGSiTF^u_V@k47WH1?}EZo=kG42rDH4zE)?E z1?}VoFVmWhM31x{T9knO;12Uh7xZ)FsHc8wnXm#ssImxbib7P9-YS_(IH(5z%q7SU zRJM4L{clw(7rc(Jm1E>#1O9A4aH`BsHo-g!*C_~E`YDR73VoSjhY3*!Ch3h9#Zejf zcR19guA7Z3-ziOI02-mNaWo9=Tj?C+nLAO;m4~o%rz@F}T8Xt(fkk!Kgr3pWDNpue zWxeC<>Xv;U|J*7vXU4PDiW*+LlhgW=U8gMSmRGPk^yWOlN@;b@WNx5x^zSYb@-!36fmDwp&k9LzYzbXtR?${l2EC2B_~9 zsH50@a&CM=^ojC=fsDTX7U7RKdyy0qf|6m2^(F_d;vT-3-Z7>OsrW|qs`IkRX6VpC zolK*Pk|t-)=|tCZ4M$9lP$6USFS4Msk=4X_pwAgeQTpWKTb6NqEM;h*JT5susoR;s)t7b2yBQPG-PR=*m3iHg;f4^$*>*yOZ-wDd* z-@)acFJ8Xj?Cns)PQ)35 zLh$`U(rmI@8X?hMg~6$^_8dK2f17hy0}u^rVWW2hnDD=1Em&4T`Ps2ztuGwU2C_B( z;NAbKO7~tp$zzo#Fz2;yU_6ED_PB(t6VfSPcZH&uIW7L2KW2>4ChQ~7_G7N$8U8Mk zI2|$T`$qg(SdZ{C)w|}XNGqEbLK|RolSGARV$u3gh(0*huuMP3M0fo0&a>pL>YnGI zv8y-6i%zetUfzCX3)02fuFWgY-p6`|7V;WO0%{Bw%B=1IVqEJ&f5`g#hW_|bSB0Q< z@z9U6q~qd4u(9F`4m{g`cSf8Rm8cL}Hy2{dLEYSx;3Q?8B|}%rf1~TtZlG)s&$c|u zk*ABdfJ-ned%tIe9ThJgNiWvfd{`i3rQAPo3h|OF4Bg}-oPWLa(s?pw{8ilS$CAbW zFlB=qSmxTLFI?%62$hAobRv=no0aaziWO=f8)30D!^vN#s6a-;5ph$`x^_CIG*XJ& z%U0OzusAf@y(@ldbeM<9N%_O!%V}Q|viPo@NnMjs$@=p;@B<7I%i8!KEHdl=<6h zurdFeMNZM=u|H%-yndijJ%+h)d_@#9Xc!`D2RZ|3`=z|vm`yGA{Rd?_iKOQFa@XKe znkW0k=i#pO4@J8052Q+rO3sF<<-)tdpF5qOkB{qzjfcl4?`$A?wBxfk-|?&YQlBwb zFTalmY+*T~!*7nB&0NA446e;wo4Y7k>PR& zDhHAIb1atnrk_#7NIM^Y{Twqd%?dN!7`Gh`Zi|TaP@makqfe^;ru+`E4&pEPld1BI zQUWy77hdtl;nT%0?W$l|Sl8Y^lnIJvKT0xDhAKLKcWRM5Wr;)8LK=6ePE@UT=|rIh zL{d3*SjrQ#yCokDV|rYl#ydGFRWA&-Iblw2f57O(DK9-?z;`5^r5N!fC{Gf2Sg3uq z7*#9o4yCo)k7wewGe&*ycv51MG7v|?mAlC`g+4V<5?U0k#MYH+S6>$X>rery za4~wC92#OCG9_3Li&)Dcq$qg75{5lST&n7O-Wwbi>%NwByFj9p1P^#8kRHRxaW|r} z7{p(+JCa*q(88MCF&EgEJ9bG*W!kb_DpW#3*5uH0wK7^}ZQ>bUa*c?TO>JLD?*0wa zT$wIRNXAt?+jX^$|NPW@(hzT+$>QYVCDgkIR-=bZh1jN3JNkwz9m=x)DY zb=~hy({8etw`b%0AifChKGg6z{P-wd&{e-3Zj*tP4J^ti)+~GHaNN%=I0B+3OQoa= z#-{3y9xEtb9>P$8L;PLq7daaYEy`&?$AiV>n+Be-%Jy%iz`h=XfqS+-o`@HvEr0Zif|;O0JR!p!8gQrH)3|)S+Mtx|8cfCtH~NJ*l&&=Bx3~)zJQVL<^F_Mu`08 zvdIF`V5SCU0t)DKrfc@Rokv>R)O8@bi)4ZeHT2L7JOu zXCQi#4q%sMDR1TcRU}S{$7(?{)0t35AZf=%`@^8$R^n7eC)E$MzfyAN0s121;%5+! zZI{aq`SZhpS#rFtBVJKK3!#uoQ&}^4#5Dk_T*nm*w`hCE5M9};a|0@CV0El(eJlpZzr<1iCg5nVQ{n!D9c`&E1XX?z!t^s!BTo3`xCgnf6@{VR7nMeuzM0?7o3Fq_)|*uMUsS0?)YXmdp7y3@{`obuR5YveZO8M?f>ORuyU0dQ>%b8X9c_NXfIC_x;80N@}B3@b3y33i4ELrVd%rzv97( zbJD8E!pPJ~5G5G%e;-Yq%N;Mg0?MoE>@DUV{|Y7Yk=xi=Sp`~P`%%qCQb|=R(x6i` z6rNF4XA*HnGo|xkB8&FTI#K$#?DVl~?zl06*lG@iRtvUqRK{bv#AB*c)dmHd`QJN( zGReiDBw0|y39*1y8{lAoNx8SEC@45G=N_~dA0+#r0PYcIQ!@o7{|<-8A2kPk@{ zv;IK-ywC^-?i_&N?M_C$+rvk^S-KFZ22dv;q;N0n3g?^_{I!cnIis0;hne?Ch(nio z9lf;c>9P##TI|L&c@CUg7j`fBhoBOOq!NRL|oB zTRi&f-I?giTPwbgf4eHJIzNw#63UX?bDMer6{9VL8)O_CrY$Q?#S8)f4t>3*dJ~T; z_nW}}{{6!K=HK*bUFT2eb@1Ew^v7#Q0zd&Hc69%NVO;+|8OFu(Zw@F&c}D?N12uQe zbgU4kQao`@bJwCWK`?@RS(0E~IRXbKu*$y&xoBW=>u7kb{b2J~TSE)_x+}n43q3w0 zf8rdQ5y^)U+tX5+;GDp~$13x6mB8}dEW6LP(7isg`&53K`zG5vmR3zos*yI# z+bh-`k7sey`%qJr4OVTW_~9T(mW}uRmwVKqA0^<4Xw3GuJ;qO?!s^xf+jLapt@{2{ z>|D${UoX+erYP(X6yXjcU|t>g`<|R?H_3#GFm595Xy?crpyTe_r~<{t#?g97p*}x+ z=LZ4USYOX)fc{~c4WNw8su#fw8YNT*$*4x_Kmu8gBMgND`#=y=z`cy@otD&K5RH$p zS@J&>{3qr)9mHm)vI@j(==1Oej!EY>OHcocwoX+pNLOr#6g+J`dLNEnAD`2k41>cM z(pHlRd)5nOpuz&-qM-L`JZNIDCVC6q7p|%#jJ0HA!-N<)7YkZvZSqC^;6jDY8lUb_ zfnwWxaZER&YNKp$%KPB*j(D1`9h!;SL7ZZGu!Iaj$$nscPG?ZmW_1BkSmoKJuXTa2 zeWu+^OVW-lqpshzG!jirvfRNQ^uvKB)R$EfAb;BeVSU3xCiUzp*`gV7rQCY>vfM{8 zHAb*hJ6V#Zv7C08@-FZxbJVfeF5n#C^knpet3Sy*n#|nP*_(Hp2O}`U1MvIUusohF z`^+0WsD&&`eMZ+|PGPY;0-PJ>9`iU^7V>`HnUym3CmH(P-LI({&GcxI`R>h!--*m( zX+tkg-TSx^d3*m&rn}F7SX!8&`M~xpfv}CjJ-vEMyxc2$8}l)VR4}J-8ZxNC!Kde% z$-!Q$aXRYQ-1^)<{E~T_SVQG#6FL~7)13qK;C0Rlc*D0>E%eQg zKCZw8;li~Lkco^B&h9_y3<-b3ASlSe{#h7B$g6m?Oqo*Ip=o-nYIpK8Q_XM}9x6Ku zo9GF?p=J2p(`WvsTIs-~$)oBtlbWqkavtjX-T$0ckO{c()9p#dF!!NFb#ltm2mU#= zsrm1OWBu0{{{KTuV_{+ZH^Ld47#JEFn3x>GIm11e3My!A##f%BY)WpQ%lCI*0gN=TIg^kXB z_(dJXyuc!j%9_ZA21`db<%${&3rmZJ>#P#RoP_d#>VqP~tn7>;!#K^?l*!Fe0`@!`BI6S`b%EChUSv37V`}}G- zt@hk{IMHxf8UOK|QrDT&)#I=8PW_D2=ot_-7dEjv{OL5Q3K*SCKB3Wl9q^of>9j<_ zX#A4TobHai0D0-B;^oCSez+ZXY)GFu+57xf1OuYobJ6_|2(kT(B>oQz6wF-z{+&cu zX~JrN5wY`w=6O|EZ9{r@Q6wSTO*sX+0xnYwE&m zmpA zTB!E-hd$Y|Qba{o`WLUM%6Q9C^eTBFz*V#=_L`G;hanN`lU8*A1q8Tg? zJwi0ocPglRQVVukZfnSPuziZo4V&Rt)9`e)#~T2h4Tj^?c9rp?e7c7NM%YtjPPmJ;u=bJ1nmlC?U1 zxxktTN*l`uFi$e!F_sCf4lx>HbV2^J^g$Fbt(gWsk9?7idhCR_2#2!xgfmJ!@4t;W z_E(6sZX-xDKvxG5U}j{lIJ%0Z5{Vm2_}mWMpg{gQ8Sf2|cMKFb1Di6syTI7MqHbU= zTW;0U7K%!>n#wePe1+j7-lYj9Nib|ajTKtV zbsgvRNc?sn>rTuNs`L_T*3=d|m_c6@J1_fae3H@Mm8FeCHl<^|vk@VR%+M5da8$)(wpK1iit~18aFMG<4X17G=>7ZR%{$lp?h*uug&kIwH$b~rdtKASC z_Swv{YbrY{ha_aV>m3U$$P1#k><=cW0`2IB8dR&HwE{WUjm* zOe2;w_sX>0NU(kiDcm|2S~hamjN>NW`f1HgTXnqDdSE}2E zk(Jcx41CrEO`$pk8c)y-j4v>v1~t9a>eqtEBkMG&ju16Pmy{}JEYOR9r!KBZpM8V* zQk6`8z(zi*4unPXs%gXT&~fBL6A$QlV@KWC{E7B$Q*dVl*<{ZDl7 z@CC~MXV>FDK+gHUQyNybf48^)KyG#LUyy&Hc3cW~rT!spMue?d5Tsw%WxoRDeCGTw zb+@(H?rw}d9SSAA)naj;j+W4J-XwTxtG}#c*@{4CkM3}1&02CDoOZ!-O~j}s0FWua z-jptrYdx&5x4IFhmRwxUzcslg(?7Doof$kbmklTC*`+E2a_fxK^eBFEqS?6=@N}20 zFZ{Ah+fK5}2yViDH|9bt`o36JhyNGo55fyd2%%SQ(8g;5_QrafFb}hnK}ZGF z@{wUot1ASVZip0z3j9Z@)jh%bnuER2e1}IhS&+$e(9q4vcllLLD?{J+N#k52b|Nl0lp#4`Z07)Kdk+ zWn?c9yoz!{p-~!7&GI&`4J+H$laqUwDrY~T0kmILJ+^1%XD8^|#;w)29}te`*VO-^ z82)wk`5*Yt!uo$T`oxWY8a-_2%>!!prh#7eV?3I0elgvn5oGE!Fs{THY>X`~KK>gqDrapRJyt6{jKil- zPU9`bLHH`aKjBFAUwtgDlY;K5z2wH3iAcYmBgxHg4lQ@i&~pdY%C35)+2YHrw+55U z*pDHKFhn(#)$TfAt&ux|#vFiCtAG>dtqstq6#Uen%dhjctfLHV8GMHN_;@_rhg-DJ zDC-i`6Vz#n<%-!vyhPTCmLz31fzQHSqW_L1$V+BPf{XOGfi!s|B{-2z>gA_i;akQ8 zW5KAZ(>m%DgIHqV)T?4^BSHM61JY}r;9}}XVyRu!s6BWIi{!#1cL-j$8l?U^Mr{Ad zE{jNsird?{{8tdf_&*0hOdK5lPUIYQI2)A}G{2c=fIQz#!fxj**sZUfXtsyiSb+%lTqalV$ z`>exCWn->^csb|A3)l99OR_4yRL)zYyIbHS_aW!O_vS1zaEd&U6U(P4h{My}H+!#a1NmCN&bg?|zN- zyavrvBd5tC-g#WhY`WdKCEgfGlt@D{*k~{LNJ^%a%y8{Vu=pro6y8-%5N;&s?r~g_ zwFaCAu-{o~=-8F^9-3*Re{$6PaI$h=&9;*qJG7LGsgGh8XE7Ik99IQQppOr87jwLa zxl=x@N@qzgFLuy}`Qj{h1v@V!2Y1_cj8M+!SCb{d{R;C(c*(<>yK=bGC&|7&?neVx zF2aFu6+b%I)0z?cZ1EGX&rLwTj+{HOPca-388DY9gEi%I4%nHcQmn8NT>*uEl@i75mvgXFf%GlRQlLUQ#UzGRuy%pGCLt{saG1OE>G@AqJ<^Dze%zFnt zjYXKRLF>e<6EU>w?&3@Fpn1$cn=bO@DzyeVU~I56SR%;Gp%~<0r0FBIz;puO^57Z` zFS}v$Lge}cK#(r}a4mtx?_`u=yp*+JPAfZxIGt75COuzY&M7HX-4EdugY~$0FuYUw zuiLKW9(7+c%@Ta{U;S~Rb-L8Sl_3YCTWzTUmTZ%?U#!iCbG3)k3YlZQ%Dqn{y+Z{4 zE0u9LhDVRfsYN)Ew^Z3tbZsc0ib2)|7Z?Wv3~}BDC0(`>57WoQq@?s#;5Wj+$(_D6 zCl@Dt%!_b>-P;Ing%QuJT%D>w)IoOTfkU(4Tbz*ROqYG|DAN2TK>A+;#J{6c^l3)e zaxs$-;z3{0zp%@)HJ{)hR0suBIM?bhg!%a;0D!YU9IWd{m?=c<&5ZL6*Gp~l^5v}SY7fP zv`f6Ms5ZA4K@ssJdS=Z^Yw{O|bBFtC=&NSc)c2W)7xqfC_HVvMKQ)H2>|}lxz5e^j zg0IP?RH`|J#6MVmxq)l%!0ER06?GI4(I-;zh$AqXq0m-nZLkOBLR9<^fbx`nX}=N0aWFkWBB)ZjXK}RLT6c>I=&pF|Cg9 z>wUffNF-%Y*iSGxrTv-pcKZOy4q>LLgIU5#`ruX&&Uo3#cBd}@5j|CRM4E#?!}N=} zp&Ok+W_Iv4FSHuh97SWH+Y`q)CRdLG$;mP6W@$<6^nR+3 z*XMMYgWvMycgzp>orR&y$TCGHN1U!ydl$cD-|9_nn+8~S+z+rbp9AIydr;MFOhOaJ zHYbjX0^mqKa8z=%)SMRSDCcSPO-Tt%`5ON)2_LYV7LxI zs;cosu;hq@0gRzoEGUBUo}q|l85Bs(>ee|Ii+hnRE!x&6-PZ$Kn=ji^_52%UH_zQS z-V6l#{akAom(4$aYMOt#y6fW`jFHfz#Jt+BwomstH`M;RJG-=5gt_}m65simr%@6L z@8j8FHH79$D=e%hvPR6JMv)KC3gfK%T}sm121x7Tk_B|?FO?h^S7 zo_21TLb6MvPR|nR=O`u;CQ!<69n`0MYUL*uCWlFb2JJ6FDrn?5fzkyHG7%18{Y`7a zKOlyj&?R67uUPfOyZWb*okw(lY#a8{24FhB~Rw30{FBIWjcp(f#a^^T{{~ z$g$)A=Yl3?`*XVB({6-+Ha9ZJ1LYV~9Q704dMr5s^YO!Uaiv+k)aTAvPa84qf z_?l666Q2FROHnr>0FuTaoTJ1u;YYkkILCg_D>3mP0kkkYiIHevFkyitCb0-HA~MDFvN90Pv=wG0g8rfbtdQM`MCzm zCHaQy+E>uY=&$U-oaDigxAsqa3V?h!`_UCu*K_fgV;^h){>h$YFdY%uds(S`6`{x| z7;{jkWK7Xa^QfA*b`*OYX5Qni)D5I~D)UG*5y1c&Vkd3WTyaf~VfH>=TeEBI< zbFY2~gz@_^x-QGya;;LUlG)^rKLmENR)~k5pZOZ}ybBii)g^v-f!H?IxCt*4JKJwL z{QC8$$8)dnSoxTapumv)NCcKOe@F;LJEMN@+~&t*A;9}b?*%Gwp zu1#Mr?TyeV*I^bVM#8nxB370BRF#OgjB}G$c%9DI3065rp%N#X&q{z8CxfF}-3i5z z7=$r8H2JoIe)RqFc~zC`lBimVxUN{=3UJ&E%X_d5+8a`L09_7UiQF-AW#z8R_tkIs z7jLz&JWmhV4Z|ipL4u&*8(hErF7X{8rB7fnKzS%?Mqq61^yN$irxo?yBb%yg?qJXA z8YD%0X9PWuDm|FOX}3g~ds!gU_O=;9J?~!NOQQX>uCn2xeS0<^B0O zkLuOnRMZgM*5U`AoJQ?l9o6dsUQ7QSk&~|(Cv>gB`*P7bG6*6DZ?G!kYkJS<5s=9U zWwRGFbgV(f656|0Ajk8{wtop~D~X}MnJ!puDLVMYJYxChW+;0g4F7I*7tfc!go`HT zSqBcdfYI2R9#fc#IyqFadHsI%pn(-B&_*e=+?FVMX){}qq3x7Zkf^n*x)9dWmxP6e z@4^<&D=`5-?nR^3U9RiX5hxbg*WKJ|Es1|YBoAZ91>OY-Cf zxa=2?@TvY(Rw@O2)W-MJGQGj15`(4jiJGeXncTi#mHZ&G@1=X^Z0AK6x&NuP#DXhwqPX zgSi&IW>DbRjP;}bZYj82*s|>UBvsWv=cOF5DPMo;4{a z_8->cwA@Rc2auTfH)l4tUaUv(5Lhp|Q`YzU=T7maI42%Pq6Hb~lW1=AR{Fsv60FYF z=I0{RySDCX{)k9tAhcih?%A^3)x6n}mnJsrC2SIFYd1!1RAg>8c1BOmi+h{GQ?C;% z@QZQ+mY!uQ%0pF?-f1Bxjr6Zr5A{l1-a3j#Zo#2p@2DRQ(x*?wpUjd^)RJ^fbJ%=E zMVguarc(=KqwI$}qoBnMe2VitJbkqao~Lx>-9FFz^@jaa{Gl+(Bkg(5?>@lw`*W~% z=I~j_#=^c5-n}aJE;>XiF6!i=Mx4;^jPSj_#Rg=mjPZsfORC}ZkY?km-M#v>S%lNf z#{n2q0G<6ce}Hf20Xubk1M}kiOs;FKxV>HDy5EBGTGD-e!SViFN724DNzU>bz33FW z7`i9tGisW@T`$xWu-fndu2@!vtD5jpqxkZ2_8#%=7t@iq=J1R&8x;C@so=2tp|GZ~ z69iGwALoH&N8h7lPruxy@R%bdQ&B{{0DCs0G-KKc_C!u4iH)prDn=3;(3rlXK z6;wWc5kDN&76wc^^Qb51Xk|h+#<&RQ*AA12R;|<3>oF7h43gpFPyd73_yS4!)tVp? z0F^Drfvx1)z}^JJuxNH3HReK6nUmY?`QTENDyy(EV0kLMgqwn)l^Oc{ynOhAvXi|5 zi)7I>gw?qDjHWT70{KH8Rk!LYn|(tT&eBZgj+|Dv@|;sTombi8loWwtcGM(Lyp>+L z%xEGaL;K_(ur?5%?hfA*(xqRjlu>?>zAQpUE&@CR0ssCn?8e;JMvuxK_j)|8q}Y*q zm56y$K;EjAz{yt5BL)j=u?ULl0?y5k+7y8K8R!(CG{;0+Z-eC0E0SdhJ$oCcT3|&$ z5QDIHP_^*EXz0L!R-;THm%7))n8Cd#NsX0qihLs~h+qME89?9w3WfRWx4c^k8E)CM zxZ0Rln>QYN{94wZ$m%e5@3>k*1Nz@^=6$zXa=-Xj#cjKvs}Ob(SL%-OAXTX9k*gNY zoSDK`gQPW_9PCVl{L8L<^GzSirjyC$d?f)d>Jb9RT|NQCIcfz+i9-g&k@M+xp zW$_asNBVp4cuMqnGD}oOwp_4mufm4l|;|;s+q?wRaMcL2*5B8v&kFQBoCO)*a`t526C@KzwE_> zAFcqjG$Pt0JwUoeAy~dF-+v&F@xytIX5f)6WSt#^a?LhG+Z7Jnl!4_yZza?6`b|me zVF*#uEby_?l#72Ecj6Nf(hd{z^Kg`#H854$G3oA~}Tmodm4|66;+CY;ee>qJ%Lw%>C-` zHa%;6eDk4oh~IB30e)OYzcI(Fb8M|Nc;C0Lvu;k{@9_P{cdk};1?mr@yAr(jw-&(W z(W~i5mLdX?7^-A=Kbjv!{AvnmMGF)alb9>t}q_+CU}6!cOSRQezhd~dX-{%Xls6C6jYam>j#Vp zN1hL&CXLlWu#z1QCBp#yK>lrSC2vh&?1{BqbLyQ?X!Hx?^S%9(e{)vJovzm08Z&?L z1PhN0Iy4)u?dtPT4@c5X&U@TE>)xl`ot{(Au*Sf~Os5-P{yQ0F1O?csN}KwP1viAT z@+*Vfn|dDy+XW7c;qHK^G7frsJ(jd4fQnHo2kVJTfm8iFbVNY+i$K#mVf_8z{X2?{ zzU$4>f9-LWfUo84BT!&?b!50G`Ug+_hYUF;7*WQ>3*sa&K}b^+q0p*9l*G@A1@pp- zoMo$C?QBeT*Hf#^hOgF|Y<*o^!A);@gy0LajE_lC;2vCiyp}t~miht$o^_@Bf}eZ; zO%)MbD~5;1C1()2^f=Gp^GLHu-EsRGK!j-Xmc^Q~BEW0o_|eN!=~I8>=EumvNXo^Y z`MU|Li}bm|k?^x*lI*A+ylZ5TZT!QA^)DC0uPb?o;+AElxURejLp_0{p%eM$OI?xv zsXu0EEu~*7K|Ax?zpa((emfjcy&HR#b!`xfx zYeS;ACgeZ!daU+RcOgWkiyRBR>+)+he^OW=P)zyZ?IQx34^M`dHv^6C*Vjwi8AmgmdKeS~@IOuk2;Dfe(vw3tn8`I7iNU+n1QJcUC2Q#iiM z&~BgN{$=vW)r@zAnlc3a)(V23TBMjLPq462Fe=Q7+_W9~eze@O+&|10_f6dAH*&k& z_F|kxwV~6>KlRRy?-t>>@6meVKTH6)Rk4dmdHKtsAZ<0LM3dZ0Gj)UKEaU`nqAJp!k<+99v66GhU@h^rr6nL0?P~Is$zTkOC}TmJm13wJb)FUAOrtVKV}m2D}|ZEF>-$gQeaoMEVXyCYiE zUcMiRy2F9Js}JViYWat*{77#H?qk|L#YpXyWo*62{;3{kd6o>gIllXw3)E`o7u#N{ z)9IaC!fJ*~8TmsAsn(l4_c%RFt~R5w|47Cw8Pt|uO{zU$t4BCyOPSMD@3Qv7;Q#2Sx1 zR`+J`!A?UF-qh<^-;!UG>+~(jy9c1<-IUzdVQoQd?IDOA0DZ!o&-Bc$`OtI92{+pc z619eqSxyrOHD~8X*2~qJAL!ID?0i7;9sNMbk>2#2;Cj7q?xy)_YvETD{K?rfYsQ{| z+2hC8>36sF#7!6m6>gF`*?{m0slIWmUE4D4Ju{f}p$EH2duQ~Fr zy&3N7q|wzDM`PslfT}a7#C1s(#wPxRPOqt@cwjjenzqFuEOA2j;rXmDW&}t4rsjUqh9q-9bkj$P=K6nSi^iM%^y~V4bi%14J=J$@qu?A2TxNEO} zui}m>~+YS>d10uVX>pcxBlE{mWzyxWq5(jJ*(<`nB=OZA6DcbPlCMQG8ZMhx7c#LJ0TK;j64RkGrk#R>gj%mDzN(c zNBb&gLVY}jA z8Nhhm=r0zFb6JPudCeRp)}^=tKKkSn6&FCo)TNo-aZ7qQ2S3SqxdZ>R6ULUwTljoK zF1#l)(t5ybZgd?ok~lOh_)=}&mi7DT$qRR>+a@wByt8d$FQJ*Iz^MG;MhC0kqPe4*R(mprGWn5*bMc4oq!do@?HMPTJmz(2h?)tb%YA zTmsAhn~M4Loa`Q;7*V)c=El%%@h3F_e=>gvh`Dp%G$aqjx#?QF=InrNOks86U_D)64d4?kM0L< z7COPc#Ll>=S8>YVn4WN1H`^(hg^qqeJqzypav&P!G50VwvlNRFm<@&ySCQtity|M@ z)u1(;%ewugW z|yT#5c{SSYNWSgii-4DU}icP4r29%Q_D|5 zTcQZW3}C31(5a?UI~{=bq1BjG+=;Pps-&lg3m0Z6kq+?Vib)~spM7`M?L@ovFEvr;mA=&T&7@&y+Wv37P9gef-udrL}^hMbcQ9nYN1#wzti!dZV* z3@H!O@ch~ofzD8G;?H~SzdTDee?+`lS_dIn+hgM0MiK;vw}eBIPCb8=rhL+s70Lh} zo#Jx^rJbZ;+ac+baXs2sZ)brx>9=}Lhf?Ne$ipVU<6A1jS40h+nBtSJthGp69jd?u zGKn%qnKybBxcFo)WQ>N1O4PJO$4K z0;F}r=eO`dE@8?rymXqD+G93m=J9RPE0IxHpX-v|6@Atb=(8l(7?W_TJ@}0Slu#er1f7crp zd1g|}Gj*SsoWJ|#7askH1*WBqoS}>Ea8kzP#U5+PW8RNF^R@vkl>`i|x*huTYc3YN z=GXn&4eCA*fhFU@lbCJz#v2`rY5c)!OUIdKKi;F*^9Hjef!=YcV zc_JFka|895AU~Osr>cG2Mg#a6r~f|iIIQ+I+0p8^pf;BKP=kg_=%Cu|bog{U^ohKS z?A;r9-K{#*vz%`J^fV%pU7N`WB^U#noMkfS7CpZ$Q8XU|5#G+ zwXhk@)WgKGTR3NOo}wmJ+qJRMs$D&9IeSd!>SmmD1#06DpE>y<-yIxQ`wy+4|4XEf z`Ct8$R{HvT@o~Fv52s(AVRkBNl^B0f3o!3K5P%s92>t)?n|~&E{@*)Je}SHVCHg|t z|8|_zp8XK-10%}1Oh?A8ZX&CL6j7?6Eww6zl7S8a^|u8yh)4vG8W057AO-hAaicG} z8-{$Szbcjut!)nTRxGQzE|j$@Wvs0J_{%0j`Y*RS{$uinQYglOMwRAJ^PaM&0oO*~)YK8rq+o^L%yTOyZbm!8KNp64X8(I7psuWJmz&e!SM@Wk|m9fT>Wh}oyFf)`#}cBALmM6|wMC_BjT$Do65ebXE=3!oMA;m8DVz-tD>9aD%lu`XS7{Ma(mT@IWU<|Ob}YK?m1l2&2`qMb_F zP(6}t63efvuc2FzqJCU@q(un%1C@cuTS8Mrr{5k!*h54EvFc(V2sa3VA`JCe z5QM78R=*jDvec!kh}n>}ep3)-#33;e1xOH~1TTaGr*M}bC|M(LBX|Cu>lXlq`e^(>T}zFzAO)o!k`e}PNLmHKqv?*K?db$TVQwsXU& z*774PK%2svL%}1_+A<&qM*k{V7C<>)?$1;@Y8MFlSV7AxKEt;jIl`DF_CPDoOUqN# z#8rQPOQ-_y?0pf=Gk)Q^X3=R@hfg^M3s4v)R zP>N!Y%#@g~y#rT?dGuCE{*2>rBc4HB;c1=NNNZgk;q*#$ z{DkLm-qcEeKuCY&0*<(w#n(ApnlYcEjEQBzTxDlGrw7_xrxBl$o7f7Ba`t2a>OuEW zNNibPPJ?3P4{XnX6=HD|BevW2!^#}Y_ue`D1x=k{MB9#LTkLDfPMA}t1NxMK`#$B3 zP_nJBu-NmRdYyYI0`zzqQ!-FMRZ=jgo*o(?_ z@H&=f-krqz$CexUw||@`GJ9OBo0hz{)X&^sD16BjUDNI>EZNOp=X!{&5jn&Za|*!S ze=1H;Lx`P)@eH!P6Z)B9ZL*yl(G>VIcuFb|l8L#>U=BdTDAY6)!(8wS?Gb1+o#0RG zYW5oGCFXv*D1*^Fww7tYtwo^BF#`kt-uF4Ny3%>xq4QGgNVNWG8}dq>N^Ho*CU! zW+iRnWw?6sj3(DTJg0#>24hoKrGO8N#>C&nk+3AAvlNR z!dEXr)eB@cWA|-~vzF{J_XOd0F~vNVNoaEX`Q5JjTx&MBit7h+T0*y~}Cw*ekSJ`z{YV1$<5+ zNce`1pv>B%6Y(9QTAzTC1C39Rr^eM8Vi}Vs97;%PCVPO&SO`0)3w`N87|TTg1h%?& zcu+W5U6(3-yroGTH<~-1c39{o$AuNx&5-RJQUD&8V6`o=sK}Osg@u73N9@0oBKrXs~s?2#K$7&a#rarD9TX@Wu0H*QX7AZssy{ zlFCq(ZcFjd#L+sg!aSaOl2H^l$V&3_$V!y#rt+GGo&?Iu_OPiD%T&Ru8*;wV(N44% zC7-e3XRm$H3*eR{B}HIlKq?d#;UgYy-42^dq$P|daQ21x@Oy1u;q5a@sQ?S-gzXcw zwdp-y!9drwZ!OnUCnk(0Z2=#9tDo{iFRDX!z$)}2e&F|eJFu@&77Pyi*LE}(M_0ki;?v)xF?9>@w`#7n|60Grr ztcs&|moX>%Z<6wIM7Bkpl7t`_t!0{5Spo+M*kSso9efZ5mwkY$m^T}jvmg9BX_xO_ zL>}ckcQk6t1NSyaQmFZ%^j#w$Wad(b@!c=2T5|awpWj=83L$gU#`V2^W6R(#{M>$b z2W9~~gNRs`g7y@5D2_--Yor4o-)@YM_65K_F@oKWB2&>WeT5?V=0M@vbiQ7xkW5g* z#JlN$WW=w~Z?|g;&1_$F6M?4d4Vkz&HQI89i8huoxi7P7pG0hSr;fJ(YC^t|6MH-S zBBEkcTJU2ZXL{HC0w~Cj|{$K0w~2dAb~nt7czY2QL}ZO>s%3Ph5DOk&TuP zW%h2m7AJE4{lJu-CIXLn;Ie`TfC}A>1*$&BiVYAnJ_@`@A~A*>GO?I-r4G z^tEfkttrbYjBpOT+QcmWGy_)Uv7g`qwa(+a3Fgs+e|u;|y?mL_Pwn06!EXhtNA6U< zCRk$=%co}Hj+2X88$EUEUIE9ZA8u}M;Q`-j`H46UVfQ-~y+%~sght?nT8>#r24v6{6f5iyfQ~m` zVQxRZjX-~3^5F8K*{*dhfWH#HnIs$_?&9gS8<70Or#_ikJ+WHyRT^#wj~4tzt&uaD z$bfK=6!hRo#o)vAY7f3eiMzCqZto}0W=U8nnW6VYiix-Ipr+XFbUR9X+smXjYlaBZMcJ+I`Fg;07-jx(LU;Jvfq zEh8c>63j7WmKiZ1S_H3hN)9Fp$?A>^p*_UShk74lJ7eG;6F^NgkIZDCp%_ri7N~qk6z$85RY0gT1L5S3Fhp`@b zM&9ALbyu0Qp~FtjyA!4k8AB>$0;!^f*sjh&?}w{`2_A75fW7aK{0W09B7>#5b$PUL zv&TIDN`HJ<+;#gjC#R}ts^SYHqP-rqt?CX?tkSAZqR;eCfsm)7o92Vy<6OJD3oD?G ztk%|a)NFGIud?y7@o{nyul6Me0q>+kS)KnPKVTbII_=FGRXBxewWZz?`~&sS$Rv7C zrVj2oKHBZOA5bZ5Dj&_GGmVmZCrvFyRVgC7iM2%8Q8mN~9~h7A*z0o>z*|^fymj(N zEDCjJuLW=k0NY7qjifF!p#bbJko*^+(zEe&2ImY%;5_W8miD!p;SrGiy`iQ6S0<9_ z3hQzFRyMPVCvjdgiAJ zdLXVGu=z9ELl9mBM*}e@)&u$C2E6pX+l|W>&xcRl^LeDNr9<4R^4sdTkK>2uAyaX! z5lxr3($Vn2nsl~h4#q3RHbiIY&#UE1$76A*d5s5?s~_EtiuN(rH9uIZ zgzCS`Kem511^f@j4W@r(|No&f^<%eu&G?gTA>0SFR7#s9 zYYlfOuqtmtMzmzl%V{+?b1hql|2SP4w_3Bn z;WRDuZyf0}cfqFXPWUd`?+!7@wIsfM{28B{Hvs#X^vD>baT?7Vw2m?hai|+#Q-;xJ zuGg#IGgF}q zAv~B->R7i{G%`uj&@ElNiy{NJP@M!5vq;_8cXry%Nxx3Ze|N(rK=992a=Da2lMG!H z5+$)ADb$M^9Qef10tDoUZ;6ppUHcl3HD3#W)9@C9n>FG|IPI{i;f>+W8kRNiNvzEf zP@U|W=_Rl=G+HpI5#gGuH8)yBh{3)(2wM16zab-KivQK{Le)C0DgQ(L+sSlGe*_0I zr%4tdHFl*)<)kWoHM#<&-gltgNHY)r4)Z2&i(k3|-@X&KerZF{v+i^H z^V72$P7$0NSqZ0M3}6gcsfAo+Rz^0SF#afzQ61tA3s)|CKMZpG#FFdn=5Z4hQE3AY zya42l>DcaE)$RHtmW7g7apV#t2j1JZ@3t7(KnjPT^dK;Zf6|ih%Kl@5zt! zY+p}IqyKFROD^>lDeSZ}c!=uHgxbRRFAIN;G&;EAhcVfB@z@uiRbVqUg&69n+pf2R zKIO$EJIp&a-%Ra9=($J7cl9s;y2EC)To_yMQq3ow3{|bO0jpy9q#N+=rLY(sN}Cz;0nO^(o{v(i|sW4p$ zfC;74huQci+N$v0Mj+k&R%#YbfqX$#6Vk61d(v3}zou$foNtCaqEMOC#e|_yB3w{G zNZJ_WJd-Zo_t!-!1Qfl7N(p~&56C8rX*$4@WzuTr?Uz{Wn+wY56h1oyo+I4Cr%}s* z?=?=@$YIcjcY3pi$Vu51nF&?CKZ;-i(0Vh$&<;M-+qSk!WAL$)v(csB&GF>(Cpe&U zv-g*~Hx;rgw~in`$6~tT{Q>HIaBAZg8EROi7+4E(#&h?1T(UV!L*Z@azj_E-0|^tj2ox+(K=u8XqKC`Ing z&Xp|H4CK8CNhtHA>5f{$Th}93BnQOUBm&>&#Kt%FR83#k{lXET*q;+)m`sAhz*b~f zj<2C09R$RsB&7JpXhqlegsHW-pZi&W_XO+CSu@3kEU*`74aMa(Hs)lw5nb$s=hAV5 z@{Re-MKRyn;+zZiOyJjnPF4Ghsy&Ps?rV6(EV@f%Hw)K`NmsxXvYku~5!Ox-i%?1+ zR6d?Mr}U|LsP;wCQ5$3JW`R#tix#&?_;fm8=y?_Zq6e7P`o&=uyvsaQa#&XV5;#2$ z_)z;y6K+%)w}6g9D}Bko`1=W4@nwZumLQV9t%J1-nR#HgePny)d?JYB$Cc}|q>7da zlI9Om#;QD6Yh&IaJPvFe;+ywL4%+BMKtz)YCE?M#&iGx5rU8BnFjn&chTd+dGDmOi z=e~-#DV#;P?QD8K^e6UX2d#bt1_E|6`CLs}*+w&b;HR9owIXS#rdcTdxGH` zcHzBb7Ef}8efRc!!}5N~c+7CxOjce*L*CMt>QXd*--BX=G-`^1{!y?Lmw298Sljlm z8?0b9#Bxovz44%juZ>(RRRjT_oEei%fHHKxub>e)lD~&6UNwVG;e`)Y+DkeG-US*E zUrTiL{XndvXD;wyY=@hjESkKhPywu-EclJNzTUB5s;Tg+q>=+bYQiho=O|XV6;C0! z)XmFQs{@c_xk(vkY3+I$Z-QT71-HVYa1P9PdKU?;g_h*tQdD58?DKId9XeI+Y6FTO z6oG$=+xuGkX3^}Tl8YN#(8datSA@@B+Y>W#;unVfw{2iQ%|OM*8?Ix@tbszvkI~1u z3bv4uk`f@y(H9)d;6+_Mns(qHxW*+jz2bpAuFv(W9eKJ7<4i!Q@B5s#vl)ixbu459V?JY8=^_<ED~; z%xg<(ab^a5Cu>UvtOv1HyV%o4RX%pu2G#R#i{H23z0ZT^m(d(U%5{yTi+BFN9=I~q z@en4bNTERX<3%rqv@L<5rjbs`P)KB`hn;5V0$52kH&bvptI2|BAY+L(NMxPKF;pK8 z5-lJ_m8G;ScrkR6UJZ?}MPO?^UtCbMx11DI zn{pts&l=?T6FH8o^Y$g|z8Ug2vFj<7Ap1-SMlKryKzl^94DlAnr~hKIS0=5DB|a>^ zAORKKvZ0sEUmSalm$fU1{blY?E{7L==`-n@d78K9i8ueUceI`u=6wDQQ0`G)JK%TP zEa=A&j1@gw=$p1>H+{sc7ae>?%OlKqq41}!3Pl`HhOw1eTBRpdn_rb)Gi>?3ekHO# zgT%}j3kM7a-N^`UB?XUc(`v$tQfza=%g2mHhJf##UYS59yUBO8rsY}m8FY(Qo`2|+>m{WmUrP(Ay<2S$cKP-3(dP#b>S9bVw-P7bp3(C+O` z4t?2H@r;UBoY+={KDObFW`gpocfOz$b>k9BD`skG{&LSl^Mxkgf$j6Zg^Iuxj*n0) zKs+{*R@gA$bo5c(i*?Af25Z4zc^3xnh}uPvwsY3jH#g?{H9gRCj{7R70DhF)6=&sg z5{ODC6fG4P*o6xvZHYHoWSg<1(+;LmA704L%pw?A8!nRFr&xFn+#00l#3iMr1$8Bz zJEFc{56m9Fw}c5QW!H;pT1a#Jb{6${C||F)y~g13Tn--qq@NiD)?$sgqdu_g%u7wO zQSP*yn6gm_Z%d0|?*uf>#yF&WSMELNLwNRUKeO}k-Tj*xuF>RpUhk$p@MBDX&P)SW z$lkREg2weP-=Z~2Rjqf^&=6o^?tj;cNrfz{m$XNMb&Zs(Zew)Cu}0DvuE?w9_JvnT z)Jj{BW@1~#Mp0>=LbHF$n#sm7h1((^i2-%Lt`3(rzy*#Wd%aZGyi~Y_mm#|3qIxg5 zU6VB}c1ku8ly`OHIVi7ns0UXUFz94}cMuC(3P~DYbz8oL#(wquQiy9!l|F+7-7TU) zwTAZdOSU!}j^`yfkeL6Q(y-aFu<-8A6e=f|9|uhzl+wy><&v!;;5m@ykR9C8dd=@% z9At-S&^@048I<_;Bl21DDpl7LnyFCVbJDhFZ#H^p=e#iF_&3;G`UazZg2hldo~svx zB_w3NoJ6zSTdV8L<2}~0g8rs9T?=%6-T|*Ou%b65=zXBzVhE8NO?|4YV)=0iz+rb7 z2^f->yO-1K*X*1fOi0x9je&Z&?|~oP6-qd=Nqs{{8ebeHk(A2uIvPsl12auR$J;om zxlu|UO)BTfo@c`)RKU!qX9p!&4GSl(HPv=PA<1q&eUl;sdd*&RwwdbIAokgbdSmk5 zRWY&07EJp`8HU?eqx;$O_9A@W?;1P|`ifHI3dUbwQk`9@yoL@1cJ*qipHv-{$bG`S zYMNzGZmg7{aydNuXMjUB?y|@_@-Y|I;!O)Q#c`AP66GDnBQs+b0o>R9xGOS@%vQ>` zG8n>@fZ{mL;0}AZH3nORwgVW=8xCW_v2FFN(ay;0_P?(yO#4Co4q9SM_TmRIy;%5X zt}TGPh(ZDy_p!&I)~A}YYJmI%?dnwePcTsg{Pwh_#UDQ7w=b5z=ip|7oMhzf;9+OP zodOEqVsj;7&KXW+VFtUV{Epf<&D5B!D)%ogA}n#%)EwEU?>0EH_)5cij-|)YGfb=Irw*Q0D{b%Fe|D>Z?{>7j7Upo4`hOkZF zPUiu<6$~7J>HzFN_yXDfxf%TbE~J@Q|F@J~7gaB1WHqc@1h9q-L(-B=+FoN6@DO3{djz%F1INdl5M@_od6B19NT386k70)#@8h(I3HTjh!swgu)Y zoU6qgn^w=#?C*NWtLv81&X-v3zhe63!?$<$g##ZpB09vJWn*C~>PVJbDa8yyLyLup zD4U*O$nZO&8mdQK8K*l=hEIxadsajeg@_R=d?=11`)KGQosqvs=7|rA*{=K#qUPJlbO2&-tHaBn!m)=HNIA5vb{l zB9+P&>m*f6W1;|sfLh(=3;9{Gm(*Lua+uhz(Lyl$0!=}QBcgJV7<7xy#9NqriI~db z21lOzWJJS2dU7fI;-cWVeJH{JyY+5r;9Z(zid*Gh4nwp7wdw(fj)8q4P@=~j$<(M) zf(goG>Ka(>#Cnw9=Ly~>kJ}N`c10#t65Ab&q&8|&2?d&rvxUJ5x1xp~3$;)FL{^lK zG><$@fE!T7LKbs!=Ilw(clm~8FbM)GaT$2?%nfM5k%bc@LV#ylni^9DX^QV?@2&Nl zM!f;D4vY>RZMBCl%y?#5#Y>Ro9acwylo?-K>P|s&d`_u0P#IX6nbT;tCLY=OaPyfu zd>_7<%b`9}bPJu6EOo6#8@k~V%ceG@DO7Z|aVE@As*&X(^bb(o1v7j4b$$rA*f@NQ z4*|i34m24h8Jal*IyYCCoS!@5XIQt0TH$JDDJG07rbQ*mgs{=LVRWH@Ut_qJG88Cj zLY`c-E=KADiyeta3zQrETF$q7&@J_qa`30tB$mNrbmO;IAHchEc!NOwK{dX{2jgYpeSUbPC40ZHANz&qz#JxZ z44(+WLfq_J;E9>XB}D26xeHFKvlmxAvR%1+A6%_cq94-L~CBySs09sj84U)Z*NdcRqvHL#96E z=#I$)mq+&d^F{;@DL*hRLb=Ljjvur=r>kCMYxb&zd94~^1#fMEVzbgVez36kl`8px z!^NM5h_Y?bCFC^Omz>s*%)hU$^oK<+0 zxScSb1}kPWx*rDowc@M!Cw62^3%`UYP)*3-1JF611*AK}&{*eph2K8fic`px}+q3zpzOy zw5>Cf!?(2T>2QDF+&{VX#9z}g=5hf>yQ;H$dSL~--uA3i&3t_xJFXpbs%_O)9b0b; zt@TtzC#|87G3vw1Osf{SgdJq(!wHsynLW~V)qU-IWMB7-HC&WSu`@#3lm1wf$lW+4 zS#BhOrNw4!{szl{!x^m-6EX5HW-OCLI~9Kg7bFr>xVgSK_$T6IjDFY%_FJ6qDn_WH zEKM`8IlUu$`x%d$9VRV+o7XNj_oRU>3LXc_m3W8Eshj*lhh2t`U548Mt;;?Hjt_AB zQt4m^7JN!Y%1&_Sf&e9xqu3(Qbu+QPA7-XvcftI(qR@?G3naqxWEKmuiO@B+xpG}DxFB>5Y;h__NYCdLh=3CV&9G>N zh?|*z%0o5$EL=KmY>(pPf635#g=&6yx41tpxQ%}M!jtH~159Nwlr?`8N&Jf6Rt;uV zhy8Fqd6&u-$m~Chd59hO8#_Rba%K%hAeg(1HbTKQeAN5Oa+7?Bi7p-nL_9i=Y?P?> zUc=N;`lS}?5nphP3wjUVm%54^T{>!OtZf7-NvqMn(oV`HTG$XyUBm&~K5*1vZJWZ- zDG*MfMaEr)@kf$=4Y>8O83y-6HVue%!&6BE?tOtS*7RL=?fT4~Gro?rvIVm;p9WEG zdYGZNUC9vvGpiGNFH@bZ8LvdK;hC(IU|>R697i>EIr{vHlc;UIsds*g8V5yg3!HJ! zpzez%10h=!pA7LtDZU)(-N|`?C$S9DI~OVM@ydbVGD;kM$W$V|cXL*oqeI`-Aa)w7 z-puMwYcn)$D$=!cd1Jw_c>||V zZx#0;6dXYS5J-T;{(X1Ge;xwykBBLHc|%JjXB#MbIRX~;f5eVEIynDFL|du=0j4K2E4mj)*2{J+xz8W7T0WGshKzKsMlOJL0+n6Csun1Mu?3PBTNFvCc5 zG-7o1(H^9jtpwAa0wD*pLM4W~>}4=dc#9hKc?XMN<_yvmc9=pCr}IOCAeF@~(o;$mU~aJlsHo`q>-QfB)3Bu+*c5MV9D1|jS~yJYH?evNrPnUgOoU6rSrfS3}WZ;H0%7&+3m2h z5AVu;qkIPy($6R8Q)@wGnMfluv+nlTIMb0G?x^UJWJeYr%CZ0aa+bqy==`iBxsh_r zSlGyjNY)o-OPUsdqsdetM1$ge^N@90&I@Vxi{+Ij_7&VZy(G*f^AbPQ{ZlbUCW_n_ zT-IIr61KhgF#BU3S>pVJxSrI*SCGMUZGD5-(4YTWL7pj9B#+oyWALIp^34n3qh8O zcf4pQ<-o^}gw;7!dbnPcDvy|aT<(&K_b>zM9Y{P$_Ewh15W^hCVenWQuHIfjCp@7x z8E2RY^WcfuL1r}0Sj@_NVD8KZ_rAE`Q00t3El*(MJahj|;KlwOk3c2kPtW_H+2o@h z!IO2Q!I*RDGtu#;_Pi_TFngVE!ZdETZnW|NpTwH)t^o7E-uiY$d)hJEUjNFm5$FVI zc{ds2p_LyUI;>ebAs;vzFlyL+wLhp?IoU&OIl*_A=$%mM;!@jAJK#;e__2f%e}hz>8%O_PR!MY z@sgx`@N0BY8_%0(Z?I}C^sp17G>PK`AE=KL>Sy;!9+4(*u-T_05PRD*dS>}jOOO`} zDSBo)Pj5?}ayj{BPoIev?xqn^(X#M?GSdyl&Uq)+^YB~W89Zga-<(VR5j)#$orQcq zJJHXFr;TT$x2LOpEOfAnAH&a5CkDUaeneE}g)B&1KZ(hVv@Ug}T6J%&VKmX(-c64I zn}mS*>X$aB#o;jmH|=?SkkYIB95==`tP>b>ZN%q3v>`B8*!p(Ty5bu*R1Vq#s{^I+ zFA$eGDq6Ex@Ijn6-*^}XA+t^&{y(91E|q$2mDL?85MNZiuAO@PyIHGOc{Xm7Pp`8H@ZF24JkqgFjGSE^D&IowM2B0KG`Maj zKmNpqv44b_r3Tt@$=)2W)TYQa3aVw1mA%kl*tnxbtqDX@ac|j4Y=*9W77u6iT59D% zi=K&4=1^s;D?5*rsh?{`7pxVUHc>4z#fz0-n7dwq0K1&ql11bW*)nC-QuJFXbJp$t%4cI6HIG zDV=TWYs%z2!@VEcC8f-ee|ZDV-DE%SPYK>Fk(k8B4t@q?VY_wt+M2@o7h-oP2S8Ib z`MD-M6h15{$Ij5sM;|Hh+M3yy>7C-r`t!f~q63aaRp`-f9CFYmSF~+XW|m;vQwS1N zqmTt2LYBY6n#9bzXY}qr62Gxu-)4TFziM^s*zy4xWG1SuF*;qwZ9V{n(sdsvMkOfr`vG z4VX~fod|Nf5IM#Zz!cAr08%`)h4S)lKm^5Ii3V<<#-d{%Yt|6wNgGf-0qUPv)P)Zx zii399Dfn8HD{ZkyUPiFlMCOFKrrU({k+lIOhYPQpRS9nPY2QMnVi*VlucsAY>6ltd z%Mb1`fMstLA(8EAOXrK-6)P4h6J_smYJYe2aB=aki?=`TTO=rb`A61;tfj#=j@w0I zbw^>HPnOl!MrSl%{&8mARt*eIsXh>Sx6H!GKc!Un1c8S4)yi6&1=m#f?bm#fnf4t;(Y~W?Ij{La|s!ir`wJ zLhHEIglbJduD68Tj62g6F?0tE*W%Ue*Do`ei@*XamQ*3ROs8Ki;X-Akl$M!5ifRBk z^v(f*JecqRI>_KrYe2wjHGeFY6|~bq44^#q*NzHTH0jbkzV7SdXnXbE)@h6XvDMdX z2bB5x0cN%aH?7j6>JtG2XZTe~$PbkEd;@->8e8LUp@N%c3OW@>;96Fr*&SnY|4*%dAgpR*&6wS34s61!}s@)Vcz#sHxg7Tec5|3eCH7HoKbOSg8vvCw+Tyc!j?B zXR@XmU~Cx>CFe=>%jf0cRK&`GpgM8lu$@cYvp%;!hVkpl=DW02&dr+_OW!A^bta~y zOa0fSZTxoA*tPk5%2Y!qt?wdQYay)`UE;P4qW2<3Zy_Un{sTIV=Sf)#$2Y#~4cqg- z*DUM|{}k-?Up0%mqKPRKy`r6+Gr>P#8AT|1Nn2Apf`810|GQSQ@H8R#o4Q4>rp5U8 z8GzaitKq@wa-Q2z@pl-6(n literal 0 HcmV?d00001 diff --git a/README.text b/README.text new file mode 100644 index 0000000..4cfda30 --- /dev/null +++ b/README.text @@ -0,0 +1,101 @@ +% MVPA by cross-validated MANOVA +% Carsten Allefeld +% v1, 2013-12-20 + +This text documents the Matlab implementation for SPM of the method introduced +by Carsten Allefeld and John-Dylan Haynes, "Searchlight-based multi-voxel +pattern analysis of fMRI by cross-validated MANOVA", NeuroImage, 89:345–357, +2014. + +# prerequisites + +The cross-validated MANOVA is based on a ("first-level") multivariate General +Linear Model. This model has to be specified and estimated in SPM before using +these functions. + +Estimation of the model is necessary in order to access SPM's estimates of +various fMRI data properties, especially the temporal correlation of the +residuals. The functions use the generated `SPM.mat` file and the data files +referenced therein, as well as the mask image (`mask.hdr` and `mask.img`). The + beta images generated during estimation can be deleted. + +The functions may work with other versions, but have been specifically +developed and tested with SPM8 and Matlab 7.11–7.14. + + +# interface + +The main interface is given by the function + +> cvManovaSearchlight(dirName, slRadius, Cs, permute) + +which computes the cross-validated MANOVA on a searchlight. `dirName` is the +name of the directory where the `SPM.mat` file is located, `slRadius` is the +radius of the searchlight in voxels, and `Cs` is a cell array whose elements +are contrast matrices. `permute` specifies whether permutation values should +be computed and defaults to `false`. + +Simple ("t-like") contrasts are specified as a *column vector*, complex +("F-like") contrasts as a matrix of several columns. Please note that this is +the transpose of the format used in the SPM user interface. + +The rows of a contrast matrix correspond to the model regressors for each +session *separately*, i.e. other than in SPM the contrast should not be +explicitly replicated for several sessions. Instead, the program performs the +replication internally, assuming that (at least the leading) regressors for +each session model the same effects. If there are fewer rows in a contrast +matrix than there are regressors for a session, the matrix is zero-padded. + +The searchlight radius $r$ is interpreted such that every voxel is included +for which the distance from the center voxel is *smaller than or equal* to the +radius. This means that $r = 0$ leads to a searchlight size of 1 voxel, +$r = 1$ to 7 voxels, $r = 2$ to 33 voxels, and so on. This definition may +differ from the one used in other implementations of MVPA algorithms and in +publications. Note that it is possible to use fractional values for $r$. + +The result of the analysis are estimates of a multivariate measure of effect +size, the pattern discriminability $D$, which is intended as a drop-in +replacement for the conventional measure of classification accuracy. +Statistical parametric maps of $D$ are written to images with filenames +of the form + +> spmD_C####_P####.nii + +enumerating all contrasts and permutations, in the same directory as the +`SPM.mat` file. Additionally, an image of the numbers of voxels contained in +each searchlight is written to `VPSL.nii`. + +To ease the specification of contrasts, the utility function `contrasts` can +be used to generate contrast matrices for all main effects and interactions of +a factorial design, in a form suitable for use with `cvManovaSearchlight`. + +For example, `Cs = contrasts([2 3])` results in + +> Cs = { [ 1 1 1 -1 -1 -1]' +> [ 1 -1 -0 1 -1 -0 ; -0 1 -1 -0 1 -1]' +> [ 1 -1 -0 -1 1 0 ; -0 1 -1 0 -1 1]' }; + +For further documentation, please refer to the help texts in the `m`-files. + + +# remarks + +— The functions are optimized for the computation of several contrasts (and +permutations) in one run. One call of `cvManovaSearchlight` with several +contrasts will take substantially less time than several calls for each +contrast separately. + +— The function reads the complete data set into memory. The analysis should +therefore be run on a computer with a sufficient amount of main memory, and +using other memory-intensive programs at the same time should be avoided. + +— The estimation of $D$ is based on the GLM residuals and therefore depends on +a properly specified model. That means that all effects that are known to +systematically occur should be included in the model. Because sub-effects +can be selected through the mechanism of constrasts, it is neither necessary +nor advisable to use different GLMs as the basis of different MVPA analyses. + +— Depending on the data set, it may be possible to perform the analysis for +searchlight radii of up to 5 or 6. However, a large radius leads to very long +computation times as well as decreased numerical precision. The recommended +radius is 3, resulting in a searchlight of 123 voxels. diff --git a/contrasts.m b/contrasts.m new file mode 100644 index 0000000..672d7aa --- /dev/null +++ b/contrasts.m @@ -0,0 +1,60 @@ +function [cMatrix, cName] = contrasts(fLevel, fName) + +% generate contrasts for main effects and interactions of a factorial design +% +% [cMatrix, cName] = contrasts(fLevel) +% [cMatrix, cName] = contrasts(fLevel, fName) +% +% fLevel: number of levels for each factor (row vector) +% fName: names of factors (cell array) +% cMatrix: contrast matrices (cell array) +% cName: contrast names (cell array) +% +% the first factor is the one being enumerated slowest. +% contrasts are not orthonormalized! +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +% number of factors +nf = size(fLevel, 2); + +if nargin < 2 + % generate generic names for factors + fName = num2cell(char((1 : nf) + '@')); +end + +% generate sorted list of contrast signatures +nc = 2 ^ nf - 1; +cs = bitget(ones(nf, 1) * (1 : nc), (1 : nf)' * ones(1, nc))'; +[~, ind] = sort(sum(cs, 2)); +cs = cs(ind, :); +cs = logical(cs); + +% compute contrast elements +e = cell(2, nf); +for fi = 1 : nf + e{1, fi} = ones(fLevel(fi), 1); + e{2, fi} = -(diff(eye(fLevel(fi)))'); +end + +% compute contrasts +cMatrix = cell(nf, 1); +cName = cell(nf, 1); +for ci = 1 : nc + % contrast matrix + cMatrix{ci} = 1; + for fi = 1 : nf + cMatrix{ci} = kron(cMatrix{ci}, e{cs(ci, fi) + 1, fi}); + end + % contrast name + l = sprintf('×%s', fName{cs(ci, :)}); + cName{ci} = l(2:end); +end diff --git a/cvManovaSearchlight.m b/cvManovaSearchlight.m new file mode 100644 index 0000000..a85addd --- /dev/null +++ b/cvManovaSearchlight.m @@ -0,0 +1,119 @@ +function cvManovaSearchlight(dirName, slRadius, Cs, permute) + +% cross-validated MANOVA on searchlight +% +% cvManovaSearchlight(dirName, slRadius, Cs, permute = false) +% +% dirName: directory where the SPM.mat file referring to an estimated +% model is located +% slRadius: radius of the searchlight sphere, in voxels +% Cs: cell array of contrast matrices +% permute: whether to compute permutation values of the test statistic +% +% Output files are written to the same directory: +% spmD_C####_P####.nii: +% images of the pattern discriminability D +% contrast and permutation are identified by numbers +% VPSL.img an image of the number of voxels within each searchlight +% cms.mat a record of the analysis parameters +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + + +fprintf('\n\ncvManovaSearchlight\n\n') + +if nargin < 4 + permute = false; +end + +nContrasts = numel(Cs); +outpattern = 'spmD_C%04d_P%04d.nii'; + +if dirName(end) ~= '/', dirName = [dirName '/']; end + +% simplify saving images +wd = cd; +cd(dirName) + +% load data, design matrix etc. +[Y, X, mask, misc] = loadDataSPM(dirName); + +% determine per-run design and data matrices +nRuns = misc.m; +Xrun = cell(nRuns, 1); +Yrun = cell(nRuns, 1); +% for each run +for ri = 1 : nRuns + Yrun{ri} = Y(misc.sRow{ri}, :); + Xrun{ri} = X(misc.sRow{ri}, misc.sCol{ri}); +end +clear Y X + +% check contrasts +for ci = 1 : nContrasts + if size(Cs{ci}, 2) > rank(Cs{ci}) + error('contrast %d is misspecified!', ci) + end + for ri = 1 : nRuns + if inestimability(Cs{ci}, Xrun{ri}) > 1e-6 + error('contrast %d is not estimable in run %d!', ci, ri) + end + end +end + +% precomputation +fprintf('\nprecomputing GLM runwise\n') +[XXs, betas, xis] = cvManova_precompute(Xrun, Yrun); +clear Xrun Yrun + +% determine voxels per searchlight, and save as image +fprintf('\ncomputing voxels per searchlight image\n') +p = runSearchlight(slRadius, mask, @(vi)(size(vi, 1))); +saveMRImage(reshape(p, size(mask)), 'VPSL.nii', misc.v2mm, ... + 'voxels per searchlight'); + +% error check +pMax = max(p(:)); +if pMax > 0.9 * (misc.m - 1) * misc.fE + error('insufficient amount of data for searchlight size %d!', pMax) +end + +% bias correction factor +factor = ((misc.m - 1) * misc.fE - p - 1) / ((misc.m - 1) * misc.n); +clear p + +% run searchlight +fprintf('\ncomputing cross-validated MANOVA on searchlight\n') +mDl = runSearchlight(slRadius, mask, @cvManova_compute, ... + XXs, betas, xis, Cs, permute); + +% separate contrast and permutation dimensions +nContrasts = numel(Cs); +nPerms = size(mDl, 2) / nContrasts; +mDl = reshape(mDl, [], nContrasts, nPerms); + +% compute the unbiased estimate of the pattern discriminability D +D = bsxfun(@times, factor, mDl); +clear mDl + +% save results +for ci = 1 : nContrasts + for pi = 1 : nPerms + fn = sprintf(outpattern, ci, pi); + saveMRImage(reshape(D(:, ci, pi), size(mask)), fn, misc.v2mm, ... + 'pattern discriminability'); + end +end + +% analysis parameters +save cms.mat slRadius Cs permute misc nPerms + +cd(wd) diff --git a/cvManova_compute.m b/cvManova_compute.m new file mode 100644 index 0000000..75c9dac --- /dev/null +++ b/cvManova_compute.m @@ -0,0 +1,139 @@ +function mDl = cvManova_compute(vi, XXs, betas, xis, Cs, permute) + +% cross-validated MANOVA, main computation +% +% mDl = cvManova_compute(vi, XXs, betas, xis, CC, permute) +% +% vi: indices of voxels to use +% XXs: cell array of per-run precomputed X' X +% betas: cell array of per-run precomputed GLM parameter estimates +% xis: cell array of per-run precomputed GLM residuals +% CC: cell array of contrast matrices +% permute: whether to compute permutations +% mDl: cross-validated MANOVA statistic, raw version +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +nRuns = size(betas, 1); + + +persistent sp % sign permutations +persistent CCs % canonical contrast matrices + +if isempty(vi) + % initialization call + + % generate sign permutations + sp = signPermutations(nRuns); + nPerms = size(sp, 2) / 2; % the two halves are equivalent + if ~permute + nPerms = 1; % identity permutation only + end + sp = sp(:, 1 : nPerms); + + % prepare canonical contrast matrices + nContrasts = numel(Cs); + CCs = cell(nContrasts, nRuns); + for ci = 1 : nContrasts + % compute canonical contrast matrix + CC = pinv(Cs{ci}') * Cs{ci}'; + + for k = 1 : nRuns + % zero-pad contrast to number of regressors + CCs{ci, k} = CC; + nReg = size(betas{k}, 1); + if size(CC, 1) < nReg + CCs{ci, k}(nReg, nReg) = 0; + end + end + end + + % answer initialization call + mDl = nan(1, nContrasts * nPerms); + return +else + nPerms = size(sp, 2); + nContrasts = size(CCs, 1); +end + + +% pre-compute per-run E +Es = cell(nRuns, 1); +for k = 1 : nRuns + x = xis{k}(:, vi); % weirdly, this is faster than + Es{k} = x' * x; % Es{k} = xis{k}(:, vi)' * xis{k}(:, vi); +end +clear x + +% pre-compute inverse of per-fold summed E +iEls = cell(nRuns, 1); +for l = 1 : nRuns + ks = 1 : nRuns; + ks(l) = []; + + El = sum(cat(3, Es{ks}), 3); + + iEls{l} = eye(size(El)) / El; % faster than inv +end +clear Es + + +mDl = zeros(nContrasts, nPerms); + +% for each contrast +for ci = 1 : nContrasts + + % pre-compute per-run betaDelta + betaDelta = cell(nRuns, 1); + for k = 1 : nRuns + betaDelta{k} = CCs{ci, k} * betas{k}(:, vi); + end + + % pre-compute per-run H + Hs = cell(nRuns, nRuns); + for k = 1 : nRuns + for l = 1 : nRuns + if l == k, continue, end + + Hs{k, l} = betaDelta{k}' * XXs{l} * betaDelta{l}; + end + end + clear betaDelta + + % for each permutation + for pi = 1 : nPerms + + % for each cross-validation fold + for l = 1 : nRuns + ks = 1 : nRuns; + ks(l) = []; + + % sign-permuted, summed H + Hl = sum(bsxfun(@times, ... + cat(3, Hs{ks, l}), ... + reshape(sp(ks, pi) * sp(l, pi)', 1, 1, nRuns - 1) ... + ), 3); + + % fold-wise D + Dl = sum(sum(Hl' .* iEls{l})); % faster than trace(Hl * iEls{l}) + + % sum across cross-validation folds + mDl(ci, pi) = mDl(ci, pi) + Dl; + end + end +end + +% sum -> mean +mDl = mDl / nRuns; + +% return row vector +mDl = mDl(:)'; + diff --git a/cvManova_precompute.m b/cvManova_precompute.m new file mode 100644 index 0000000..f90b2ff --- /dev/null +++ b/cvManova_precompute.m @@ -0,0 +1,40 @@ +function [XXs, betas, xis] = cvManova_precompute(Xrun, Yrun) + +% cross-validated MANOVA, GLM precomputation +% +% [XXs, betas, xis] = cvManova_precompute(Xrun, Yrun) +% +% Xrun: cell array of per-run design matrices X +% Yrun: cell array of per-run data matrices Y +% XXs: cell array of per-run X' X +% betas: cell array of per-run GLM parameter estimates +% xis: cell array of per-run GLM residuals +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +Xrun = Xrun(:); +Yrun = Yrun(:); +nRuns = size(Xrun, 1); + +betas = cell(nRuns, 1); +xis = cell(nRuns, 1); +XXs = cell(nRuns, 1); + +% for each run +for ri = 1 : nRuns + % estimate GLM + beta = pinv(Xrun{ri}) * Yrun{ri}; + xi = Yrun{ri} - Xrun{ri} * beta; + % store precomputed matrices + betas{ri} = beta; + xis{ri} = xi; + XXs{ri} = Xrun{ri}' * Xrun{ri}; +end diff --git a/inestimability.m b/inestimability.m new file mode 100644 index 0000000..ddbc0c8 --- /dev/null +++ b/inestimability.m @@ -0,0 +1,41 @@ +function ie = inestimability(C, X) + +% degree of inestimability of a contrast wrt a design matrix +% +% ie = inestimability(C, X) +% +% C: contrast matrix +% X: design matrix +% +% For an estimable contrast, ie should be 0 save for numerical error. +% A number smaller than 1 indicates that the contrast has an estimable +% part. For a completely inestimable contrast, ie = 1. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +% maximum "0" observed so far: 4.79*eps + +if size(C, 1) < size(X, 2) + C(size(X, 2), end) = 0; +end + +% determine orthonormal basis of the range of the contrast matrix +RC = orth(C); + +% determine orthonormal basis of the null space of the design matrix +NX = null(X); + +% NX' * RC is a matrix of the inner products ("correlations") +% of all the C-range vectors with all the X-null vectors. +% Squared values indicate the proportion of variance of an effect described +% by a C-range vector that falls within the null-space of X. +% The worst-case proportion across the range of C is given by the matrix norm: +ie = norm(NX' * RC); diff --git a/loadDataSPM.m b/loadDataSPM.m new file mode 100644 index 0000000..2a3683d --- /dev/null +++ b/loadDataSPM.m @@ -0,0 +1,121 @@ +function [Y, X, mask, misc] = loadDataSPM(dirName) + +% load fMRI data via SPM.mat +% +% [Y, X, mask, misc] = loadDataSPM(dirName) +% +% dirName: name of directory that contains SPM.mat +% Y: MR data (within mask), samples x voxels +% X: design matrix, samples x regressors +% mask: brain mask, logical 3d-volume +% misc: struct with additional data: +% v2mm voxels to mm transformation matrix +% sRow rows for each session +% sCol columns for each session +% m number of sessions +% n number of images per session +% fE residual degrees of freedom per session +% +% Y & X and are high-pass filtered and whitened. +% Y includes only those voxels selected through mask. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +SPMname = [dirName 'SPM.mat']; + +fprintf('loading data\n') +fprintf(' via %s\n', SPMname) + +% load SPM.mat +load(SPMname, 'SPM'); + +% recreate volumes, for backwards compatibility (with SPM2?) +fnames = {SPM.xY.VY.fname}; +if ~exist(fnames{1}, 'file') + % and to enable reading moved data files + fnames = patchPath(fnames, dirName); +end +fprintf(' reading volume information\n') +VY = cell2mat(spm_vol(fnames)); +% but copy scaling information +[VY(:).pinfo] = SPM.xY.VY(:).pinfo; +nImages = numel(VY); +nVoxels = prod(VY(1).dim); + +% check dimensions and orientations of all images +spm_check_orientations(VY); + +% read mask image +if isfield(SPM, 'VM') + VM = spm_vol([dirName SPM.VM.fname]); +else + VM = spm_vol([dirName 'mask.img']); +end +mask = logical(spm_read_vols(VM)); +nMaskVoxels = sum(mask(:)); + +% check memory +memNeed = nMaskVoxels * nImages * 8 / 1024 / 1024; % assuming double precision +memFree = systemFree / 1024; +fprintf(' memory needed %7.1f MiB\n', memNeed) +fprintf(' free %7.1f MiB\n', memFree) +if memFree < memNeed, warning('not enough memory!'), end + +% read and mask data +Y = nan(nImages, nMaskVoxels); +fprintf(' reading images\n') +for i = 1 : nImages % slow, but saves memory + y = spm_read_vols(VY(i)); + Y(i, :) = y(mask)'; % apply mask + if mod(i, 50) == 0 + fprintf(' %d of %d images read\n', i, nImages) + end +end +fprintf('\n') + +% get whitening matrix +if isfield(SPM.xX, 'W') + W = SPM.xX.W; +else + fprintf(' * SPM.mat does not define whitening matrix!\n') + W = 1; +end + +% whiten and high-pass filter data +fprintf(' filtering and whitening data\n') +Y = spm_filter(SPM.xX.K, W * Y); + +% get whitened and high-pass filtered design matrix +if isfield(SPM.xX, 'xKXs') + X = SPM.xX.xKXs.X; +else + fprintf(' * SPM.mat does not contain whitened and high-pass filtered design matrix!\n') + X = spm_filter(SPM.xX.K, W * SPM.xX.X); +end + +% degrees of freedom +Tdf = sum(SPM.nscan); % total +Kdf = sum(arrayfun(@(x)(size(x.X0, 2)), SPM.xX.K)); % loss from high-pass filter +Xdf = rank(X); % loss from regressors +Rdf = Tdf - Kdf - Xdf; % residual +fprintf(' df: %d - %d - %d = %d\n', Tdf, Kdf, Xdf, Rdf); + +% output +misc.v2mm = VY(1).mat; +misc.m = size(SPM.nscan, 2); +misc.sRow = {SPM.Sess.row}; +misc.sCol = {SPM.Sess.col}; +for si = 1 : misc.m % add constant regressors + misc.sCol{si} = [misc.sCol{si}, SPM.xX.iB(si)]; + % NOT SURE WHETHER THIS IS GENERALLY CORRECT +end +misc.fE = Rdf / misc.m; % if not consistent across sessions, +misc.n = Tdf / misc.m; % then this is an approximation diff --git a/patchPath.m b/patchPath.m new file mode 100644 index 0000000..c6725e6 --- /dev/null +++ b/patchPath.m @@ -0,0 +1,44 @@ +function fnames = patchPath(fnames, dirName) + +% patch the path portion of filenames to access them after being moved +% +% fnames = patchPath(fnames, dirname) +% +% Assumption: SPM.mat and image files resided in subfolders of some common +% folder and were moved together. In moving, at least one element of the +% folder hierarchy common to SPM.mat and image files was preserved. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +% find common beginning of all filenames +f = cell2mat(fnames(:)); +ind = find(var(f) > 0, 1, 'first'); +b = f(1, 1 : ind - 1); + +% decompose into folder hierarchy +bp = regexp(b, '[^/]*/', 'match'); +dp = regexp(dirName, '[^/]*/', 'match'); + +% use last common element to determine from and to root folders +[~, bpi, dpi] = intersect(bp, dp); +if ~isempty(bpi) + from = [bp{1 : max(bpi)}]; + to = [dp{1 : max(dpi)}]; +else + from = [bp{:}]; + to = [dp{:}]; +end + +% patch folders +fprintf(' patching paths\n from %s\n to %s\n', from, to) +f = [repmat(to, size(f, 1), 1), f(:, size(from, 2) + 1 : end)]; +fnames = cellstr(f); + diff --git a/runSearchlight.m b/runSearchlight.m new file mode 100644 index 0000000..b5415f9 --- /dev/null +++ b/runSearchlight.m @@ -0,0 +1,93 @@ +function res = runSearchlight(slRadius, mask, fun, varargin) + +% general purpose searchlight: +% apply function to data contained in a sliding spherical window +% +% res = runSearchlight(slRadius, mask, fun, ...) +% +% slRadius: radius of searchlight, in voxels +% mask: 3-dimensional logical array indicating which voxels to use +% fun: function to call with voxel indices within window +% – additional arguments are passed through to the function – +% res: results, 2-dimensional matrix, voxels x output values +% +% The function has to be of the form r = fun(mvi, ...) +% mvi: column vector of linear indices into the mask voxels +% r: row vector of results +% The output of fun([]) is used to initialize res. +% +% A voxel is included in the searchlight if its distance from the center is +% *smaller than or equal to* the radius. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +dim = size(mask); % volume dimensions +nVolumeVoxels = prod(dim); +nMaskVoxels = sum(mask(:)); + +% determine searchlight voxel offsets relative to center voxel +% prototype searchlight +[dxi, dyi, dzi] = ndgrid(-ceil(slRadius) : ceil(slRadius)); +PSL = (dxi .^ 2 + dyi .^ 2 + dzi .^ 2 <= slRadius .^ 2); +% spatial offsets +dxi = dxi(PSL); +dyi = dyi(PSL); +dzi = dzi(PSL); +% index offsets +PSL(dim(1), dim(2), dim(3)) = 0; +di = find(PSL); +cInd = find((dxi == 0) & (dyi == 0) & (dzi == 0)); +di = di - di(cInd); %#ok +clear PSL cInd + +% mapping from volume to mask voxel indices +vvi2mvi = nan(nVolumeVoxels, 1); +vvi2mvi(mask) = 1 : nMaskVoxels; + +% initialize result volume(s) +res = fun(nan(0, 1), varargin{:}); +res = repmat(reshape(res, 1, []), nVolumeVoxels, 1); + +fprintf(' running searchlight\n') +fprintf(' searchlight size: %d\n', size(di, 1)) +tic +t = 0; +cmvi = 0; +for cvvi = 1 : nVolumeVoxels % searchlight center volume voxel index + % process only if center is within mask + if mask(cvvi) + cmvi = cmvi + 1; % searchlight center mask voxel index + + % searchlight center coordinates + [xi, yi, zi] = ind2sub(dim, cvvi); + % searchlight voxel coordinates; limit to volume boundaries + ind = (xi + dxi >= 1) & (xi + dxi <= dim(1)) & ... + (yi + dyi >= 1) & (yi + dyi <= dim(2)) & ... + (zi + dzi >= 1) & (zi + dzi <= dim(3)); + % searchlight voxel volume indices + vvi = cvvi + di(ind); + % discard out-of-mask voxels + vvi = vvi(mask(vvi) == 1); + % translate to mask voxel indices + mvi = vvi2mvi(vvi); + + % call function and store output + res(cvvi, :) = fun(mvi, varargin{:}); + end + + % progress + nt = toc; + if (nt - t > 30) || (cvvi == nVolumeVoxels) + t = nt; + fprintf(' %6.1f min %6d voxels %5.1f %%\n', ... + t / 60, cmvi, cmvi / nMaskVoxels * 100) + end +end diff --git a/saveMRImage.m b/saveMRImage.m new file mode 100644 index 0000000..a03a4b6 --- /dev/null +++ b/saveMRImage.m @@ -0,0 +1,66 @@ +function V = saveMRImage(data, filename, mat, descrip) + +% saves a 3d-matrix to an MR image file, using SPM routines +% +% V = saveMRImage(data, filename, mat = eye(4), descrip = 'saveMRImage') +% +% data: image data to be written, 3d matrix +% filename: name of the file to write to +% mat: transformation from voxel indices to mm coordinates +% descrip: string describing contents +% V: spm volume struct +% +% The file format is determined by the filename extension. +% The data type is is determined by the class of data. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +if nargin < 3 + mat = eye(4); +end + +if nargin < 4 + descrip = 'saveMRImage'; +end + +if size(data, 4) > 1 + warning('only the first volume is written') +end + +if isfloat(data) && ~isreal(data) + warning('only real part of complex data is written') + data = real(data); +end + +% determine data type +dt = class(data); +if isa(data, 'single') + dt = 'float32'; +elseif isa(data, 'double') + dt = 'float64'; +end +dt = spm_type(dt); +if isnan(dt) + error('data type is not supported') +end + +V = struct; +V.fname = filename; +V.dim = size(data); +V.dt = [dt 0]; +V.mat = mat; +V.descrip = descrip; + +V = spm_write_vol(V, data); + +if nargout == 0 + clear V +end diff --git a/signPermutations.m b/signPermutations.m new file mode 100644 index 0000000..7f97fe0 --- /dev/null +++ b/signPermutations.m @@ -0,0 +1,35 @@ +function [perms, nPerms] = signPermutations(n, maxPerms) + +% generate sign permutations +% +% perms = signPermutations(n, maxPerms = 5000) +% +% Permutations are randomly selected if the full enumeration +% is larger than maxPerms. +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +if nargin < 2 + maxPerms = 5000; % adapted for z-values up to +-3 +end + +% determine permutations; first one is always the identity permutation +if 2 ^ n <= maxPerms + % full enumeration of permutations + nPerms = 2 ^ n; + perms = bitget(ones(n, 1) * (0 : nPerms - 1), (1 : n)' * ones(1, nPerms)); +else + % random (Monte Carlo) selection of permutations + nPerms = maxPerms; + perms = [zeros(n, 1), (rand(n, nPerms - 1) > 0.5)]; +end + +perms = (-1) .^ perms; % permuted signs diff --git a/systemFree.m b/systemFree.m new file mode 100644 index 0000000..10ccc7c --- /dev/null +++ b/systemFree.m @@ -0,0 +1,26 @@ +function kbytes = systemFree + +% returns the amount of free memory in units of 1024 bytes +% +% kbytes = systemFree() +% +% Copyright (C) 2013 Carsten Allefeld +% +% 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. + +if ispc + % Windows: largest possible size for a single new array. + user = memory; + kbytes = user.MaxPossibleArrayBytes / 1024; +else + % Works on Linux and probably other Unices. Macs, too? + % Memory that is available, possibly by flushing cache & buffers. + [~, out] = system('free -k'); + kbytes = str2double(out(178:188)); +end