diff --git a/.gitignore b/.gitignore
deleted file mode 100644
index 06fd450243217f4db170db1c64876470590cddcb..0000000000000000000000000000000000000000
--- a/.gitignore
+++ /dev/null
@@ -1,13 +0,0 @@
-.vscode
-.DS_Store
-
-exp/ICASSP/*.txt
-
-# Julia's .gitignore
-Manifest.toml
-*.jl.cov
-*.jl.*.cov
-*.jl.mem
-/docs/build/
-/docs/Manifest.toml
-
diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index ae2ef67bd2af50eb7977cf964d43db46c9feca6c..0000000000000000000000000000000000000000
--- a/LICENSE
+++ /dev/null
@@ -1,661 +0,0 @@
- GNU AFFERO GENERAL PUBLIC LICENSE
- Version 3, 19 November 2007
-
- Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
- Preamble
-
- The GNU Affero General Public License is a free, copyleft license for
-software and other kinds of works, specifically designed to ensure
-cooperation with the community in the case of network server software.
-
- The licenses for most software and other practical works are designed
-to take away your freedom to share and change the works. By contrast,
-our General Public Licenses are intended to guarantee your freedom to
-share and change all versions of a program--to make sure it remains free
-software for all its users.
-
- When we speak of free software, we are referring to freedom, not
-price. Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-them if you wish), that you receive source code or can get it if you
-want it, that you can change the software or use pieces of it in new
-free programs, and that you know you can do these things.
-
- Developers that use our General Public Licenses protect your rights
-with two steps: (1) assert copyright on the software, and (2) offer
-you this License which gives you legal permission to copy, distribute
-and/or modify the software.
-
- A secondary benefit of defending all users' freedom is that
-improvements made in alternate versions of the program, if they
-receive widespread use, become available for other developers to
-incorporate. Many developers of free software are heartened and
-encouraged by the resulting cooperation. However, in the case of
-software used on network servers, this result may fail to come about.
-The GNU General Public License permits making a modified version and
-letting the public access it on a server without ever releasing its
-source code to the public.
-
- The GNU Affero General Public License is designed specifically to
-ensure that, in such cases, the modified source code becomes available
-to the community. It requires the operator of a network server to
-provide the source code of the modified version running there to the
-users of that server. Therefore, public use of a modified version, on
-a publicly accessible server, gives the public access to the source
-code of the modified version.
-
- An older license, called the Affero General Public License and
-published by Affero, was designed to accomplish similar goals. This is
-a different license, not a version of the Affero GPL, but Affero has
-released a new version of the Affero GPL which permits relicensing under
-this license.
-
- The precise terms and conditions for copying, distribution and
-modification follow.
-
- TERMS AND CONDITIONS
-
- 0. Definitions.
-
- "This License" refers to version 3 of the GNU Affero General Public License.
-
- "Copyright" also means copyright-like laws that apply to other kinds of
-works, such as semiconductor masks.
-
- "The Program" refers to any copyrightable work licensed under this
-License. Each licensee is addressed as "you". "Licensees" and
-"recipients" may be individuals or organizations.
-
- To "modify" a work means to copy from or adapt all or part of the work
-in a fashion requiring copyright permission, other than the making of an
-exact copy. The resulting work is called a "modified version" of the
-earlier work or a work "based on" the earlier work.
-
- A "covered work" means either the unmodified Program or a work based
-on the Program.
-
- To "propagate" a work means to do anything with it that, without
-permission, would make you directly or secondarily liable for
-infringement under applicable copyright law, except executing it on a
-computer or modifying a private copy. Propagation includes copying,
-distribution (with or without modification), making available to the
-public, and in some countries other activities as well.
-
- To "convey" a work means any kind of propagation that enables other
-parties to make or receive copies. Mere interaction with a user through
-a computer network, with no transfer of a copy, is not conveying.
-
- An interactive user interface displays "Appropriate Legal Notices"
-to the extent that it includes a convenient and prominently visible
-feature that (1) displays an appropriate copyright notice, and (2)
-tells the user that there is no warranty for the work (except to the
-extent that warranties are provided), that licensees may convey the
-work under this License, and how to view a copy of this License. If
-the interface presents a list of user commands or options, such as a
-menu, a prominent item in the list meets this criterion.
-
- 1. Source Code.
-
- The "source code" for a work means the preferred form of the work
-for making modifications to it. "Object code" means any non-source
-form of a work.
-
- A "Standard Interface" means an interface that either is an official
-standard defined by a recognized standards body, or, in the case of
-interfaces specified for a particular programming language, one that
-is widely used among developers working in that language.
-
- The "System Libraries" of an executable work include anything, other
-than the work as a whole, that (a) is included in the normal form of
-packaging a Major Component, but which is not part of that Major
-Component, and (b) serves only to enable use of the work with that
-Major Component, or to implement a Standard Interface for which an
-implementation is available to the public in source code form. A
-"Major Component", in this context, means a major essential component
-(kernel, window system, and so on) of the specific operating system
-(if any) on which the executable work runs, or a compiler used to
-produce the work, or an object code interpreter used to run it.
-
- The "Corresponding Source" for a work in object code form means all
-the source code needed to generate, install, and (for an executable
-work) run the object code and to modify the work, including scripts to
-control those activities. However, it does not include the work's
-System Libraries, or general-purpose tools or generally available free
-programs which are used unmodified in performing those activities but
-which are not part of the work. For example, Corresponding Source
-includes interface definition files associated with source files for
-the work, and the source code for shared libraries and dynamically
-linked subprograms that the work is specifically designed to require,
-such as by intimate data communication or control flow between those
-subprograms and other parts of the work.
-
- The Corresponding Source need not include anything that users
-can regenerate automatically from other parts of the Corresponding
-Source.
-
- The Corresponding Source for a work in source code form is that
-same work.
-
- 2. Basic Permissions.
-
- All rights granted under this License are granted for the term of
-copyright on the Program, and are irrevocable provided the stated
-conditions are met. This License explicitly affirms your unlimited
-permission to run the unmodified Program. The output from running a
-covered work is covered by this License only if the output, given its
-content, constitutes a covered work. This License acknowledges your
-rights of fair use or other equivalent, as provided by copyright law.
-
- You may make, run and propagate covered works that you do not
-convey, without conditions so long as your license otherwise remains
-in force. You may convey covered works to others for the sole purpose
-of having them make modifications exclusively for you, or provide you
-with facilities for running those works, provided that you comply with
-the terms of this License in conveying all material for which you do
-not control copyright. Those thus making or running the covered works
-for you must do so exclusively on your behalf, under your direction
-and control, on terms that prohibit them from making any copies of
-your copyrighted material outside their relationship with you.
-
- Conveying under any other circumstances is permitted solely under
-the conditions stated below. Sublicensing is not allowed; section 10
-makes it unnecessary.
-
- 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
-
- No covered work shall be deemed part of an effective technological
-measure under any applicable law fulfilling obligations under article
-11 of the WIPO copyright treaty adopted on 20 December 1996, or
-similar laws prohibiting or restricting circumvention of such
-measures.
-
- When you convey a covered work, you waive any legal power to forbid
-circumvention of technological measures to the extent such circumvention
-is effected by exercising rights under this License with respect to
-the covered work, and you disclaim any intention to limit operation or
-modification of the work as a means of enforcing, against the work's
-users, your or third parties' legal rights to forbid circumvention of
-technological measures.
-
- 4. Conveying Verbatim Copies.
-
- You may convey verbatim copies of the Program's source code as you
-receive it, in any medium, provided that you conspicuously and
-appropriately publish on each copy an appropriate copyright notice;
-keep intact all notices stating that this License and any
-non-permissive terms added in accord with section 7 apply to the code;
-keep intact all notices of the absence of any warranty; and give all
-recipients a copy of this License along with the Program.
-
- You may charge any price or no price for each copy that you convey,
-and you may offer support or warranty protection for a fee.
-
- 5. Conveying Modified Source Versions.
-
- You may convey a work based on the Program, or the modifications to
-produce it from the Program, in the form of source code under the
-terms of section 4, provided that you also meet all of these conditions:
-
- a) The work must carry prominent notices stating that you modified
- it, and giving a relevant date.
-
- b) The work must carry prominent notices stating that it is
- released under this License and any conditions added under section
- 7. This requirement modifies the requirement in section 4 to
- "keep intact all notices".
-
- c) You must license the entire work, as a whole, under this
- License to anyone who comes into possession of a copy. This
- License will therefore apply, along with any applicable section 7
- additional terms, to the whole of the work, and all its parts,
- regardless of how they are packaged. This License gives no
- permission to license the work in any other way, but it does not
- invalidate such permission if you have separately received it.
-
- d) If the work has interactive user interfaces, each must display
- Appropriate Legal Notices; however, if the Program has interactive
- interfaces that do not display Appropriate Legal Notices, your
- work need not make them do so.
-
- A compilation of a covered work with other separate and independent
-works, which are not by their nature extensions of the covered work,
-and which are not combined with it such as to form a larger program,
-in or on a volume of a storage or distribution medium, is called an
-"aggregate" if the compilation and its resulting copyright are not
-used to limit the access or legal rights of the compilation's users
-beyond what the individual works permit. Inclusion of a covered work
-in an aggregate does not cause this License to apply to the other
-parts of the aggregate.
-
- 6. Conveying Non-Source Forms.
-
- You may convey a covered work in object code form under the terms
-of sections 4 and 5, provided that you also convey the
-machine-readable Corresponding Source under the terms of this License,
-in one of these ways:
-
- a) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by the
- Corresponding Source fixed on a durable physical medium
- customarily used for software interchange.
-
- b) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by a
- written offer, valid for at least three years and valid for as
- long as you offer spare parts or customer support for that product
- model, to give anyone who possesses the object code either (1) a
- copy of the Corresponding Source for all the software in the
- product that is covered by this License, on a durable physical
- medium customarily used for software interchange, for a price no
- more than your reasonable cost of physically performing this
- conveying of source, or (2) access to copy the
- Corresponding Source from a network server at no charge.
-
- c) Convey individual copies of the object code with a copy of the
- written offer to provide the Corresponding Source. This
- alternative is allowed only occasionally and noncommercially, and
- only if you received the object code with such an offer, in accord
- with subsection 6b.
-
- d) Convey the object code by offering access from a designated
- place (gratis or for a charge), and offer equivalent access to the
- Corresponding Source in the same way through the same place at no
- further charge. You need not require recipients to copy the
- Corresponding Source along with the object code. If the place to
- copy the object code is a network server, the Corresponding Source
- may be on a different server (operated by you or a third party)
- that supports equivalent copying facilities, provided you maintain
- clear directions next to the object code saying where to find the
- Corresponding Source. Regardless of what server hosts the
- Corresponding Source, you remain obligated to ensure that it is
- available for as long as needed to satisfy these requirements.
-
- e) Convey the object code using peer-to-peer transmission, provided
- you inform other peers where the object code and Corresponding
- Source of the work are being offered to the general public at no
- charge under subsection 6d.
-
- A separable portion of the object code, whose source code is excluded
-from the Corresponding Source as a System Library, need not be
-included in conveying the object code work.
-
- A "User Product" is either (1) a "consumer product", which means any
-tangible personal property which is normally used for personal, family,
-or household purposes, or (2) anything designed or sold for incorporation
-into a dwelling. In determining whether a product is a consumer product,
-doubtful cases shall be resolved in favor of coverage. For a particular
-product received by a particular user, "normally used" refers to a
-typical or common use of that class of product, regardless of the status
-of the particular user or of the way in which the particular user
-actually uses, or expects or is expected to use, the product. A product
-is a consumer product regardless of whether the product has substantial
-commercial, industrial or non-consumer uses, unless such uses represent
-the only significant mode of use of the product.
-
- "Installation Information" for a User Product means any methods,
-procedures, authorization keys, or other information required to install
-and execute modified versions of a covered work in that User Product from
-a modified version of its Corresponding Source. The information must
-suffice to ensure that the continued functioning of the modified object
-code is in no case prevented or interfered with solely because
-modification has been made.
-
- If you convey an object code work under this section in, or with, or
-specifically for use in, a User Product, and the conveying occurs as
-part of a transaction in which the right of possession and use of the
-User Product is transferred to the recipient in perpetuity or for a
-fixed term (regardless of how the transaction is characterized), the
-Corresponding Source conveyed under this section must be accompanied
-by the Installation Information. But this requirement does not apply
-if neither you nor any third party retains the ability to install
-modified object code on the User Product (for example, the work has
-been installed in ROM).
-
- The requirement to provide Installation Information does not include a
-requirement to continue to provide support service, warranty, or updates
-for a work that has been modified or installed by the recipient, or for
-the User Product in which it has been modified or installed. Access to a
-network may be denied when the modification itself materially and
-adversely affects the operation of the network or violates the rules and
-protocols for communication across the network.
-
- Corresponding Source conveyed, and Installation Information provided,
-in accord with this section must be in a format that is publicly
-documented (and with an implementation available to the public in
-source code form), and must require no special password or key for
-unpacking, reading or copying.
-
- 7. Additional Terms.
-
- "Additional permissions" are terms that supplement the terms of this
-License by making exceptions from one or more of its conditions.
-Additional permissions that are applicable to the entire Program shall
-be treated as though they were included in this License, to the extent
-that they are valid under applicable law. If additional permissions
-apply only to part of the Program, that part may be used separately
-under those permissions, but the entire Program remains governed by
-this License without regard to the additional permissions.
-
- When you convey a copy of a covered work, you may at your option
-remove any additional permissions from that copy, or from any part of
-it. (Additional permissions may be written to require their own
-removal in certain cases when you modify the work.) You may place
-additional permissions on material, added by you to a covered work,
-for which you have or can give appropriate copyright permission.
-
- Notwithstanding any other provision of this License, for material you
-add to a covered work, you may (if authorized by the copyright holders of
-that material) supplement the terms of this License with terms:
-
- a) Disclaiming warranty or limiting liability differently from the
- terms of sections 15 and 16 of this License; or
-
- b) Requiring preservation of specified reasonable legal notices or
- author attributions in that material or in the Appropriate Legal
- Notices displayed by works containing it; or
-
- c) Prohibiting misrepresentation of the origin of that material, or
- requiring that modified versions of such material be marked in
- reasonable ways as different from the original version; or
-
- d) Limiting the use for publicity purposes of names of licensors or
- authors of the material; or
-
- e) Declining to grant rights under trademark law for use of some
- trade names, trademarks, or service marks; or
-
- f) Requiring indemnification of licensors and authors of that
- material by anyone who conveys the material (or modified versions of
- it) with contractual assumptions of liability to the recipient, for
- any liability that these contractual assumptions directly impose on
- those licensors and authors.
-
- All other non-permissive additional terms are considered "further
-restrictions" within the meaning of section 10. If the Program as you
-received it, or any part of it, contains a notice stating that it is
-governed by this License along with a term that is a further
-restriction, you may remove that term. If a license document contains
-a further restriction but permits relicensing or conveying under this
-License, you may add to a covered work material governed by the terms
-of that license document, provided that the further restriction does
-not survive such relicensing or conveying.
-
- If you add terms to a covered work in accord with this section, you
-must place, in the relevant source files, a statement of the
-additional terms that apply to those files, or a notice indicating
-where to find the applicable terms.
-
- Additional terms, permissive or non-permissive, may be stated in the
-form of a separately written license, or stated as exceptions;
-the above requirements apply either way.
-
- 8. Termination.
-
- You may not propagate or modify a covered work except as expressly
-provided under this License. Any attempt otherwise to propagate or
-modify it is void, and will automatically terminate your rights under
-this License (including any patent licenses granted under the third
-paragraph of section 11).
-
- However, if you cease all violation of this License, then your
-license from a particular copyright holder is reinstated (a)
-provisionally, unless and until the copyright holder explicitly and
-finally terminates your license, and (b) permanently, if the copyright
-holder fails to notify you of the violation by some reasonable means
-prior to 60 days after the cessation.
-
- Moreover, your license from a particular copyright holder is
-reinstated permanently if the copyright holder notifies you of the
-violation by some reasonable means, this is the first time you have
-received notice of violation of this License (for any work) from that
-copyright holder, and you cure the violation prior to 30 days after
-your receipt of the notice.
-
- Termination of your rights under this section does not terminate the
-licenses of parties who have received copies or rights from you under
-this License. If your rights have been terminated and not permanently
-reinstated, you do not qualify to receive new licenses for the same
-material under section 10.
-
- 9. Acceptance Not Required for Having Copies.
-
- You are not required to accept this License in order to receive or
-run a copy of the Program. Ancillary propagation of a covered work
-occurring solely as a consequence of using peer-to-peer transmission
-to receive a copy likewise does not require acceptance. However,
-nothing other than this License grants you permission to propagate or
-modify any covered work. These actions infringe copyright if you do
-not accept this License. Therefore, by modifying or propagating a
-covered work, you indicate your acceptance of this License to do so.
-
- 10. Automatic Licensing of Downstream Recipients.
-
- Each time you convey a covered work, the recipient automatically
-receives a license from the original licensors, to run, modify and
-propagate that work, subject to this License. You are not responsible
-for enforcing compliance by third parties with this License.
-
- An "entity transaction" is a transaction transferring control of an
-organization, or substantially all assets of one, or subdividing an
-organization, or merging organizations. If propagation of a covered
-work results from an entity transaction, each party to that
-transaction who receives a copy of the work also receives whatever
-licenses to the work the party's predecessor in interest had or could
-give under the previous paragraph, plus a right to possession of the
-Corresponding Source of the work from the predecessor in interest, if
-the predecessor has it or can get it with reasonable efforts.
-
- You may not impose any further restrictions on the exercise of the
-rights granted or affirmed under this License. For example, you may
-not impose a license fee, royalty, or other charge for exercise of
-rights granted under this License, and you may not initiate litigation
-(including a cross-claim or counterclaim in a lawsuit) alleging that
-any patent claim is infringed by making, using, selling, offering for
-sale, or importing the Program or any portion of it.
-
- 11. Patents.
-
- A "contributor" is a copyright holder who authorizes use under this
-License of the Program or a work on which the Program is based. The
-work thus licensed is called the contributor's "contributor version".
-
- A contributor's "essential patent claims" are all patent claims
-owned or controlled by the contributor, whether already acquired or
-hereafter acquired, that would be infringed by some manner, permitted
-by this License, of making, using, or selling its contributor version,
-but do not include claims that would be infringed only as a
-consequence of further modification of the contributor version. For
-purposes of this definition, "control" includes the right to grant
-patent sublicenses in a manner consistent with the requirements of
-this License.
-
- Each contributor grants you a non-exclusive, worldwide, royalty-free
-patent license under the contributor's essential patent claims, to
-make, use, sell, offer for sale, import and otherwise run, modify and
-propagate the contents of its contributor version.
-
- In the following three paragraphs, a "patent license" is any express
-agreement or commitment, however denominated, not to enforce a patent
-(such as an express permission to practice a patent or covenant not to
-sue for patent infringement). To "grant" such a patent license to a
-party means to make such an agreement or commitment not to enforce a
-patent against the party.
-
- If you convey a covered work, knowingly relying on a patent license,
-and the Corresponding Source of the work is not available for anyone
-to copy, free of charge and under the terms of this License, through a
-publicly available network server or other readily accessible means,
-then you must either (1) cause the Corresponding Source to be so
-available, or (2) arrange to deprive yourself of the benefit of the
-patent license for this particular work, or (3) arrange, in a manner
-consistent with the requirements of this License, to extend the patent
-license to downstream recipients. "Knowingly relying" means you have
-actual knowledge that, but for the patent license, your conveying the
-covered work in a country, or your recipient's use of the covered work
-in a country, would infringe one or more identifiable patents in that
-country that you have reason to believe are valid.
-
- If, pursuant to or in connection with a single transaction or
-arrangement, you convey, or propagate by procuring conveyance of, a
-covered work, and grant a patent license to some of the parties
-receiving the covered work authorizing them to use, propagate, modify
-or convey a specific copy of the covered work, then the patent license
-you grant is automatically extended to all recipients of the covered
-work and works based on it.
-
- A patent license is "discriminatory" if it does not include within
-the scope of its coverage, prohibits the exercise of, or is
-conditioned on the non-exercise of one or more of the rights that are
-specifically granted under this License. You may not convey a covered
-work if you are a party to an arrangement with a third party that is
-in the business of distributing software, under which you make payment
-to the third party based on the extent of your activity of conveying
-the work, and under which the third party grants, to any of the
-parties who would receive the covered work from you, a discriminatory
-patent license (a) in connection with copies of the covered work
-conveyed by you (or copies made from those copies), or (b) primarily
-for and in connection with specific products or compilations that
-contain the covered work, unless you entered into that arrangement,
-or that patent license was granted, prior to 28 March 2007.
-
- Nothing in this License shall be construed as excluding or limiting
-any implied license or other defenses to infringement that may
-otherwise be available to you under applicable patent law.
-
- 12. No Surrender of Others' Freedom.
-
- If conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License. If you cannot convey a
-covered work so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you may
-not convey it at all. For example, if you agree to terms that obligate you
-to collect a royalty for further conveying from those to whom you convey
-the Program, the only way you could satisfy both those terms and this
-License would be to refrain entirely from conveying the Program.
-
- 13. Remote Network Interaction; Use with the GNU General Public License.
-
- Notwithstanding any other provision of this License, if you modify the
-Program, your modified version must prominently offer all users
-interacting with it remotely through a computer network (if your version
-supports such interaction) an opportunity to receive the Corresponding
-Source of your version by providing access to the Corresponding Source
-from a network server at no charge, through some standard or customary
-means of facilitating copying of software. This Corresponding Source
-shall include the Corresponding Source for any work covered by version 3
-of the GNU General Public License that is incorporated pursuant to the
-following paragraph.
-
- Notwithstanding any other provision of this License, you have
-permission to link or combine any covered work with a work licensed
-under version 3 of the GNU General Public License into a single
-combined work, and to convey the resulting work. The terms of this
-License will continue to apply to the part which is the covered work,
-but the work with which it is combined will remain governed by version
-3 of the GNU General Public License.
-
- 14. Revised Versions of this License.
-
- The Free Software Foundation may publish revised and/or new versions of
-the GNU Affero General Public License from time to time. Such new versions
-will be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
- Each version is given a distinguishing version number. If the
-Program specifies that a certain numbered version of the GNU Affero General
-Public License "or any later version" applies to it, you have the
-option of following the terms and conditions either of that numbered
-version or of any later version published by the Free Software
-Foundation. If the Program does not specify a version number of the
-GNU Affero General Public License, you may choose any version ever published
-by the Free Software Foundation.
-
- If the Program specifies that a proxy can decide which future
-versions of the GNU Affero General Public License can be used, that proxy's
-public statement of acceptance of a version permanently authorizes you
-to choose that version for the Program.
-
- Later license versions may give you additional or different
-permissions. However, no additional obligations are imposed on any
-author or copyright holder as a result of your choosing to follow a
-later version.
-
- 15. Disclaimer of Warranty.
-
- THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
-APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
-HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
-OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
-THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
-IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
-ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
-
- 16. Limitation of Liability.
-
- IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
-THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
-GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
-USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
-DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
-PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
-EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
-SUCH DAMAGES.
-
- 17. Interpretation of Sections 15 and 16.
-
- If the disclaimer of warranty and limitation of liability provided
-above cannot be given local legal effect according to their terms,
-reviewing courts shall apply local law that most closely approximates
-an absolute waiver of all civil liability in connection with the
-Program, unless a warranty or assumption of liability accompanies a
-copy of the Program in return for a fee.
-
- END OF TERMS AND CONDITIONS
-
- How to Apply These Terms to Your New Programs
-
- If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
- To do so, attach the following notices to the program. It is safest
-to attach them to the start of each source file to most effectively
-state the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
- BnB-screening : A Branch-and-Bound algorithm enhanced with node-screening tests
- Copyright (C) 2021 Theo Guyard
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU Affero General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Affero General Public License for more details.
-
- You should have received a copy of the GNU Affero General Public License
- along with this program. If not, see <https://www.gnu.org/licenses/>.
-
-Also add information on how to contact you by electronic and paper mail.
-
- If your software can interact with users remotely through a computer
-network, you should also make sure that it provides a way for users to
-get its source. For example, if your program is a web application, its
-interface could display a "Source" link that leads users to an archive
-of the code. There are many ways you could offer source, and different
-solutions will be better for different programs; see section 13 for the
-specific requirements.
-
- You should also get your employer (if you work as a programmer) or school,
-if any, to sign a "copyright disclaimer" for the program, if necessary.
-For more information on this, and how to apply and follow the GNU AGPL, see
-<https://www.gnu.org/licenses/>.
diff --git a/Project.toml b/Project.toml
deleted file mode 100644
index fcd9fed7b2b1d01cedcca1c8e2d876d0de901d0f..0000000000000000000000000000000000000000
--- a/Project.toml
+++ /dev/null
@@ -1,5 +0,0 @@
-[deps]
-CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
diff --git a/README.md b/README.md
index 5233a3caf23161043d3a187a2c9998e30b981bd3..1a5eb202404ed5d5db528931a107ac8064e3812b 100644
--- a/README.md
+++ b/README.md
@@ -1,82 +1 @@
-# Node-screening tests for the L0-penalized least-square problem
-
-This repository contains numerical procedures linked to the paper
-
-> Guyard, T., Herzet, C., & Elvira, C. (2022, May). Node-Screening Tests For The L0-Penalized Least-Squares Problem. In *ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* (pp. 5448-5452). IEEE.
-
-available [here](https://arxiv.org/pdf/2110.07308.pdf). If you encounter a bug or something unexpected, please let me know by raising an issue on the project page or by contacting me by [mail](mailto:theo.guyard@insa-rennes.fr).
-
-## Requirements
-
-This repository works with Julia v1.5+. Please refer to the [Julia download page](https://julialang.org/downloads/) for install instructions.
-
-**Dependencies :** [JuMP.jl](https://github.com/jump-dev/JuMP.jl), [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl), [Distributions.jl](https://github.com/JuliaStats/Distributions.jl), [StatsBase.jl](https://github.com/JuliaStats/StatsBase.jl)
-
-Our code works with the [CPLEX](https://www.ibm.com/fr-fr/analytics/cplex-optimizer) solver through [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl). You cannot use CPLEX without having purchased and installed a copy of CPLEX Optimization Studio from [IBM](https://www.ibm.com). However, CPLEX is available for free to [academics and students](https://community.ibm.com/community/user/datascience/blogs/xavier-nodet1/2020/07/09/cplex-free-for-students). We recommend you to follow carefully the install instructions described at [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl) to install the package.
-
-## Install instructions
-
-1. Clone the repository
-```bash
-git clone https://gitlab.insa-rennes.fr/Theo.Guyard/bnb-screening.git
-```
-
-2. Enter in the project root folder
-```bash
-cd bnb-screening
-```
-
-3. Ensure that Julia v1.5+ is installed
-```bash
-julia -v
-> julia version 1.5.3 (or higher)
-```
-
-4. Install dependencies using [Julia's Pkg](https://docs.julialang.org/en/v1/stdlib/Pkg/) in REPL mode
-```julia
-(@v1.5) pkg> activate .
-(bnb-screening) pkg> instantiate
-```
-It is possible that you also have to build [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl) manually with the command
-```julia
-(bnb-screening) pkg> build CPLEX
-```
-
-## ICASSP experiments
-
-Experiments presented in the paper submitted to [ICASSP 2022](https://2022.ieeeicassp.org) are located in the folder `bnb-screening/exp/ICASSP/`. The available experiments are
-* `performances.jl` corresponding to results of table 1
-
-Run experiments directly from the project root folder as follows.
-```bash
-julia --project=. -t 1 exp/ICASSP/performances.jl
-```
-
-The `--project=.` flag installs and builds dependencies if needed. The `-t 1` flag restricts computations to only 1 CPU core in order to avoid bias due to parallelization capabilities. Logs of the experiment are saved in the same folder as `performances.jl`. You can edit the experimental setup directly in the script, which is well commented.
-
-## Running demo examples
-
-The `bnb-screening/examples/` folder contains portions of code explaining how to use our package. You can run these scripts as follows.
-```bash
-julia --project=. examples/solve_bnb.jl
-julia --project=. examples/solve_mip.jl
-```
-
-## Licence
-
-This software is distributed under the [MIT Licence](https://mit-license.org).
-
-## Cite this work
-
-If you use this package for your own work, please consider citing it as :
-
-```bibtex
-@inproceedings{guyard2022node,
- title={Node-Screening Tests For The L0-Penalized Least-Squares Problem},
- author={Guyard, Th{\'e}o and Herzet, C{\'e}dric and Elvira, Cl{\'e}ment},
- booktitle={ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)},
- pages={5448--5452},
- year={2022},
- organization={IEEE}
-}
-```
+This repo has been migrated to [BnbScreening.jl](https://github.com/TheoGuyard/BnbScreening.jl) on GitHub.
diff --git a/examples/solve_bnb.jl b/examples/solve_bnb.jl
deleted file mode 100644
index 4f04cb56d527206221e54d5f0325aa800242c076..0000000000000000000000000000000000000000
--- a/examples/solve_bnb.jl
+++ /dev/null
@@ -1,32 +0,0 @@
-include("../src/BnbScreening.jl")
-
-# ---------------------------------------------------------------------------- #
-# Example : Solving a L0-penalized least-square problem using a BnB algorithm #
-# ---------------------------------------------------------------------------- #
-
-# Generate problem data (gaussian setup or toeplitz setup)
-k, m, n, snr = 5, 500, 300, 10*log10(10)
-A, y, λ, M = data_toeplitz(k,m,n,snr)
-#k, m, n, snr = 5, 500, 1000, 10*log10(10)
-#A, y, λ, M = data_gaussian(k,m,n,snr)
-
-# BnB parameters
-bnbparams = BnbParams(
- lb_method=ACTIVESET,
- # Relaxation solving method
- # - ACTIVESET : the relaxation is solved using an Active-Set algorithm
- # - LBMIP : the relaxation is solved using CPLEX convex solver
- ub_method=BVLS,
- # Heuristic method
- # - BVLS : the heuristic is solved using the Bounded-Variable least-square algorithm
- # - LBMIP : the heuristic is solved using CPLEX convex solver
- maxtime=1000, # Maximum solution time allowed
- tol=1e-8, # Integer tolerance, ie. |x|<=tol <=> x=0
- maxiter=1000, # Maximum number of iterations allowed in the Active-Set method
- screening=false, # Toogle branch-screening rules
- verbosity=true, # Toogle verbosity
- showevery=1, # If verbosity=true, number of nodes between two displays
-)
-
-# Solve the problem using the BnB algorithm with appropriate parameters
-result_bnb = solve_bnb(A,y,λ,M,bnbparams=bnbparams)
diff --git a/examples/solve_mip.jl b/examples/solve_mip.jl
deleted file mode 100644
index 5af1f4e23974c747cf1404902365a6e6a912a493..0000000000000000000000000000000000000000
--- a/examples/solve_mip.jl
+++ /dev/null
@@ -1,27 +0,0 @@
-include("../src/BnbScreening.jl")
-
-# ---------------------------------------------------------------------------- #
-# Example : Solving a L0-penalized least-square problem using a MIP solver #
-# ---------------------------------------------------------------------------- #
-
-# Generate problem data (gaussian setup or toeplitz setup)
-k, m, n, snr = 5, 500, 300, 10*log10(10)
-A, y, λ, M = data_toeplitz(k,m,n,snr)
-#k, m, n, snr = 5, 500, 1000, 10*log10(10)
-#A, y, λ, M = data_gaussian(k,m,n,snr)
-
-# MIP solver parameters
-mipparams = MipParams(
- maxtime = 1000, # Maximum solution time in seconds
- silent = false, # Toogle solver verbosity
-)
-
-# Solve the problem using a MIP solver with appropriate parameters
-result_mip = solve_mip(A,y,λ,M,mipparams=mipparams)
-
-# Display solver results and statistics
-println("\nMIP")
-println(" Status : $(result_mip.termination_status)")
-println(" Objective : $(result_mip.objective_value)")
-println(" Nodes : $(result_mip.node_count)")
-println(" Timer : $(result_mip.solve_time)")
\ No newline at end of file
diff --git a/exp/ICASSP/performances.jl b/exp/ICASSP/performances.jl
deleted file mode 100644
index 771abcbccd1dedb95562b4b9f7998c74cbae3c95..0000000000000000000000000000000000000000
--- a/exp/ICASSP/performances.jl
+++ /dev/null
@@ -1,173 +0,0 @@
-# ICASSP experient - Performances
-
-include("../../src/BnbScreening.jl")
-using StatsBase
-
-################################ Configuration #################################
-
-# This is the only part of code you may have to edit
-
-setup = "toeplitz" # Setup to use : "gaussian" or "toeplitz"
-k = 5 # Sparsity level
-m, n = 500, 300 # Dictionnary dimension
-snr = 10. * log10(10.) # SNR of the observation noise
-maxtime = 1000 # Maximum solution time (in second) allowed
-maxiter = 1000 # Maximum Active-Set iterations allowed per relaxation resolution
-tol = 1.e-8 # Integer tolerance (ie, |x|<tol <=> x=0)
-repeats = 10 # Number of repeats in the experiment
-
-
-################################## Experiment ##################################
-
-function print_and_log(s::String, f::IOStream)
- print(s)
- write(f,s)
-end
-
-function precompilation(
- setup::String,
- k::Int,
- m::Int,
- n::Int,
- snr::Float64,
- maxtime::Int,
- maxiter::Int,
- tol::Float64,
- )
- print_and_log("Precompiling methods ...\n",logfile)
- try
- if setup == "gaussian"
- A, y, λ, M = data_gaussian(k,m,n,snr)
- elseif setup == "toeplitz"
- A, y, λ, M = data_toeplitz(k,m,n,snr)
- else
- error("Parameter `setup` must be either \"gaussian\" or \"toeplitz\".")
- end
- result_mip = solve_mip(A,y,λ,M,mipparams=MipParams(maxtime=maxtime))
- result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false))
- result_bnbs = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false,screening=true))
- catch e
- print_and_log("Precompilation failed with the following error\n",logfile)
- print_and_log("$e\n",logfile)
- print_and_log("Stopping experiment\n",logfile)
- rethrow(e)
- end
- print_and_log("\n",logfile)
- return nothing
-end
-
-logfile_name = string(
- @__DIR__,
- "/",
- Dates.format(Dates.now(),"yyyy-mm-dd-HH:MM:SS"),
- ".txt"
- )
-logfile = open(logfile_name, "w")
-
-print_and_log("--------------------------------\n",logfile)
-print_and_log("ICASSP experiment - Performances\n\n",logfile)
-print_and_log("Configuration\n",logfile)
-print_and_log(" setup : $(setup)\n",logfile)
-print_and_log(" k : $(k)\n",logfile)
-print_and_log(" m, n : $(m), $(n)\n",logfile)
-print_and_log(" snr : $(snr)\n",logfile)
-print_and_log(" maxtime : $(maxtime)\n",logfile)
-print_and_log(" maxiter : $(maxiter)\n",logfile)
-print_and_log(" tol : $(tol)\n",logfile)
-print_and_log(" repeats : $(repeats)\n",logfile)
-print_and_log("\n",logfile)
-
-fail_mip = Vector{Union{Bool,Missing}}(missing,repeats)
-fail_bnb = Vector{Union{Bool,Missing}}(missing,repeats)
-fail_bnbs = Vector{Union{Bool,Missing}}(missing,repeats)
-node_mip = Vector{Union{Int64,Missing}}(missing,repeats)
-node_bnb = Vector{Union{Int64,Missing}}(missing,repeats)
-node_bnbs = Vector{Union{Int64,Missing}}(missing,repeats)
-time_mip = Vector{Union{Float64,Missing}}(missing,repeats)
-time_bnb = Vector{Union{Float64,Missing}}(missing,repeats)
-time_bnbs = Vector{Union{Float64,Missing}}(missing,repeats)
-
-precompilation(setup,k,m,n,snr,maxtime,maxiter,tol)
-
-for i in 1:repeats
- print_and_log("Repeat $i/$repeats\n",logfile)
- try
- # Generate data
- if setup == "gaussian"
- A, y, λ, M = data_gaussian(k,m,n,snr)
- elseif setup == "toeplitz"
- A, y, λ, M = data_toeplitz(k,m,n,snr)
- end
-
- # Run the tree methods
- result_mip = solve_mip(A,y,λ,M,mipparams=MipParams(maxtime=maxtime))
- result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false))
- result_bnbs = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false,screening=true))
-
- # Gather results
- if result_mip.termination_status == MOI.OPTIMAL
- fail_mip[i] = false
- node_mip[i] = result_mip.node_count
- time_mip[i] = result_mip.solve_time
- else
- fail_mip[i] = true
- end
- if result_bnb.termination_status == MOI.OPTIMAL
- fail_bnb[i] = false
- node_bnb[i] = result_bnb.node_count
- time_bnb[i] = result_bnb.solve_time
- else
- fail_bnb[i] = true
- end
- if result_bnbs.termination_status == MOI.OPTIMAL
- fail_bnbs[i] = false
- node_bnbs[i] = result_bnbs.node_count
- time_bnbs[i] = result_bnbs.solve_time
- else
- fail_bnbs[i] = true
- end
-
- # Log results
- print_and_log(" Fails\n",logfile)
- print_and_log(" MIP : $(fail_mip[i])\n",logfile)
- print_and_log(" BnB : $(fail_bnb[i])\n",logfile)
- print_and_log(" SBnB : $(fail_bnbs[i])\n",logfile)
- print_and_log(" Node count\n",logfile)
- print_and_log(" MIP : $(node_mip[i])\n",logfile)
- print_and_log(" BnB : $(node_bnb[i])\n",logfile)
- print_and_log(" SBnB : $(node_bnbs[i])\n",logfile)
- print_and_log(" Solve time\n",logfile)
- print_and_log(" MIP : $(round(time_mip[i],digits=3)) seconds\n",logfile)
- print_and_log(" BnB : $(round(time_bnb[i],digits=3)) seconds\n",logfile)
- print_and_log(" SBnB : $(round(time_bnbs[i],digits=3)) seconds\n",logfile)
-
- catch e
- if e == InterruptException()
- close(logfile)
- rethrow(e)
- end
- print_and_log("warning : Repeat $i failed with the following error\n",logfile)
- print_and_log("$e\n",logfile)
- print_and_log("Skipping this repeat\n",logfile)
- end
-end
-
-
-################################### Results ####################################
-
-print_and_log("\n",logfile)
-print_and_log("Results\n",logfile)
-print_and_log(" Fails (mean over $repeats repeats)\n",logfile)
-print_and_log(" MIP : $(100 * mean(skipmissing(fail_mip)))%\n",logfile)
-print_and_log(" BnB : $(100 * mean(skipmissing(fail_bnb)))%\n",logfile)
-print_and_log(" SBnB : $(100 * mean(skipmissing(fail_bnbs)))%\n",logfile)
-print_and_log(" Node count (mean over $repeats repeats)\n",logfile)
-print_and_log(" MIP : $(mean(skipmissing(node_mip)))\n",logfile)
-print_and_log(" BnB : $(mean(skipmissing(node_bnb)))\n",logfile)
-print_and_log(" SBnB : $(mean(skipmissing(node_bnbs)))\n",logfile)
-print_and_log(" Solve time (mean over $repeats repeats)\n",logfile)
-print_and_log(" MIP : $(round(mean(skipmissing(time_mip)),digits=3)) seconds\n",logfile)
-print_and_log(" BnB : $(round(mean(skipmissing(time_bnb)),digits=3)) seconds\n",logfile)
-print_and_log(" SBnB : $(round(mean(skipmissing(time_bnbs)),digits=3)) seconds\n",logfile)
-print_and_log("--------------------------------",logfile)
-close(logfile)
diff --git a/src/BnbScreening.jl b/src/BnbScreening.jl
deleted file mode 100644
index 9ebc7e51b4e4e7657d325a9b8774f41287606f78..0000000000000000000000000000000000000000
--- a/src/BnbScreening.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-# Standard libraries
-using Dates
-using Printf
-using Random
-
-# External libraries
-using CPLEX
-using Distributions
-using JuMP
-using LinearAlgebra
-
-version() = "v1.0"
-author() = "Theo Guyard"
-contact() = "theo.guyard@insa-rennes.fr"
-license() = "AGPL 3.0 License (https://www.gnu.org/licenses/agpl-3.0.en.html)"
-
-# Module files
-include("types.jl")
-include("activeset.jl")
-include("bnb.jl")
-include("bvls.jl")
-include("data.jl")
-include("displays.jl")
-include("mips.jl")
-include("screen.jl")
diff --git a/src/activeset.jl b/src/activeset.jl
deleted file mode 100644
index 5365243c4520bbfc1ae40ca148b1df33ed1c48c6..0000000000000000000000000000000000000000
--- a/src/activeset.jl
+++ /dev/null
@@ -1,520 +0,0 @@
-tolsign(x::Float64, tol::Float64) = sign(x) * (abs(x) > tol)
-
-function symbi_app(Ainv::Matrix{Float64}, a::Vector{Float64})
- Ainv_updated = Matrix{Float64}(undef,length(a),length(a))
- @views begin
- d = a[1:end-1]
- g = a[end]
- dw = BLAS.symv('U',Ainv,-d)
- gw = 1. / (g + dw' * d)
- gw_dw = gw .* dw
- BLAS.syr!('U',gw,dw,Ainv)
- Ainv_updated[1:end-1,1:end-1] = Ainv
- Ainv_updated[1:end-1,end] = gw_dw
- Ainv_updated[end,end] = gw
- end
- return Ainv_updated
-end
-
-function symbi_rem(Ainv::Matrix{Float64}, i::Int)
- n = size(Ainv)[1]
- Ainv_updated = Matrix{Float64}(undef,n-1,n-1)
- @views begin
- q = -1. / Ainv[i,i]
- c1 = 1:i-1
- r1 = Ainv[c1,i]
- c2 = i+1:n
- r2 = Ainv[i,c2]
- c2red = i:n-1
- Ainv_updated[c1,c1] = BLAS.syr!('U',q,r1,Ainv[c1,c1])
- Ainv_updated[c1,c2red] = BLAS.ger!(q,r1,r2,Ainv[c1,c2])
- Ainv_updated[c2red,c2red] = BLAS.syr!('U',q,r2,Ainv[c2,c2])
- end
- return Ainv_updated
-end
-
-function symbi_ro!(Ainv::Matrix{Float64}, r::Vector{Float64}, g::Float64)
- Ainv_r = BLAS.symv('U',Ainv,r)
- c = -g / (1. + g * (r' * Ainv_r))
- BLAS.syr!('U',c,Ainv_r,Ainv)
- return nothing
-end
-
-function swap_Sbarin_to_Sbar0!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- i::Int
- )
- iSbarin = findfirst(j->j==i,node.support.Sbarin)
- node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
- push!(node.support.Sbar0,i)
- deleteat!(node.support.Sbarin,iSbarin)
- deleteat!(asviews.xSbarin,iSbarin)
- deleteat!(asviews.sSbarin,iSbarin)
- deleteat!(asviews.xnewSbarin,iSbarin)
- asviews.ASbarin = A[:,node.support.Sbarin]
- return nothing
-end
-
-function swap_Sbarin_to_Sbarbox!(node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- i::Int
- )
- iSbarin = findfirst(j->j==i,node.support.Sbarin)
- node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
- push!(node.support.Sbarbox,i)
- push!(asviews.xSbarbox,popat!(asviews.xSbarin,iSbarin))
- push!(asviews.sSbarbox,popat!(asviews.sSbarin,iSbarin))
- deleteat!(asviews.xnewSbarin,iSbarin)
- deleteat!(node.support.Sbarin,iSbarin)
- asviews.ASbarin = A[:,node.support.Sbarin]
- asviews.ASbarbox = A[:,node.support.Sbarbox]
- return nothing
-end
-
-function swap_S1in_to_S1box!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- i::Int
- )
- iS1in = findfirst(j->j==i,node.support.S1in)
- push!(node.support.S1box,i)
- push!(asviews.xS1box,popat!(asviews.xS1in,iS1in))
- push!(asviews.sS1box,popat!(asviews.sS1in,iS1in))
- deleteat!(asviews.xnewS1in,iS1in)
- deleteat!(node.support.S1in,iS1in)
- asviews.AS1in = A[:,node.support.S1in]
- asviews.AS1box = A[:,node.support.S1box]
-
- @views begin
- a = A[:,i]
- d = vcat(node.GS1in_inv[1:(iS1in-1),iS1in], node.GS1in_inv[iS1in,(iS1in+1):end])
- e = node.GS1in_inv[iS1in,iS1in]
- Ad = asviews.AS1in * d
- r1 = (Ad./e) + a
- BLAS.syr!('U',e,r1,B)
- node.GS1in_inv = symbi_rem(node.GS1in_inv,iS1in)
- symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1, e)
- end
-
- return nothing
-end
-
-function swap_Sbar0_to_Sbarin!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- i::Int
- )
- @views begin
- a = A[:,i]
- Ba = BLAS.symv('U',B,a)
- ABa = asviews.ASbarin' * Ba
- aBa = a' * Ba
- v = vcat(ABa,aBa)
- node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
- end
- push!(node.support.Sbarin,i)
- filter!(j->j!=i,node.support.Sbar0)
- push!(asviews.xSbarin, 0.)
- push!(asviews.sSbarin, 0.)
- push!(asviews.xnewSbarin, 0.)
- asviews.ASbarin = A[:,node.support.Sbarin]
- return nothing
-end
-
-function swap_Sbarbox_to_Sbarin!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- i::Int
- )
- iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
- @views begin
- a = A[:,i]
- Ba = BLAS.symv('U',B,a)
- ABa = asviews.ASbarin' * Ba
- aBa = a' * Ba
- v = vcat(ABa,aBa)
- node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
- end
- push!(node.support.Sbarin, popat!(node.support.Sbarbox,iSbarbox))
- push!(asviews.xSbarin, popat!(asviews.xSbarbox,iSbarbox))
- push!(asviews.sSbarin, popat!(asviews.sSbarbox,iSbarbox))
- push!(asviews.xnewSbarin, asviews.xSbarin[end])
- asviews.ASbarin = A[:,node.support.Sbarin]
- asviews.ASbarbox = A[:,node.support.Sbarbox]
- return nothing
-end
-
-function swap_S1box_to_Sbarin!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- i::Int
- )
- iS1box = findfirst(j->j==i,node.support.S1box)
- @views begin
- a = A[:,i]
- Ba = BLAS.symv('U',B,a)
- ABa = asviews.ASbarin' * Ba
- aBa = a' * Ba
- v = vcat(ABa,aBa)
- node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
- end
- push!(node.support.Sbarin, popat!(node.support.S1box,iS1box))
- push!(asviews.xSbarin, popat!(asviews.xS1box,iS1box))
- push!(asviews.sSbarin, popat!(asviews.sS1box,iS1box))
- push!(asviews.xnewSbarin, asviews.xSbarin[end])
- asviews.ASbarin = A[:,node.support.Sbarin]
- asviews.AS1box = A[:,node.support.S1box]
- return nothing
-end
-
-function swap_Sbar_to_S0!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- i::Int
- )
- iSbarin = findfirst(j->j==i,node.support.Sbarin)
- if !isa(iSbarin,Nothing)
- node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
- deleteat!(node.support.Sbarin,iSbarin)
- deleteat!(asviews.xSbarin,iSbarin)
- deleteat!(asviews.sSbarin,iSbarin)
- deleteat!(asviews.xnewSbarin,iSbarin)
- asviews.ASbarin = A[:,node.support.Sbarin]
- else
- iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
- if !isa(iSbarbox,Nothing)
- deleteat!(node.support.Sbarbox,iSbarbox)
- deleteat!(asviews.xSbarbox,iSbarbox)
- deleteat!(asviews.sSbarbox,iSbarbox)
- asviews.ASbarbox = A[:,node.support.Sbarbox]
- else
- filter!(j->j!=i,node.support.Sbar0)
- end
- end
- return nothing
-end
-
-function swap_Sbar_to_S1!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- i::Int
- )
- iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
- if !isa(iSbarbox,Nothing)
- push!(node.support.S1box, popat!(node.support.Sbarbox,iSbarbox))
- push!(asviews.xS1box, popat!(asviews.xSbarbox,iSbarbox))
- push!(asviews.sS1box, popat!(asviews.sSbarbox,iSbarbox))
- asviews.ASbarbox = A[:,node.support.Sbarbox]
- asviews.AS1box = A[:,node.support.S1box]
- else
- @views begin
- asviews.AS1in = A[:,node.support.S1in]
- a = A[:,i]
- Ata = asviews.AS1in' * a
- ata = a' * a
- d = - BLAS.symv('U',node.GS1in_inv,Ata)
- g = -1. / (ata + d'*Ata)
- r1 = asviews.AS1in * d + a
- v = vcat(Ata,ata)
- node.GS1in_inv = symbi_app(node.GS1in_inv,v)
- BLAS.syr!('U',g,r1,B)
- end
- iSbarin = findfirst(j->j==i,node.support.Sbarin)
- if !isa(iSbarin,Nothing)
- node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
- push!(node.support.S1in, popat!(node.support.Sbarin,iSbarin))
- push!(asviews.xS1in, popat!(asviews.xSbarin,iSbarin))
- push!(asviews.sS1in, popat!(asviews.sSbarin,iSbarin))
- push!(asviews.xnewS1in, popat!(asviews.xnewSbarin,iSbarin))
- asviews.ASbarin = A[:,node.support.Sbarin]
- asviews.AS1in = A[:,node.support.S1in]
- symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1,g)
- else
- iSbar0 = findfirst(j->j==i,node.support.Sbar0)
- push!(node.support.S1in, popat!(node.support.Sbar0,iSbar0))
- push!(asviews.xS1in, 0.)
- push!(asviews.sS1in, 0.)
- push!(asviews.xnewS1in, 0.)
- asviews.AS1in = A[:,node.support.S1in]
- symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1,g)
- end
- end
- return nothing
-end
-
-function solve_kkt_system!(
- node::Node,
- asviews::ActivesetViews,
- B::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64
- )
- Axbox = asviews.ASbarbox * asviews.xSbarbox + asviews.AS1box * asviews.xS1box
- r = y - Axbox
- rhsSbarin = asviews.ASbarin' * BLAS.symv('U',B,r) - (λ/M) .* asviews.sSbarin
- asviews.xnewSbarin = BLAS.symv('U',node.GSbarin_inv,rhsSbarin)
- rhsS1in = asviews.AS1in' * (r - asviews.ASbarin * asviews.xnewSbarin)
- asviews.xnewS1in = BLAS.symv('U',node.GS1in_inv,rhsS1in)
- return Axbox
-end
-
-function discrete_line_search!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- Axbox::Vector{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams,
- )
-
- # TODO : optimize
-
- xdiffSbarin = asviews.xnewSbarin - asviews.xSbarin
- xdiffS1in = asviews.xnewS1in - asviews.xS1in
- Ax = asviews.ASbarin * asviews.xSbarin + asviews.AS1in * asviews.xS1in + Axbox
- Axdiff = asviews.ASbarin * xdiffSbarin + asviews.AS1in * xdiffS1in
- u = y - Ax
- λ_M = λ / M
- λ_S1 = λ * length(node.branch.S1)
-
- a = Axdiff' * Axdiff
- b_base = u' * Axdiff
-
- if a == 0.
- Atu = A' * u
- lb = 0.5*u'*u + λ_M*asviews.xSbarin'*asviews.sSbarin + λ_S1
- return u, Atu, lb
- end
-
- tflip = []
- tmax = 1.
- for i in 1:length(node.support.Sbarin)
- if xdiffSbarin[i] < 0.
- tmax = min(tmax, -(M+asviews.xSbarin[i])/xdiffSbarin[i])
- else
- tmax = min(tmax, (M-asviews.xSbarin[i])/xdiffSbarin[i])
- end
- push!(tflip, -asviews.xSbarin[i]/xdiffSbarin[i])
- end
- for i in 1:length(node.support.S1in)
- if xdiffS1in[i] < 0.
- tmax = min(tmax, -(M+asviews.xS1in[i])/xdiffS1in[i])
- else
- tmax = min(tmax, (M-asviews.xS1in[i])/xdiffS1in[i])
- end
- push!(tflip, -asviews.xS1in[i]/xdiffS1in[i])
- end
-
- filter!(ti -> 0. < ti < tmax, sort!(unique!(tflip)))
- pushfirst!(tflip,0.)
- push!(tflip,tmax)
-
- topt = 0.
- uopt = u
- vopt = Inf
-
- for i in 1:(length(tflip)-1)
-
- tl = tflip[i]
- tu = tflip[i+1]
- tmed = (tl+tu)/2.
-
- thetaSbarin = sign.(asviews.xSbarin + tmed .* xdiffSbarin)
-
- b = b_base - λ_M * thetaSbarin' * xdiffSbarin
-
- tinf = max(tl,min(tu,b/a))
- uinf = u - tinf .* Axdiff
- vinf = 0.5 * norm(uinf,2)^2 + λ_M * norm(asviews.xSbarin + tinf .* xdiffSbarin,1)
-
- if vinf < vopt
- topt = tinf
- uopt = uinf
- vopt = vinf
- end
- end
-
- asviews.xSbarin += topt .* xdiffSbarin
- asviews.xS1in += topt .* xdiffS1in
- asviews.sSbarin = tolsign.(asviews.xSbarin,bnbparams.tol)
- asviews.sS1in = tolsign.(asviews.xS1in,bnbparams.tol)
- Atuopt = A' * uopt
- lb = vopt + (λ_S1 + λ * sum(node.support.Sbarbox))
-
- return uopt, Atuopt, lb
-end
-
-function detect_blocking_behaviour(swapped::Vector{Int}, imax::Int)
- return (length(swapped)==1) & (imax in swapped)
-end
-
-function active_optimality(
- asviews::ActivesetViews,
- M::Float64,
- bnbparams::BnbParams
- )
- cond = (
- all(asviews.sSbarin .== tolsign.(asviews.xnewSbarin,bnbparams.tol)) &
- all(asviews.sS1in .== tolsign.(asviews.xnewS1in,bnbparams.tol)) &
- (norm(asviews.xnewSbarin,Inf) - M <= bnbparams.tol) &
- (norm(asviews.xnewS1in,Inf) - M <= bnbparams.tol)
- )
- return cond
-end
-
-function inactive_optimality(
- node::Node,
- asviews::ActivesetViews,
- Atu::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams
- )
- vmax = 0.
- imax = nothing
- smax = ""
- λ_M = λ / M
- λ_M_m_tol = λ_M - bnbparams.tol
- λ_M_p_tol = λ_M + bnbparams.tol
- m_tol = -bnbparams.tol
- for i in node.support.Sbar0
- v = abs(Atu[i]) - λ_M_p_tol
- if v > vmax
- vmax = v
- imax = i
- smax = "Sbar0"
- end
- end
- for (iSbarbox,i) in enumerate(node.support.Sbarbox)
- v = λ_M_m_tol - Atu[i] * asviews.sSbarbox[iSbarbox]
- if v > vmax
- vmax = v
- imax = i
- smax = "Sbarbox"
- end
- end
- for (iS1box,i) in enumerate(node.support.S1box)
- v = m_tol - Atu[i] * asviews.sS1box[iS1box]
- if v > vmax
- vmax = v
- imax = i
- smax = "S1box"
- end
- end
- return imax, smax
-end
-
-function new_tentative_sign!(
- node::Node,
- asviews::ActivesetViews,
- imax::Int,
- smax::String,
- A::Matrix{Float64},
- B::Matrix{Float64}
- )
- if smax == "Sbar0"
- swap_Sbar0_to_Sbarin!(node,asviews,A,B,imax)
- elseif smax == "Sbarbox"
- swap_Sbarbox_to_Sbarin!(node,asviews,A,B,imax)
- elseif smax == "S1box"
- swap_S1box_to_Sbarin!(node,asviews,A,B,imax)
- end
- return nothing
-end
-
-function update_sets!(
- node::Node,
- asviews::ActivesetViews,
- A::Matrix{Float64},
- B::Matrix{Float64},
- M::Float64,
- bnbparams::BnbParams
- )
-
- swapped = Vector{Int}()
- for (iS1in,i) in enumerate(node.support.S1in)
- absxi = abs(asviews.xS1in[iS1in])
- if abs(M - absxi) <= bnbparams.tol
- swap_S1in_to_S1box!(node,asviews,A,B,i)
- push!(swapped,i)
- end
- end
- for (iSbarin,i) in enumerate(node.support.Sbarin)
- absxi = abs(asviews.xSbarin[iSbarin])
- if absxi <= bnbparams.tol
- swap_Sbarin_to_Sbar0!(node,asviews,A,i)
- push!(swapped,i)
- elseif abs(M - absxi) <= bnbparams.tol
- swap_Sbarin_to_Sbarbox!(node,asviews,A,i)
- push!(swapped,i)
- end
- end
- return swapped
-end
-
-function solve_relax_activeset!(
- tree::Tree,
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams,
- )
-
- asviews = ActivesetViews(node,A,bnbparams)
- B = I(size(A)[1]) - BLAS.symm('R','U',node.GS1in_inv,asviews.AS1in) * asviews.AS1in'
- node.branch_on[2] == 0 && swap_Sbar_to_S0!(node,asviews,A,node.branch_on[1])
- node.branch_on[2] == 1 && swap_Sbar_to_S1!(node,asviews,A,B,node.branch_on[1])
-
- it = 0
- lb = Inf
-
- while true
- it += 1
- it >= bnbparams.maxiter && error("activetset : iteration limit")
-
- Axbox = solve_kkt_system!(node,asviews,B,y,λ,M)
- u, Atu, lb = discrete_line_search!(node,asviews,A,Axbox,y,λ,M,bnbparams)
- swapped = update_sets!(node,asviews,A,B,M,bnbparams)
-
- if bnbparams.screening
- prune = screen!(tree,node,asviews,u,Atu,B,A,y,λ,M)
- prune && break
- end
-
- if active_optimality(asviews,M,bnbparams)
- imax, smax = inactive_optimality(node,asviews,Atu,λ,M,bnbparams)
- isa(imax,Nothing) && break
- new_tentative_sign!(node,asviews,imax,smax,A,B)
- detect_blocking_behaviour(swapped,imax) && break
- end
- end
-
- node.lb = lb
- node.x[node.support.Sbar0] .= 0.
- node.x[node.support.Sbarin] = asviews.xSbarin
- node.x[node.support.Sbarbox] = asviews.xSbarbox
- node.x[node.support.S1in] = asviews.xS1in
- node.x[node.support.S1box] = asviews.xS1box
-
- return nothing
-end
-
diff --git a/src/bnb.jl b/src/bnb.jl
deleted file mode 100644
index ec8819f20831f879dcf801adbb6fbdffdea8833a..0000000000000000000000000000000000000000
--- a/src/bnb.jl
+++ /dev/null
@@ -1,191 +0,0 @@
-function best_lb(tree::Tree)
- isempty(tree.active) && return 0.
- return maximum([node.lb for node in tree.active])
-end
-
-function branch!(tree::Tree, node::Node)
- prune(tree,node) && return nothing
- j = branching_index(node)
- isa(j,Nothing) && return nothing
- node_j0 = Node(node,j,0)
- node_j1 = Node(node,j,1)
- push!(tree.queue,node_j0)
- push!(tree.queue,node_j1)
-end
-
-function branching_index(node::Node)
- isempty(node.branch.Sbar) && return nothing
- return node.branch.Sbar[argmax(abs.(node.x[node.branch.Sbar]))]
-end
-
-function depth(node::Node)
- depth = 0
- while !isa(node,Nothing)
- depth += 1
- node = node.parent
- end
- return depth
-end
-
-function elapsed_time(tree::Tree)
- return (Dates.now() - tree.start_time).value / 1000.
-end
-
-function fixto!(node::Node, j::Int, jval::Int)
- jpos = findfirst(i->i==j,node.branch.Sbar)
- isa(jpos,Nothing) && error("Branching index is already fixed")
- if jval == 0
- push!(node.branch.S0,popat!(node.branch.Sbar,jpos))
- node.x[j] = 0.
- elseif jval == 1
- push!(node.branch.S1,popat!(node.branch.Sbar,jpos))
- end
- return nothing
-end
-
-function global_gap(tree::Tree)
- return (tree.ub - best_lb(tree)) / tree.ub
-end
-
-function has_perfect_relaxation(node::Node, M::Float64, bnbparams::BnbParams)
- return !(any(bnbparams.tol .< abs.(node.x[node.branch.Sbar]) .< M - bnbparams.tol))
-end
-
-function lower_bound!(
- tree::Tree,
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams
- )
- if bnbparams.lb_method == ACTIVESET
- solve_relax_activeset!(tree,node,A,y,λ,M,bnbparams)
- elseif bnbparams.lb_method == LBMIP
- solve_relax_mip!(node,A,y,λ,M)
- else
- error("Paramter `bnbparams.lb_method` not understood.")
- end
- return nothing
-end
-
-function next_node!(tree::Tree)
- return pop!(tree.queue)
-end
-
-function processed!(tree::Tree, node::Node)
- if !prune(tree,node)
- push!(tree.active,node)
- end
- tree.nexpl += 1
- return nothing
-end
-
-function prune(tree::Tree, node::Node)
- return (tree.ub < node.lb)
-end
-
-function terminated!(tree::Tree, bnbparams::BnbParams)
- if isempty(tree.queue)
- tree.termination_status = MOI.OPTIMAL
- elseif elapsed_time(tree) > bnbparams.maxtime
- tree.termination_status = MOI.TIME_LIMIT
- end
- return (tree.termination_status != MOI.OPTIMIZE_NOT_CALLED)
-end
-
-function update_upper_bound!(tree::Tree, ub::Float64, inc::Vector{Float64})
- if ub < tree.ub
- tree.ub = ub
- tree.inc = inc
- return true
- end
- return false
-end
-
-function update_active_nodes!(tree::Tree)
- filter!(node->!prune(tree,node),tree.active)
- return nothing
-end
-
-function upper_bound!(
- tree::Tree,
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams
- )
-
- perfect = false
-
- if isempty(node.branch.S1)
- ub = 0.5*y'*y
- inc = zeros(Float64,size(A)[2])
- elseif has_perfect_relaxation(node,M,bnbparams)
- ub = node.lb
- inc = node.x
- perfect = true
- else
- if bnbparams.ub_method == BVLS
- ub, inc = solve_heuristic_bvls(node,A,y,λ,M,bnbparams)
- elseif bnbparams.ub_method == UBMIP
- ub, inc = solve_heuristic_mip(node,A,y,λ,M)
- else
- error("Paramter `bnbparams.ub_method` not understood.")
- end
- end
- updated = update_upper_bound!(tree,ub,inc)
- updated && update_active_nodes!(tree)
- return perfect
-end
-
-# ============================================================================ #
-
-"""
- solve_bnb(
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64;
- bnbparams::BnbParams=BnbParams(),
- )
-
-Solve the L0-penalized least-square problem
-
- min_x (1/2) norm(y-Ax,2)^2 + λ norm(x,0)
- st. norm(x,Inf) <= M
-
-using a Branch-and-Bound algorithm. This algorithm can use branch-screening
-tests to enhance its efficiency. Returns a [`@BnbResults`](@ref) struct
-containing the optimization results and statistics. See [`@BnbParams`](@ref) for
-informations about additional method parameters.
-"""
-function solve_bnb(
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64;
- bnbparams::BnbParams=BnbParams(),
- )
-
- validate_data(A,y,λ,M)
- tree = Tree(A,y)
- bnbparams.verbosity && display_head()
-
- while !terminated!(tree,bnbparams)
- node = next_node!(tree)
- lower_bound!(tree,node,A,y,λ,M,bnbparams)
- perfect = upper_bound!(tree,node,A,y,λ,M,bnbparams)
- processed!(tree,node)
- perfect || branch!(tree,node)
- bnbparams.verbosity && display_node(tree,node)
- end
-
- results = BnbResults(tree)
- bnbparams.verbosity && display_tail(results)
-
- return results
-end
\ No newline at end of file
diff --git a/src/bvls.jl b/src/bvls.jl
deleted file mode 100644
index b7a3aaad14758555c36746ac0f44cbf40e4c7091..0000000000000000000000000000000000000000
--- a/src/bvls.jl
+++ /dev/null
@@ -1,132 +0,0 @@
-function compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector{Int})
- length(g) == 0 && return 0.
- @views begin
- g_kkt = g .* on_bound
- free_set = (on_bound .== 0)
- g_kkt[free_set] = abs.(g[free_set])
- end
- return maximum(g_kkt)
-end
-
-function solve_heuristic_bvls(
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- bnbparams::BnbParams,
- )
-
- AS1 = A[:,node.branch.S1]
- n = size(AS1)[2]
- lb = fill(-M,(n,))
- ub = fill(M,(n,))
-
- x = zeros(Float64,n)
- on_bound = zeros(Int,n)
- free_set = (on_bound .== 0)
- active_set = .!free_set
- free_set = findall(free_set)
-
- r = AS1 * x - y
- cost = r' * r
- g = AS1' * r
- cost_change = Inf
-
- # ----- Initialization loop ----- #
-
- while !isempty(free_set)
-
- AS1_free = AS1[:, free_set]
- y_free = y - AS1 * (x .* active_set)
- z = AS1_free \ y_free
- lbv = (z .< lb[free_set])
- ubv = (z .> ub[free_set])
- v = lbv .| ubv
-
- if any(lbv)
- ind = free_set[lbv]
- x[ind] = lb[ind]
- active_set[ind] .= true
- on_bound[ind] .= -1
- end
- if any(ubv)
- ind = free_set[ubv]
- x[ind] = ub[ind]
- active_set[ind] .= true
- on_bound[ind] .= 1
- end
-
- ind = free_set[.!v]
- x[ind] = z[.!v]
- r = AS1 * x - y
- g = AS1' * r
- any(v) || break
- free_set = free_set[.!v]
- end
-
- # ----- Main loop ----- #
-
- optimality = compute_kkt_optimality(g,on_bound)
-
- while true
-
- (optimality < bnbparams.tol) && break
- move_to_free = argmax(g .* on_bound)
- on_bound[move_to_free] = 0
-
- while true
-
- free_set = (on_bound .== 0)
- active_set = .!free_set
- free_set = findall(free_set)
- x_free = x[free_set]
- lb_free = lb[free_set]
- ub_free = ub[free_set]
- AS1_free = AS1[:,free_set]
- y_free = y - AS1 * (x .* active_set)
- z = AS1_free \ y_free
- lbv = findall(z .< lb_free)
- ubv = findall(z .> ub_free)
- v = vcat(lbv,ubv)
-
- if !isempty(v)
-
- alphas = vcat(
- lb_free[lbv] - x_free[lbv],
- ub_free[ubv] - x_free[ubv]
- ) ./ (z[v] - x_free[v])
- i = argmin(alphas)
- i_free = v[i]
- alpha = alphas[i]
- x_free *= 1. - alpha
- x_free += alpha * z
- x[free_set] = x_free
-
- if i < length(lbv)
- on_bound[free_set[i_free]] = -1
- else
- on_bound[free_set[i_free]] = 1
- end
- else
- x_free = z
- x[free_set] = x_free
- break
- end
- end
-
- r = AS1 * x - y
- cost_new = 0.5 * (r' * r)
- cost_change = cost - cost_new
-
- (cost_change < bnbparams.tol * cost) && break
-
- cost = cost_new
- g = AS1' * r
- optimality = compute_kkt_optimality(g,on_bound)
- end
-
- obj = 0.5 * norm(r,2)^2 + λ * n
-
- return obj, x
-end
\ No newline at end of file
diff --git a/src/data.jl b/src/data.jl
deleted file mode 100644
index f3d730788600ebec214d5f62c60d8ad3880225fe..0000000000000000000000000000000000000000
--- a/src/data.jl
+++ /dev/null
@@ -1,89 +0,0 @@
-function validate_data(
- A::Array{Float64,2},
- y::Array{Float64,1},
- λ::Float64,
- M::Float64
- )
- size(A)[1] == size(y)[1] || throw(ArgumentError("Dimension missmatch"))
- λ > 0 || throw(ArgumentError("λ must be positive"))
- M > 0 || throw(ArgumentError("M must be positive"))
-end
-
-# ============================================================================ #
-
-"""
- data_gaussian(
- k::Int=5,
- m::Int=500,
- n::Int=1000,
- snr::Float64=10.
- )
-
-Generate a Gaussian setup.
-"""
-function data_gaussian(
- k::Int=5,
- m::Int=500,
- n::Int=1000,
- snr::Float64=10.
- )
-
- x = zeros(n)
- a = rand(Normal(0,1),k)
- a += sign.(a)
- x[randperm(n)[1:k]] = a
-
- A = rand(Normal(0,1),m,n)
- for i in 1:n
- A[:,i] /= norm(A[:,i],2)
- end
-
- σ = norm(A*x,2) / sqrt(m * 10^(snr/10))
- ϵ = rand(Normal(0,σ),m)
- y = A * x .+ ϵ
- λ = 2 * σ * log(n/k - 1)
- M = 1.5 * norm(A'*y,Inf)
-
- return A, y, λ, M
-end
-
-"""
- data_gaussian(
- k::Int=5,
- m::Int=500,
- n::Int=1000,
- snr::Float64=10.
- )
-
-Generate a Toeplitz setup.
-"""
-function data_toeplitz(
- k::Int=5,
- m::Int=500,
- n::Int=300,
- snr::Float64=10.
- )
-
- x = zeros(n)
- a = rand(Normal(0,1),k)
- a += sign.(a)
- x[randperm(n)[1:k]] = a
-
- tmax = n/2
- t = collect(1:m)
- A = zeros(m,n)
- offset = (n-m)/2
- for i in 1:n
- A[:,i] = @. sinc((offset + (t - i)) * (tmax/m))
- A[:,i] /= norm(A[:,i],2)
- end
- A += rand(m,n)./10000
-
- σ = norm(A*x,2) / sqrt(m * 10^(snr/10))
- ϵ = rand(Normal(0,σ),m)
- y = A * x .+ ϵ
- λ = 2 * σ * log(n/k - 1)
- M = 1.5 * norm(A'*y,Inf)
-
- return A, y, λ, M
-end
\ No newline at end of file
diff --git a/src/displays.jl b/src/displays.jl
deleted file mode 100644
index 391dde8faacd9380d404ba754c5bc0467358bca8..0000000000000000000000000000000000000000
--- a/src/displays.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-WIDTH = 67
-
-function display_head()
- println(repeat("*",WIDTH))
- println("B&B-screening solver for the L0-penalized least-square problem")
- println("Version : $(version())")
- println("Author : $(author())")
- println("Contact : $(contact())")
- println("License : $(license())")
- println()
- println(repeat("-",WIDTH))
- println("| Node | Bounds | Branch |")
- println("| nexpl depth queue| lb ub gap| S0 S1 Sbar|")
- println(repeat("-",WIDTH))
- return nothing
-end
-
-function display_node(tree::Tree, node::Node)
-
- @printf(
- "|%7d%7d%7d|%7.3f%7.3f%7.0e|%7d%7d%7d|",
- tree.nexpl,
- depth(node),
- length(tree.queue),
- best_lb(tree),
- tree.ub,
- global_gap(tree),
- length(node.branch.S0),
- length(node.branch.S1),
- length(node.branch.Sbar),
- )
- println()
- return nothing
-end
-
-function diplay_results(results::BnbResults)
- println("B&B results")
- println(" Status : $(results.termination_status)")
- println(" Solve time : $(results.solve_time)")
- println(" Nodes : $(results.node_count)")
- println(" Objective : $(results.objective_value)")
- return nothing
-end
-
-function display_tail(results::BnbResults)
- println(repeat("-",WIDTH))
- println()
- diplay_results(results)
- println()
- println(repeat("*",WIDTH))
- return nothing
-end
-
diff --git a/src/mips.jl b/src/mips.jl
deleted file mode 100644
index 7199214ade0ae0923b596283cfcbfd19eb6f5697..0000000000000000000000000000000000000000
--- a/src/mips.jl
+++ /dev/null
@@ -1,117 +0,0 @@
-function solve_relax_mip!(
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- )
-
- relax = Model(CPLEX.Optimizer)
- set_silent(relax)
-
- gram = A'*A
- Aty = A'*y
-
- n = size(A)[2]
- @variable(relax, x[1:n])
- @variable(relax, xabs[1:n])
- @objective(relax, Min,
- 0.5*(y'*y) - Aty'*x + (0.5*x)'*(gram*x) +
- (λ/M)*sum(xabs[node.branch.Sbar]) +
- λ*length(node.branch.S1)
- )
- @constraint(relax, box, -M .<= x .<= M)
- @constraint(relax, abs_pos, xabs .>= x)
- @constraint(relax, abs_neg, xabs .>= -x)
- @constraint(relax, setzero, x[node.branch.S0] .== 0)
- @constraint(relax, abs_setzero, xabs[node.branch.S0] .== 0)
-
- for i in 1:n
- set_start_value(x[i], node.x[i])
- set_start_value(xabs[i], abs(node.x[i]))
- end
-
- optimize!(relax)
-
- node.lb = objective_value(relax)
- node.x = value.(x)
-end
-
-function solve_heuristic_mip(
- node::Node,
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- )
-
- heur = Model(CPLEX.Optimizer)
- set_silent(heur)
-
- nS1 = length(node.branch.S1)
- AS1 = A[:,node.branch.S1]
- xS1 = node.x[node.branch.S1]
-
- @variable(heur, x[1:nS1])
- @objective(heur, Min, 0.5*(y-AS1*x)'*(y-AS1*x)+ λ*length(node.branch.S1))
- @constraint(heur, box, -M .<= x .<= M)
-
- for i in 1:nS1
- set_start_value(x[i], xS1[i])
- end
-
- optimize!(heur)
-
- inc = zeros(Float64,size(A)[2])
- inc[node.branch.S1] = value.(x)
-
- return objective_value(heur), inc
-end
-
-# ============================================================================ #
-
-"""
- solve_mip(
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64;
- mipparams::MipParams=MipParams(),
- )
-
-Solve the L0-penalized least-square problem
-
- min_x (1/2) norm(y-Ax,2)^2 + λ norm(x,0)
- st. norm(x,Inf) <= M
-
-using a MIP solver (CPLEX). Returns a [`@MipResults`](@ref) struct containing
-the optimization results and statistics. See [`@MipParams`](@ref) for
-informations about additional method parameters.
-"""
-function solve_mip(
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64;
- mipparams::MipParams=MipParams(),
- )
-
- problem = Model(CPLEX.Optimizer)
- set_optimizer_attribute(problem,"CPX_PARAM_TILIM",mipparams.maxtime)
- mipparams.silent && set_silent(problem)
-
- gram = A'*A
- Aty = A'*y
-
- n = size(A)[2]
- @variable(problem, x[1:n])
- @variable(problem, z[1:n], Bin)
- @objective(problem, Min, 0.5*(y'*y) - Aty'*x + (0.5*x)'*(gram*x) + λ*sum(z))
- @constraint(problem, box_lower, -M.*z .<= x)
- @constraint(problem, box_upper, x .<= M.*z)
-
- optimize!(problem)
- result = MipResults(problem)
-
- return result
-end
diff --git a/src/screen.jl b/src/screen.jl
deleted file mode 100644
index 569c17dda60efb119b9debce1edf826cbc48464d..0000000000000000000000000000000000000000
--- a/src/screen.jl
+++ /dev/null
@@ -1,52 +0,0 @@
-function screen!(
- tree::Tree,
- node::Node,
- asviews::ActivesetViews,
- u::Vector{Float64},
- Atu::Vector{Float64},
- B::Matrix{Float64},
- A::Matrix{Float64},
- y::Vector{Float64},
- λ::Float64,
- M::Float64,
- )
-
- @views begin
- pivotSbar = @. M * (abs(Atu[node.branch.Sbar]) - λ/M)
- pivotSbar0 = max.(pivotSbar,0)
- pivotSbar1 = max.(-pivotSbar,0)
- pivotS1 = @. M * (abs(Atu[node.branch.S1]) - λ/M)
- dobj = 0.5 * (y'y - norm(y-u,2)^2) - sum(pivotSbar0) - sum(pivotS1)
- end
-
- Sbar_to_S0 = Vector{Int}()
- Sbar_to_S1 = Vector{Int}()
- for (iSbar,i) in enumerate(node.branch.Sbar)
- if pivotSbar1[iSbar] > tree.ub - dobj
- push!(Sbar_to_S0,i)
- tree.scr0 += 1
- end
- if pivotSbar0[iSbar] > tree.ub - dobj
- push!(Sbar_to_S1,i)
- tree.scr1 += 1
- end
- end
-
- if !isempty(intersect(Sbar_to_S0,Sbar_to_S1))
- node.lb = Inf
- node.x = rand(size(A)[2])
- return true
- end
-
- for i in Sbar_to_S0
- fixto!(node,i,0)
- swap_Sbar_to_S0!(node,asviews,A,i)
- end
-
- for i in Sbar_to_S1
- fixto!(node,i,1)
- swap_Sbar_to_S1!(node,asviews,A,B,i)
- end
-
- return false
-end
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
deleted file mode 100644
index 16e48d4a6dbc4243770455188e8a4d647aeeadf3..0000000000000000000000000000000000000000
--- a/src/types.jl
+++ /dev/null
@@ -1,295 +0,0 @@
-# ========================================== #
-# Active-Set support configuration and views #
-# ========================================== #
-
-mutable struct Support
- Sbar0::Vector{Int} # Indices verifying 0 = xi with i in Sbar
- Sbarin::Vector{Int} # Indices verifying 0 < xi < M with i in Sbar
- Sbarbox::Vector{Int} # Indices verifying xi = M with i in Sbar
- S1in::Vector{Int} # Indices verifying xi < M with i in S1
- S1box::Vector{Int} # Indices verifying xi = M with i in S1
-end
-
-mutable struct ActivesetViews
- # Current iterate
- xSbarin::Vector{Float64}
- xSbarbox::Vector{Float64}
- xS1in::Vector{Float64}
- xS1box::Vector{Float64}
-
- # Sign of x
- sSbarin::Vector{Float64}
- sSbarbox::Vector{Float64}
- sS1in::Vector{Float64}
- sS1box::Vector{Float64}
-
- # New iterate
- xnewSbarin::Vector{Float64}
- xnewS1in::Vector{Float64}
-
- # Submatrices
- ASbarin::Matrix{Float64}
- ASbarbox::Matrix{Float64}
- AS1in::Matrix{Float64}
- AS1box::Matrix{Float64}
-end
-
-
-# =========================== #
-# Branch-and-Bound components #
-# =========================== #
-
-mutable struct Branch
- S0::Vector{Int} # Indices fixed to zero
- S1::Vector{Int} # Indices fixed to non-zero
- Sbar::Vector{Int} # Indices not fixed
-end
-
-mutable struct Node
- parent::Union{Nothing,Node} # Parent node
- branch::Branch # Current branch (fixed variables)
- branch_on::Tuple{Int,Int} # Branching index above node
- x::Vector{Float64} # Relaxation solution
- lb::Float64 # Relaxation objective value
-
- # Active set data cashed in each node
- support::Support
- GS1in_inv::Matrix{Float64}
- GSbarin_inv::Matrix{Float64}
-end
-
-mutable struct Tree
- termination_status::MOI.TerminationStatusCode # Optimization status
- start_time::DateTime # Start time
- nexpl::Int # Number of nodes explored
- queue::Vector{Node} # Nodes to explore
- active::Vector{Node} # Active nodes (with node.lb < tree.ub)
- ub::Float64 # Best upper bound
- inc::Vector{Float64} # Incumbent solution
- scr0::Int # Indices fixed to 0 with screening tests
- scr1::Int # Indices fixed to 1 with screening tests
-end
-
-
-# ======================================= #
-# Branch-and-Bound parameters and outputs #
-# ======================================= #
-
-# Available solving methods for the relaxed problem (lower bounding step)
-@enum LbMethod begin
- LBMIP # MIP solver
- ACTIVESET # Active-Set algorihtm
-end
-
-# Available heuristic methods (upper bounding step)
-@enum UbMethod begin
- UBMIP # MIP solver
- BVLS # Bounded Variable Least-Square algorithm
-end
-
-struct BnbParams
- lb_method::LbMethod # Method used to solve node relaxations
- ub_method::UbMethod # Heuristic to compute node uppper bounds
- maxtime::Int # Maximum solving time in seconds
- tol::Float64 # Integer tolerance
- maxiter::Int # Max. iter. of the relaxation solution algorithm
- screening::Bool # Toogle screening
- verbosity::Bool # Whether to enable iterations displays
- showevery::Int # Node display interval
-end
-
-struct BnbResults
- termination_status::MOI.TerminationStatusCode # Optimization status
- solve_time::Float64 # Solution time in seconds
- node_count::Int # Number of nodes explored
- objective_value::Float64 # Objective value (if optimal)
- x::Vector{Float64} # Optimizer (if optimal)
- scr0::Int # Number of screning tests passed (xi = 0)
- scr1::Int # Number of screning tests passed (xi != 0)
-end
-
-
-# ================================= #
-# MIP solver parameters and outputs #
-# ================================= #
-
-struct MipParams
- maxtime::Int # Maximum solution time in seconds
- silent::Bool # Set solver silent
-end
-
-struct MipResults
- termination_status::MOI.TerminationStatusCode # Optimization status
- solve_time::Float64 # Solution time in seconds
- node_count::Int # Number of nodes explored
- objective_value::Float64 # Objective value (if optimal)
- x::Vector{Float64} # Optimizer (if optimal)
-end
-
-
-# ============ #
-# Constructors #
-# ============ #
-
-function Support(n::Int)
- return Support(
- collect(1:n),
- Vector{Int}(),
- Vector{Int}(),
- Vector{Int}(),
- Vector{Int}()
- )
-end
-
-function ActivesetViews(node::Node, A::Matrix{Float64}, bnbparams::BnbParams)
- xSbarin = node.x[node.support.Sbarin]
- xSbarbox = node.x[node.support.Sbarbox]
- xS1in = node.x[node.support.S1in]
- xS1box = node.x[node.support.S1box]
- sSbarin = tolsign.(xSbarin, bnbparams.tol)
- sSbarbox = tolsign.(xSbarbox, bnbparams.tol)
- sS1in = tolsign.(xS1in, bnbparams.tol)
- sS1box = tolsign.(xS1box, bnbparams.tol)
- xnewSbarin = copy(xSbarin)
- xnewS1in = copy(xS1in)
- ASbarin = A[:,node.support.Sbarin]
- ASbarbox = A[:,node.support.Sbarbox]
- AS1in = A[:,node.support.S1in]
- AS1box = A[:,node.support.S1box]
- return ActivesetViews(
- xSbarin,
- xSbarbox,
- xS1in,
- xS1box,
- sSbarin,
- sSbarbox,
- sS1in,
- sS1box,
- xnewSbarin,
- xnewS1in,
- ASbarin,
- ASbarbox,
- AS1in,
- AS1box,
- )
-end
-
-function Branch(n::Int)
- return Branch(Vector{Int}(), Vector{Int}(), collect(1:n))
-end
-
-function Node(A::Matrix{Float64})
- m, n = size(A)
- return Node(
- nothing,
- Branch(n),
- (-1,-1),
- zeros(Float64,n),
- 0.,
- Support(n),
- zeros(Float64,0,0),
- zeros(Float64,0,0),
- )
-end
-
-function Node(parent::Node, j::Int, jval::Int)
- node = Node(
- parent,
- copy(parent.branch),
- (j,jval),
- copy(parent.x),
- parent.lb,
- copy(parent.support),
- copy(parent.GS1in_inv),
- copy(parent.GSbarin_inv),
- )
- fixto!(node,j,jval)
- return node
-end
-
-function Tree(A::Array{Float64,2}, y::Vector{Float64})
- return Tree(
- MOI.OPTIMIZE_NOT_CALLED,
- Dates.now(),
- 0,
- Vector{Node}([Node(A)]),
- Vector{Node}(),
- 0.5 * norm(y, 2)^2,
- zeros(Float64,size(A)[2]),
- 0,
- 0,
- )
-end
-
-function BnbParams(;
- lb_method=ACTIVESET,
- ub_method=BVLS,
- maxtime=1000,
- tol=1e-8,
- maxiter=1000,
- screening=false,
- verbosity=true,
- showevery=1,
- )
- maxtime > 0 || error("Parameter `maxtime` must be positive")
- tol > 0 || error("Parameter `tol` must be positive")
- maxiter > 0 || error("Parameter `maxiter` must be positive")
- showevery > 0 || error("Parameter `showevery` must be positive")
- return BnbParams(
- lb_method,
- ub_method,
- maxtime,
- tol,
- maxiter,
- screening,
- verbosity,
- showevery
- )
-end
-
-function BnbResults(tree::Tree)
- return BnbResults(
- tree.termination_status,
- elapsed_time(tree),
- tree.nexpl,
- tree.ub,
- tree.inc,
- tree.scr0,
- tree.scr1,
- )
-end
-
-function MipParams(;maxtime=1000,silent=true)
- maxtime > 0 || error("Parameter `maxtime` must be positive")
- return MipParams(maxtime, silent)
-end
-
-function MipResults(problem::JuMP.Model)
- return MipResults(
- termination_status(problem),
- solve_time(problem),
- convert(Int,node_count(problem)),
- objective_value(problem),
- value.(problem[:x])
- )
-end
-
-
-# ====================== #
-# Base module extensions #
-# ====================== #
-
-
-function Base.copy(support::Support)
- return Support(
- copy(support.Sbar0),
- copy(support.Sbarin),
- copy(support.Sbarbox),
- copy(support.S1in),
- copy(support.S1box)
- )
-end
-
-function Base.copy(branch::Branch)
- return Branch(copy(branch.S0), copy(branch.S1), copy(branch.Sbar))
-end
diff --git a/tests/test.jl b/tests/test.jl
deleted file mode 100644
index b204ae1817c4350e9d82e6f74e2db72f90f32257..0000000000000000000000000000000000000000
--- a/tests/test.jl
+++ /dev/null
@@ -1,31 +0,0 @@
-include("../src/BnbScreening.jl")
-
-k, m, n, snr = 5, 500, 300, 10*log10(10)
-A, y, λ, M = data_toeplitz(k,m,n,snr)
-#k, m, n, snr = 5, 500, 1000, 10*log10(10)
-#A, y, λ, M = data_gaussian(k,m,n,snr)
-
-result_mip = solve_mip(A,y,λ,M)
-println("MIP")
-println(" Status : $(result_mip.termination_status)")
-println(" Objective : $(result_mip.objective_value)")
-println(" Nodes : $(result_mip.node_count)")
-println(" Timer : $(result_mip.solve_time)")
-
-result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(verbosity=false))
-println("BnB")
-println(" Status : $(result_bnb.termination_status)")
-println(" Objective : $(result_bnb.objective_value)")
-println(" Nodes : $(result_bnb.node_count)")
-println(" Timer : $(result_bnb.solve_time)")
-println("Error : $(abs(result_mip.objective_value - result_bnb.objective_value))")
-
-result_bnb_scr = solve_bnb(A,y,λ,M,bnbparams=BnbParams(verbosity=false,screening=true))
-println("BnB + scr")
-println(" Status : $(result_bnb_scr.termination_status)")
-println(" Objective : $(result_bnb_scr.objective_value)")
-println(" Nodes : $(result_bnb_scr.node_count)")
-println(" Timer : $(result_bnb_scr.solve_time)")
-println(" Scr-0 : $(result_bnb_scr.scr0)")
-println(" Scr-1 : $(result_bnb_scr.scr1)")
-println("Error : $(abs(result_mip.objective_value - result_bnb_scr.objective_value))")
\ No newline at end of file