M .gitignore => .gitignore +2 -0
@@ 4,3 4,5 @@
!simulationplot.png
plotsim.gp
cmocean_balance.pal
+debug_*
+*.swp
A LICENSE.txt => LICENSE.txt +336 -0
@@ 0,0 1,336 @@
+GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 2022, Chris Geoga
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Lesser General Public License instead.) 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
+this service 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 make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. 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.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute 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 and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+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
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the 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 a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, 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.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE 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.
+
+ 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
+convey 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 2 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, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision 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, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This 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.
M Manifest.toml => Manifest.toml +16 -96
@@ 1,106 1,26 @@
# This file is machine-generated - editing it directly is not advised
-[[AbstractFFTs]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
-uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "0.5.0"
+julia_version = "1.7.1"
+manifest_format = "2.0"
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-[[FFTW]]
-deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
-git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf"
-uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.2.2"
-
-[[FFTW_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b"
-uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
-version = "3.3.9+5"
-
-[[IntelOpenMP_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6"
-uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2018.0.3+0"
-
-[[InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[LibGit2]]
-deps = ["Printf"]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[Libdl]]
+[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-[[LinearAlgebra]]
-deps = ["Libdl"]
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[MKL_jll]]
-deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08"
-uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2020.1.216+0"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[Reexport]]
-deps = ["Pkg"]
-git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
-uuid = "189a3867-3050-52da-a836-e630ba90ab69"
-version = "0.2.0"
-
-[[SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
M Project.toml => Project.toml +1 -3
@@ 1,9 1,7 @@
name = "HalfSpectral"
uuid = "28fda99a-6fd1-4789-996a-91dd14997ef0"
authors = ["Chris Geoga <cgeoga@protonmail.com>"]
-version = "0.1.0"
+version = "0.2.0"
[deps]
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
M README.md => README.md +17 -17
@@ 4,31 4,31 @@
A convenient interface for working with half-spectral covariance functions in
the space-time domain. This is part of the software companion to [Flexible
nonstationary spatio-temporal modeling of high-frequency monitoring
-data](http://arxiv.org/abs/2007.11418).
+data](https://onlinelibrary.wiley.com/doi/10.1002/env.2670).
See the paper for a discussion of the covariance function itself, and see the
-example file `./example/basicmodel.jl` for an example specification of the
-model, `./example/basicexample.jl` for an example of creating a kernel function
-and then simulating a process with it, and finally `./example/fitsimulation.jl`
-for an example of computing an MLE for the simulated data. With very few lines
-of boilerplate, you can specify a complicated covariance function corresponding
-to fields that look like this:
+example file `./example/basicexample.jl` for an example of creating a kernel
+that you can then call as a regular function. The example file
+`./example/withvecchia.jl` gives an example of the wrapper type that you can use
+to conveniently use `ForwardDiff.jl` to get derivatives of the kernel function
+as well as an example of how to efficiently compose this package with
+`Vecchia.jl` for scalable parameter estimation.
+
+With very few lines of boilerplate, you can specify a complicated covariance
+function corresponding to fields that look like this:
<p align="center">
<img src="https://git.sr.ht/~cgeoga/HalfSpectral.jl/blob/master/simulationplot.png">
</p>
-For the code used to fit models in the paper, see `./examples/paperscripts/`. I
-do not plan to modify these files if this code evolves, so if you want to run
-those scripts, please use the appropriate tagged release.
-
-For more information about estimation, please see another of my software
-packages, [GPMaxlik.jl](https://git.sr.ht/~cgeoga/GPMaxlik.jl). All the
-necessary code to try simulating and estimating is here, but several more
-complete examples are provided in that repository.
+**If you use this software (HalfSpectral.jl) in a manuscript or published work,
+please cite the paper linked above and not the software package itself.**
-**If you use this software (HalfSpectral.jl), please cite the paper and not
-the software package itself.**
+This release (0.2+) is a complete re-write from the version that was used in the
+paper (0.1), so please go to an earlier commit if you're interested in
+re-creating paper results. But unless you're doing that, I'd really strongly
+suggest using the new version of the code, because it is much cleaner and has
+complete AD compatibility.
# Installation
A example/Manifest.toml => example/Manifest.toml +583 -0
@@ 0,0 1,583 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.1"
+manifest_format = "2.0"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.3.3"
+
+[[deps.ArgCheck]]
+git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
+uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
+version = "2.3.0"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.ArrayInterface]]
+deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "Static"]
+git-tree-sha1 = "ec8a5e8528995f2cec48c53eb834ab0d58f8bd99"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "6.0.14"
+
+[[deps.ArrayInterfaceCore]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "d0f59ebfe8d3ea2799fb3fb88742d69978e5843e"
+uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
+version = "0.1.10"
+
+[[deps.ArrayInterfaceOffsetArrays]]
+deps = ["ArrayInterface", "OffsetArrays", "Static"]
+git-tree-sha1 = "7dce0e2846e7496622f5d2742502d7e029693458"
+uuid = "015c0d05-e682-4f19-8f0a-679ce4c54826"
+version = "0.1.5"
+
+[[deps.ArrayInterfaceStaticArrays]]
+deps = ["Adapt", "ArrayInterface", "LinearAlgebra", "Static", "StaticArrays"]
+git-tree-sha1 = "d7dc30474e73173a990eca86af76cae8790fa9f2"
+uuid = "b0d46f97-bff5-4637-a19a-dd75974142cd"
+version = "0.1.2"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.BangBang]]
+deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
+git-tree-sha1 = "b15a6bc52594f5e4a3b825858d1089618871bf9d"
+uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
+version = "0.3.36"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.Baselet]]
+git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
+uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
+version = "0.1.1"
+
+[[deps.BitTwiddlingConvenienceFunctions]]
+deps = ["Static"]
+git-tree-sha1 = "28bbdbf0354959db89358d1d79d421ff31ef0b5e"
+uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
+version = "0.1.3"
+
+[[deps.CPUSummary]]
+deps = ["CpuId", "IfElse", "Static"]
+git-tree-sha1 = "0eaf4aedad5ccc3e39481db55d72973f856dc564"
+uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
+version = "0.1.22"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "9489214b993cd42d17f44c36e359bf6a7c919abf"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.15.0"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.3"
+
+[[deps.CloseOpenIntervals]]
+deps = ["ArrayInterface", "Static"]
+git-tree-sha1 = "eb61d6b97041496058245821e3bb7eba2b2cf4db"
+uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
+version = "0.1.8"
+
+[[deps.CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
+[[deps.Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "9be8be1d8a6f44b96482c8af52238ea7987da3e3"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.45.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.CompositionsBase]]
+git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
+uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
+version = "0.1.1"
+
+[[deps.ConstructionBase]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.3.0"
+
+[[deps.ContextVariablesX]]
+deps = ["Compat", "Logging", "UUIDs"]
+git-tree-sha1 = "8ccaa8c655bc1b83d2da4d569c9b28254ababd6e"
+uuid = "6add18c4-b38d-439d-96f6-d6bc489c04c5"
+version = "0.1.2"
+
+[[deps.CpuId]]
+deps = ["Markdown"]
+git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406"
+uuid = "adafc99b-e345-5852-983c-f28acb93d879"
+version = "0.3.1"
+
+[[deps.DataAPI]]
+git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.10.0"
+
+[[deps.DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DefineSingletons]]
+git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
+uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
+version = "0.1.2"
+
+[[deps.DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[deps.DiffResults]]
+deps = ["StaticArrays"]
+git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.0.3"
+
+[[deps.DiffRules]]
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.11.0"
+
+[[deps.Distances]]
+deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
+uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+version = "0.10.7"
+
+[[deps.Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.FLoops]]
+deps = ["BangBang", "Compat", "FLoopsBase", "InitialValues", "JuliaVariables", "MLStyle", "Serialization", "Setfield", "Transducers"]
+git-tree-sha1 = "4391d3ed58db9dc5a9883b23a0578316b4798b1f"
+uuid = "cc61a311-1640-44b5-9fba-1b764f453329"
+version = "0.2.0"
+
+[[deps.FLoopsBase]]
+deps = ["ContextVariablesX"]
+git-tree-sha1 = "656f7a6859be8673bf1f35da5670246b923964f7"
+uuid = "b9860ae5-e623-471e-878b-f6a53c775ea6"
+version = "0.1.1"
+
+[[deps.ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "2f18915445b248731ec5db4e4a17e451020bf21e"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.30"
+
+[[deps.Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[deps.HostCPUFeatures]]
+deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
+git-tree-sha1 = "18be5268cf415b5e27f34980ed25a7d34261aa83"
+uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
+version = "0.1.7"
+
+[[deps.IfElse]]
+git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
+uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+version = "0.1.1"
+
+[[deps.InitialValues]]
+git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3"
+uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
+version = "0.3.1"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "b3364212fb5d870f724876ffcd34dd8ec6d98918"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.7"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.JuliaVariables]]
+deps = ["MLStyle", "NameResolution"]
+git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
+uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec"
+version = "0.2.4"
+
+[[deps.LayoutPointers]]
+deps = ["ArrayInterface", "ArrayInterfaceOffsetArrays", "ArrayInterfaceStaticArrays", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"]
+git-tree-sha1 = "a575de5a424a395217930fea6d0934ea853d0158"
+uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
+version = "0.1.9"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.15"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.LoopVectorization]]
+deps = ["ArrayInterface", "ArrayInterfaceCore", "ArrayInterfaceOffsetArrays", "ArrayInterfaceStaticArrays", "CPUSummary", "ChainRulesCore", "CloseOpenIntervals", "DocStringExtensions", "ForwardDiff", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "SIMDDualNumbers", "SIMDTypes", "SLEEFPirates", "SpecialFunctions", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"]
+git-tree-sha1 = "b11116837200f222c00d2b14a594a98ab3c97a0d"
+uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
+version = "0.12.116"
+
+[[deps.MLStyle]]
+git-tree-sha1 = "2041c1fd6833b3720d363c3ea8140bffaf86d9c4"
+uuid = "d8e11817-5142-5d16-987a-aa16d5891078"
+version = "0.4.12"
+
+[[deps.MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.9"
+
+[[deps.ManualMemory]]
+git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
+uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
+version = "0.1.8"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.MicroCollections]]
+deps = ["BangBang", "InitialValues", "Setfield"]
+git-tree-sha1 = "6bb7786e4f24d44b4e29df03c69add1b63d88f01"
+uuid = "128add7d-3638-4c79-886c-908ea0c25c34"
+version = "0.1.2"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NaNMath]]
+git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.0"
+
+[[deps.NameResolution]]
+deps = ["PrettyPrint"]
+git-tree-sha1 = "1a0fa0e9613f46c9b8c11eee38ebb4f590013c5e"
+uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391"
+version = "0.1.5"
+
+[[deps.HalfSpectral]]
+deps = ["LinearAlgebra"]
+path = "/home/cg/Personal/Repos/HalfSpectral.jl/"
+uuid = "a4cc074a-5bcd-4a73-94df-f94ee43b423f"
+version = "0.1.0"
+
+[[deps.NearestNeighbors]]
+deps = ["Distances", "StaticArrays"]
+git-tree-sha1 = "ded92de95031d4a8c61dfb6ba9adb6f1d8016ddd"
+uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+version = "0.4.10"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "b4975062de00106132d0b01b5962c09f7db7d880"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.12.5"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[deps.PolyesterWeave]]
+deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
+git-tree-sha1 = "7e597df97e46ffb1c8adbaddfa56908a7a20194b"
+uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad"
+version = "0.1.5"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.PrettyPrint]]
+git-tree-sha1 = "632eb4abab3449ab30c5e1afaa874f0b98b586e4"
+uuid = "8162dcfd-2161-5ef2-ae6c-7681170c5f98"
+version = "0.2.0"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.SIMDDualNumbers]]
+deps = ["ForwardDiff", "IfElse", "SLEEFPirates", "VectorizationBase"]
+git-tree-sha1 = "dd4195d308df24f33fb10dde7c22103ba88887fa"
+uuid = "3cdde19b-5bb0-4aaf-8931-af3e248e098b"
+version = "0.1.1"
+
+[[deps.SIMDTypes]]
+git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
+uuid = "94e857df-77ce-4151-89e5-788b33177be4"
+version = "0.1.0"
+
+[[deps.SLEEFPirates]]
+deps = ["IfElse", "Static", "VectorizationBase"]
+git-tree-sha1 = "ac399b5b163b9140f9c310dfe9e9aaa225617ff6"
+uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
+version = "0.6.32"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.Setfield]]
+deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
+git-tree-sha1 = "38d88503f695eb0301479bc9b0d4320b378bafe5"
+uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+version = "0.8.2"
+
+[[deps.SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "a9e798cae4867e3a41cae2dd9eb60c047f1212db"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.1.6"
+
+[[deps.SplittablesBase]]
+deps = ["Setfield", "Test"]
+git-tree-sha1 = "39c9f91521de844bad65049efd4f9223e7ed43f9"
+uuid = "171d559e-b47b-412a-8079-5efa626c420e"
+version = "0.1.14"
+
+[[deps.Static]]
+deps = ["IfElse"]
+git-tree-sha1 = "5d2c08cef80c7a3a8ba9ca023031a85c263012c5"
+uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+version = "0.6.6"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "383a578bdf6e6721f480e749d503ebc8405a0b22"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.4.6"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.StatsAPI]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "2c11d7290036fe7aac9038ff312d3b3a2a5bf89e"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.4.0"
+
+[[deps.SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[deps.Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
+git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.7.0"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.ThreadingUtilities]]
+deps = ["ManualMemory"]
+git-tree-sha1 = "f8629df51cab659d70d2e5618a430b4d3f37f2c3"
+uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
+version = "0.5.0"
+
+[[deps.Transducers]]
+deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
+git-tree-sha1 = "c76399a3bbe6f5a88faa33c8f8a65aa631d95013"
+uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
+version = "0.4.73"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.UnPack]]
+git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
+uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+version = "1.0.2"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Vecchia]]
+deps = ["BangBang", "FLoops", "LinearAlgebra", "LoopVectorization", "MicroCollections", "NearestNeighbors", "SparseArrays", "StaticArrays"]
+path = "/home/cg/Personal/Repos/Vecchia.jl/"
+uuid = "8d73829f-f4b0-474a-9580-cecc8e084068"
+version = "0.6.1"
+
+[[deps.VectorizationBase]]
+deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"]
+git-tree-sha1 = "7d3de169cd221392082a5abc7f363726e1a30628"
+uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
+version = "0.21.36"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.ZygoteRules]]
+deps = ["MacroTools"]
+git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
+uuid = "700de1a5-db45-46bc-99cf-38207098b444"
+version = "0.2.2"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
A example/Project.toml => example/Project.toml +5 -0
@@ 0,0 1,5 @@
+[deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HalfSpectral = "a4cc074a-5bcd-4a73-94df-f94ee43b423f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Vecchia = "8d73829f-f4b0-474a-9580-cecc8e084068"
A example/basic_example.jl => example/basic_example.jl +28 -0
@@ 0,0 1,28 @@
+
+using HalfSpectral, StaticArrays
+
+# Simple SDF and a simple nontrivial coherence, which make an integrand:
+sdf(f,x,p) = p[1]*(p[2]^2 + f^2)^(-p[3] - 1/2)
+coh(f,x,y,p) = exp(-p[4]*(1+(p[5]*f)^2)*abs(x-y))
+integrand(f,x,y,p) = coh(f,x,y,p)*sqrt(sdf(f,x,p)*sdf(f,y,p))
+
+# Let's make some basic points. Say, 5 measurements at x=1:10 at time points
+# 1:100 each. Importantly, the time has to come first.
+if !(@isdefined pts)
+ const xpts = 1:30
+ const timelen = 200
+ const pts = reduce(vcat, map(x->[@SVector [t, x] for t in 1:timelen], xpts))
+end
+
+# Now let's make the kernel:
+if !(@isdefined kernel)
+ const params = ones(5)
+ const kernel = HalfSpectral.HSKernel(params, xpts, integrand, timelen)
+end
+
+# Now just call it like a function! Note that you can't call it for locations
+# that aren't in pts, though.
+sample_call = kernel(pts[1], pts[3])
+
+
+
D example/basicexample.jl => example/basicexample.jl +0 -52
@@ 1,52 0,0 @@
-
-using HalfSpectral, LinearAlgebra
-
-# Load a basic specification of a nonstationary half-spectral model.
-include("basicmodel.jl")
-
-# Create the kernel for spatial locations 1:30 and 100 unit time measurements.
-const tg = timegrid(1:100)
-const xpts = 1:30
-const params = (log(50.0),log(100.0),log(0.5),log(5.0),0.5) # sample parameter values.
-kernel = tftkernel(integrand, tg, xpts, params)
-
-# Build a second time to get a timing example.
-println("Building the kernel takes this long")
-@time kernel = tftkernel(integrand, tg, xpts, params)
-
-# Simulate a realization of the field with this function. The points here should
-# be ordered as a tuple (time, spatial coordinate), and the spatial coordinate
-# does not need to be a real number. You can treat the kernel like any normal
-# function, calling it with either two points x and y like kernel(x,y) or
-# additionally with a parameter collection, and if the parameter collection
-# disagrees with the internal one it will rebuild the dictionary of function
-# values. In general, I suggest including the parameter vector for code safety.
-# But to avoid runtime checks, you can use the unchecked kernel(x,y).
-pts = vec(reverse.(collect(Iterators.product(xpts, tg.Xt))))
-K = Symmetric([kernel(x,y) for x in pts, y in pts])
-simv = cholesky(K).L*randn(length(pts))
-simm = reshape(simv, length(xpts), length(tg.Xt))
-
-# Optionally, save the simulation file.
-
-#= If you just want a csv file to plot.
-using DelimitedFiles
-writedlm("sim.csv", simm, ',')
-=#
-
-# If you want to save enough output to also fit your simulated data:
-using JLD
-JLD.save("simulation_output.jld",
- "true_params", params,
- "pts", pts,
- "tgrid", tg,
- "xpts", xpts,
- "data", simv)
-# =#
-# =#
-
-#= Optionally, visualize it.
-using Plots
-heatmap(simm, xlabel="time", ylabel="altitude")
-=#
-
D example/basicmodel.jl => example/basicmodel.jl +0 -31
@@ 1,31 0,0 @@
-
-# Marginal spectral density, with the twist of making a few parameters depend on
-# the spatial index.
-function Sf(w,x,p)
- (scale, rate, smoothness) = (exp(p[1]), exp(p[2]), exp(p[3]))
- rate += log(x)
- smoothness *= (1+x/10)
- scale*(rate*sinpi(w)^2 + 1)^(-smoothness-1/2)
-end
-
-# A coherence function that does happen to be stationary. See the paperscripts
-# model for code implementing the Paciorek-Schervish nonstationary coherence.
-function Cw(w,x,y,p)
- gamxy = exp(p[4])/(1 + (sinpi(w)/0.3)^2)^4
- exp(-abs(x-y)/gamxy)
-end
-
-# A simple phase function and it's derivative with respect to kernel parameters.
-phase(w,x,y,p) = exp(2.0*pi*im*p[5]*sin(sinpi(w))*(x-y))
-dphase(w,x,y,p)= 2.0*pi*im*sin(sinpi(w))*(x-y)*phase(w,x,y,p)
-
-# The integrand in the form that HalfSpectral requests, splitting up the phase
-# part so we can use AD on the real-valued parts, which seems to avoid a lot of
-# finnicky issues with ForwardDiff.jl. This is not strictly necessary, but as
-# you'll see in fitsimulation.jl, writing it in this slightly tricky way means
-# that with just a few more lines of code, you can use very powerful
-# quasi-Newton optimization methods.
-integrand_nophase(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)
-integrand(w,x,y,p) = integrand_nophase(w,x,y,p)*phase(w,x,y,p)
-integrand_dphase(w,x,y,p) = integrand_nophase(w,x,y,p)*dphase(w,x,y,p)
-
D example/fitsimulation.jl => example/fitsimulation.jl +0 -58
@@ 1,58 0,0 @@
-
-using HalfSpectral, GPMaxlik, JLD, ForwardDiff
-
-# Load in the model used to simulate the data.
-include("basicmodel.jl")
-
-# A lazy version of kernel derivatives: use AD on the integrand. This would
-# obviously be much faster if you hand-wrote your derivative functions for the
-# integrand, but as you'll see here this is good enough for prototyping.
-function coord_deriv(f,p,j)
- ForwardDiff.derivative(z->f(vcat(p[1:(j-1)],z,p[(j+1):end])), p[j])
-end
-dfuns = convert(Vector{Function},
- [(w,x,y,p)->coord_deriv(z->integrand_nophase(w,x,y,z),p,j)*phase(w,y,x,p)
- for j in 1:4])
-push!(dfuns, integrand_dphase)
-
-# Load the saved data in:
-const pts = JLD.load("simulation_output.jld", "pts")
-const tgrid = JLD.load("simulation_output.jld", "tgrid")
-const xpts = JLD.load("simulation_output.jld", "xpts")
-const dat = JLD.load("simulation_output.jld", "data")
-const tru = JLD.load("simulation_output.jld", "true_params")
-
-# Choose some initial conditions, then create the kernel function and the
-# derivative kernel functions.
-ini = vcat(log.([75.0, 75.0, 1.0, 4.0]), 0.75)
-kfun = tftkernel(integrand, tgrid, xpts, ini)
-dkfuns = [tftkernel(integrand_dj, tgrid, xpts, ini) for integrand_dj in dfuns]
-
-# Create some SAA vectors for faster derivatives.
-const saa = rand([-1.0, 1.0], length(pts), 100)
-
-# Create objective only and objective+grad+fisher functions for GPMaxlik, using
-# the optimized nll function it provides.
-obj(p) = nll(pts, dat, kfun, dkfuns, p, saa=saa, nll=true, grad=false, fish=false).nll
-function objgradfish(p)
- nll_call = nll(pts, dat, kfun, dkfuns, p, saa=saa, nll=true, grad=true, fish=true)
- (nll_call.nll, nll_call.grad, nll_call.fish)
-end
-
-# Perform the optimization using the trust region method (see GPMaxlik.jl for
-# details) using the stochastic expected Fisher matrix as a Hessian
-# approximation. If you want to try BFGS to compare, just change "fish=true" to
-# "fish=false" in line 38 above.
-maxlik_result = trustregion(obj, objgradfish, ini, vrb=true, maxit=100, rtol=1.0e-5)
-
-# Consider the resulting MLE, std approximations from the stochastic expected
-# Fisher matrix, and compare with the true parameters for the simulation.
-mle = maxlik_result.p
-stds = map(x->"("*string(x)*")",
- round.(sqrt.(diag(inv(maxlik_result.h))).*1.96, digits=2))
-table = hcat(pushfirst!(string.(round.(collect(tru), digits=4)), "true"),
- pushfirst!(string.(round.(collect(mle), digits=4)), "mle"),
- pushfirst!(stds, "±95%"))
-display(table)
-
-
D example/paperscripts/README.txt => example/paperscripts/README.txt +0 -12
@@ 1,12 0,0 @@
-
-These files are most of the code used to fit the model described in the paper.
-They will not "just run" on your computer as I have removed several hard paths
-for storing output and reading data files. Nonetheless, I think they are helpful
-in illustrating a more advanced workflow for several data sources and a
-reasonably complex model.
-
-I will NOT be updating these files to reflect potential breaking changes in
-HalfSpectral.jl, so please refer to the commit tagged "paper" to see a version
-of HalfSpectral.jl that was the source code that these scripts actually ran when
-generating results for the paper.
-
D => +0 -24
@@ 1,24 0,0 @@
using HalfSpectral
# UNITS:
# time (ix)
# altitudes (ix)
# data (m/s)
function extract_block(trange, times, gates, dopp, ints, fn)
qj = findall(x->trange[1]<=x<=trange[2], times)
minr = minimum.(eachrow(ints[:,qj]))
snrj = filter(x->in(x, gates), findall(x->x>1.008, minr))
tgrd = timegrid(times[qj])
tgrd isa TimeFFTKernels.TimeGridFailure && return (tgrd, nothing)
pts = reverse.(vec(collect(Iterators.product(snrj.*30.0, tgrd.Xt))))
dat = vec(dopp[snrj,qj])
if fn isa Nothing
mldat = (pts=pts, dat=dat)
else
println("Transforming data...")
mldat = (pts=pts, dat=fn(dat))
end
(tgrd, snrj.*30.0, mldat) # (timegrid object, altitude (m), named tuple)
end
D example/paperscripts/fit_trustregion.jl => example/paperscripts/fit_trustregion.jl +0 -107
@@ 1,107 0,0 @@
-
-using HDF5, HalfSpectral, LinearAlgebra, DelimitedFiles, GPMaxlik, JLD, Random
-
-# Set the seed:
-Random.seed!(123456)
-
-# At this point, do five iterations, then start using an exact gradient:
-GPMaxlik.SETTINGS.GRAD_INFNORM_TOL = 0.5
-
-include("model.jl")
-include("extractor.jl")
-include("writer.jl")
-include("nsscale.jl")
-
-# no build indices: re-building the integrand kernel object isn't necessary for
-# these parameter indices.
-const nob = (1,2,3,4,24,25)
-
-# Create a NonstationaryScale object:
-NS = NonstationaryScale{Float64,false,0,4}((156, 311, 466, 621), (1,2,3,4))
-dNS = [NonstationaryScale{Float64,true,j,4}((156, 311, 466, 621), (1,2,3,4)) for j in 1:4]
-
-# Some simple wrapper functions:
-function nll_grad_fish(pts, data, kfn, dfns, p, saa)
- nll_call = nll(pts, data, kfn, dfns, p, saa=saa, nll=true, grad=true, fish=true)
- (nll_call.nll, nll_call.grad, nll_call.fish)
-end
-
-# Convert decimal hours to minutes:
-decimalhours_to_minutes(dhours) = (dhours .- 14.0).*60.0
-
-# An extra function to print the Newton step every five iterations:
-function print_ns(c, fx, xv, gx, hx, dl; kwargs...)
- if iszero(rem(c, 5))
- println("NS: $(round.(hx\gx, digits=4))")
- end
- return (false, nothing)
-end
-
-function fit_gates_trustregion(day, timewin, gates, its, ini, rt, gt,
- trfun=nothing, sim=false, dc=1.0e-4)
- # Load in the data:
- println("Working with day $day")
- rdata = read(h5open(__YOUR_PATH_TO_RAW_DATA__)) # Put your path here
- (times, dopp, ints) = (rdata[x] for x in ("Time", "Doppler", "Intensity"))
- (tgrid, alts, data) = extract_block(timewin, times, gates, dopp, ints, trfun)
- if unique(diff(alts)) != [30.0]
- @warn "You are missing gates in the range you've requested."
- end
- # Create the necessary objects:
- saav = rand([-1.0, 1.0], length(data.pts), 250) # was 125 for first stage of fitting
- kfun = tftkernel(integrand, tgrid, alts, ini, addfn=nugfn, mulfn=NS, nob=nob)
- dfns = vcat(
- [tftkernel(integrand, tgrid, alts, ini, mulfn=dNSj, nob=nob) for dNSj in dNS],
- [tftkernel(dk, tgrid, alts, ini, mulfn=NS, nob=nob) for dk in dfuns],
- dnugfn1, dnugfn2
- )
- # Do the optimization:
- obj = p-> nll(data.pts, data.dat, kfun, dfns, p, saa=saav,
- nll=true, grad=false, fish=false).nll
- objgradfish = p -> nll_grad_fish(data.pts, data.dat, kfun, dfns, p, saav)
- maxlik_result = trustregion(obj, objgradfish, ini, vrb=true, maxit=its, gtol=gt,
- rtol=rt, dcut=dc, iter_funs=(GPMaxlik.convg_info, print_ns))
- println("Estimation result: $(maxlik_result.status)")
- println()
- if sim
- L = cholesky!(Symmetric([kfun(x, y) for x in data.pts, y in data.pts])).L
- d = reshape(data.dat, length(alts), length(tgrid.Tv))
- s = reshape(L*randn(length(data.pts)), length(alts), length(tgrid.Tv))
- gnuplot_save_matrix!(__YOUR_PATH__*"sim_"*day*"_14.csv", s, alts, # Put your path here
- decimalhours_to_minutes(tgrid.Tv))
- gnuplot_save_matrix!(__YOUR_PATH__*"tru_"*day*"_14.csv", d, alts, # Put your path here
- decimalhours_to_minutes(tgrid.Tv))
- end
- return maxlik_result
-end
-
-
-
-
-
-
-const FITTED_DAYS = tuple() # Satisfactorily fitted days
-const tw = (14.02, 14.24)
-
-for (day, gates, mit) in zip(("02", "03", "06", "20", "24", "28"),
- (7:23, 7:28, 7:23, 7:28, 7:28, 7:28),
- (30, 30, 30, 30, 30, 30))
- if in(day, FITTED_DAYS)
- println("Estimate for day $day was marked acceptable, moving on.")
- println()
- continue
- end
- init = load(__YOUR_PATH__*"fit_"*day*"_results.jld")["result"].p
- try
- result = fit_gates_trustregion(day, tw, gates, mit, init,
- 1.0e-6, 1.0e-3, nothing, true, 1.0e-7)
- # Deal with results:
- GPMaxlik.reset_settings!()
- JLD.save(__YOUR_PATH__*"fit_"*day*"_results.jld", "result", result)
- catch
- println()
- println("Day $day broke somehow. Check out the logs.")
- println()
- end
-end
-
D example/paperscripts/fitted_estimates_untransformed.jl => example/paperscripts/fitted_estimates_untransformed.jl +0 -65
@@ 1,65 0,0 @@
-
-using HDF5, TimeFFTKernels, LinearAlgebra, DelimitedFiles, GPMaxlik, JLD, Random
-
-# Set the seed:
-Random.seed!(123456)
-
-# At this point, do five iterations, then start using an exact gradient:
-GPMaxlik.SETTINGS.GRAD_INFNORM_TOL = 0.5
-
-include("model_untransformed.jl")
-include("extractor.jl")
-include("writer.jl")
-include("nsscale.jl")
-include("untransformer.jl")
-
-# no build indices: re-building the integrand kernel object isn't necessary for
-# these parameter indices.
-const nob = (1,2,3,4,23,24,25)
-
-# Create a NonstationaryScale object:
-NS = NonstationaryScale{Float64,false,0,4}((156, 311, 466, 621), (1,2,3,4))
-dNS = [NonstationaryScale{Float64,true,j,4}((156, 311, 466, 621), (1,2,3,4)) for j in 1:4]
-
-# Some simple wrapper functions:
-function nll_grad_fish(pts, data, kfn, dfns, p, saa)
- nll_call = nll(pts, data, kfn, dfns, p, saa=saa, nll=true, grad=true, fish=true)
- (nll_call.nll, nll_call.grad, nll_call.fish)
-end
-
-function untransformed_uncertainty(day, timewin, gates, trfun=nothing)
- # Load in the data:
- println("Working with day $day")
- rdata = read(h5open(__YOUR_PATH_TO_RAW_DATA)) # your path to the raw files
- (times, dopp, ints) = (rdata[x] for x in ("Time", "Doppler", "Intensity"))
- (tgrid, alts, data) = extract_block(timewin, times, gates, dopp, ints, trfun)
- # Load in the MLE and un-transform it:
- mle_result = load(__YOUR_PATH__*"fit_"*day*"_results.jld")["result"]
- mle = transf(mle_result.p)
- # Create the necessary objects:
- saav = rand([-1.0, 1.0], length(data.pts), 125)
- kfun = tftkernel(integrand, tgrid, alts, mle, addfn=nugfn, mulfn=NS, nob=nob)
- dfns = vcat(
- [tftkernel(integrand, tgrid, alts, mle, mulfn=dNSj, nob=nob) for dNSj in dNS],
- [tftkernel(dk, tgrid, alts, mle, mulfn=NS, nob=nob) for dk in dfuns],
- dnugfn1, dnugfn2
- )
- # Obtain the fisher information matrix for the untransformed model at the MLE:
- out = nll_grad_fish(data.pts, data.dat, kfun, dfns, mle, saav)
- println("Likelihood at MLE for un-transformed data: $(round(out[1], digits=5))")
- (mle, out)
-end
-
-const tw = (14.02, 14.24)
-const gater = (7:23, 7:28, 7:23, 7:28, 7:28, 7:28)
-for (day, gates) in zip(("02", "03", "06", "20", "24", "28"), gater)
- # Obtain the un-transformed MLE and uncertainty:
- untransf_result = untransformed_uncertainty(day, tw, gates, nothing)
- # Deal with results:
- GPMaxlik.reset_settings!()
- JLD.save(__YOUR_PATH__*"fit_"*day*"_results_untransformed.jld",
- "mle", untransf_result[1],
- "nllresult", untransf_result[2])
-end
-
-
D example/paperscripts/generate_nls_init.jl => example/paperscripts/generate_nls_init.jl +0 -169
@@ 1,169 0,0 @@
-
-using DelimitedFiles, Statistics, HDF5, SMultitaper, Optim, ForwardDiff, HalfSpectral
-
-# Load in the model and the extractor function.
-include("model.jl")
-include("extractor.jl")
-
-# A function that reads in the data and parses it into the format I want. In
-# this case, a Dict, indexed by gate _number_, not altitude, since indexing a
-# hashmap by a floating point number is bad style (to m
-function extract_day_slice(day, timewin, gates)
- # Load in the data:
- rdata = read(h5open("../../data/raw/Vert_2015_06_"*day*"_14.h5"))
- (times, dopp, ints) = (rdata[x] for x in ("Time", "Doppler", "Intensity"))
- (tgrid, alts, data) = extract_block(timewin, times, gates, dopp, ints,
- nothing, asmatrix=true)
- if unique(diff(alts)) != [30.0]
- @warn "You are missing gates in the range you've requested."
- end
- # Reshape to a dict for easier reading in the fitting:
- datadict = Dict{Int64, Vector{Float64}}()
- for (j, gatej) in enumerate(gates)
- datadict[gatej] = data[2][j,:]
- end
- datadict
-end
-
-function cross_specs(datadict; coherence=true)
- out_type = coherence ? Complex{Float64} : Float64
- crosses = Dict{Tuple{Int64, Int64}, Vector{out_type}}()
- gates = sort(collect(keys(datadict)))
- for gj in gates
- for gk in (gj+1):gates[end]
- if coherence
- cross = real(multitaper(datadict[gj], datadict[gk], nw=15.0, k=25,
- coherence=true, shift=true, center=true))
- else
- cross = multitaper(datadict[gj], datadict[gk], nw=5.0, k=8,
- aweight=true, shift=true, center=true)
- end
- crosses[(gj, gk)] = cross
- end
- end
- crosses
-end
-
-function marginal_specs(datadict)
- specs = Dict{Int64, Vector{Float64}}()
- gates = sort(collect(keys(datadict)))
- for gj in gates
- specs[gj] = multitaper(datadict[gj], nw=5.0, k=8,
- aweight=true, shift=true, center=true)
- end
- specs
-end
-
-# For coherences. A very simple sum of squares. But honestly, these parameters
-# are hard enough to estimate even for maximum likelihood that you might be
-# better off just plugging in random values that don't make obviously wrong
-# curves.
-function cross_objective(parms, crosses, gates)
- out = 0.0
- len = length(crosses[(gates[1], gates[2])])
- fgd = collect(range(-0.5,0.5,length=len+1)[1:len])
- for gj in gates
- for gk in (gj+1):gates[end]
- model_cross = [Cw(f, gj*30, gk*30, parms) for f in fgd]
- out += sum(abs2, model_cross .- crosses[(gj,gk)])
- end
- end
- out
-end
-
-function marginal_objective(parms, specs, gates)
- out = 0.0
- len = length(specs[gates[1]])
- fgd = collect(range(-0.5,0.5,length=len+1)[1:len])
- for gj in gates
- model_spec = [Sf(f,gj*30.0,parms) + parms[24]^2 for f in fgd]
- out += sum(log.(model_spec) .+ specs[gj]./model_spec)
- end
- out
-end
-
-
-# So, the INIT here can't be truly naive, like just a bunch of ones or
-# something, because then things like \beta will never come through. But it CAN
-# be pretty naive, as you can see here. Here is some justification for the
-# things that aren't just one:
-#
-# (i) 13.0: this is a range parameter for the coherence, which is exponentiated.
-# They spatial locations are coded as 30,60,90,etc, so if I set that to one the
-# coherence functions will be zero for x \neq y. No optimization software could
-# overcome that. This value is chosen so that the coherence at the zero
-# frequency for adjacent gates is about 0.95.
-#
-# (ii, iii) -5.0, 0.0: again an exponentiated parameter for the coherence.
-# These are chosen to be what they are just to make the initial coherence
-# function have an even vaguely reasonable shape based on how the model is
-# paramterized. They're set to the same above and below the boundary layer.
-#
-# (iv) 500: boundary layer height. If this were one, every altitude would be
-# completely below the boundary layer and so the gradient of this objective with
-# respect to beta would be zero.
-#
-# (v) 0.0: a phase term. I just don't try to capture this at all here. Zero is a
-# more natural naive init than one for a parameter like this.
-#
-# (vi) 1.0e-3: a pretty small but nontrivial nugget. 1.0 is obviously too large for
-# this data and 0 is not a good init.
-#
-const INIT = vcat(ones(4), # ns scale parameters. Not using those here because of the extra low-frequency multiplier, and because maxlik is very good at picking these parms up.
- ones(4), # sdf parameters, will come through.
- 13.0, # first coherence parameter (scale-dependent). (i)
- -5.0, # second coherence parameter (ii)
- 0.0, # third coherence parameter (iii)
- 13.0, # first coherence parameter (scale-dependent). (i)
- -5.0, # second coherence parameter (ii)
- 0.0, # third coherence parameter (iii)
- 500.0, # beta height (iv)
- 1.0, # tau parameter
- ones(4), # extra low frequency power, will come through
- ones(2), # scale decay above boundary layer
- 0.0, # phase term, won't come through. (v)
- ones(2).*1.0e-3 # nuggets, won't come through. (vi)
- )
-
-
-# An incremental approach to the incremental generation of initializations for
-# the maximum likelihood estimation: fit marginal spectra first, then fit
-# cross-spectra. You need to use real derivatives here or else the generic first
-# order methods will drive the irrelevant parameters up to ridiculous values.
-# This is not particularly speedy on my box. But it does serve as nice
-# reproducible proof that you could start from an impossibly stupid
-# initialization and bootstrap yourself to a good MLE for the half-spectral model.
-function fit_day(day, timewin, gates, init=INIT)
- datadict = extract_day_slice(day, timewin, gates)
- out = Dict{String, Vector{Float64}}()
- out["init"] = deepcopy(init)
- for (objfn, spfn, nam) in ((marginal_objective, marginal_specs, "marginal"),
- (cross_objective, cross_specs, "coherence"))
- println("Fitting $nam spectra...")
- sp = spfn(datadict)
- obj = p -> objfn(p, sp, gates)
- gr = (g,p) -> ForwardDiff.gradient!(g, obj, p)
- hs = (h,p) -> ForwardDiff.hessian!(h, obj, p)
- ini = (nam == "marginal") ? out["init"] : out["marginal"]
- opts = Optim.Options(iterations=3_000, g_tol=1.0e-5,
- show_trace=true, show_every=10)
- est = optimize(obj, gr, hs, ini, NewtonTrustRegion(), opts)
- display(est)
- out[nam] = deepcopy(est.minimizer)
- end
- out
-end
-
-for (day, gates) in zip(("02", "03", "06", "20", "24", "28"),
- (7:23, 7:28, 7:23, 7:28, 7:28, 7:28))
- mle_init = fit_day(day, (14.02, 14.24), gates)
- println("Working with day $day...")
- outfilename = string("../../data/initialization/init_", day, "_14.csv")
- if haskey(mle_init, "coherence")
- writedlm(outfilename, mle_init["coherence"], ',')
- else
- writedlm(outfilename, mle_init["marginal"], ',')
- end
- println("Finished with day $day.\n\n\n")
-end
-
D example/paperscripts/model.jl => example/paperscripts/model.jl +0 -88
@@ 1,88 0,0 @@
-
-using SpecialFunctions, ForwardDiff, DelimitedFiles
-
-# For the derivatives of the kernel:
-function coord_deriv(f, p, j)
- ForwardDiff.derivative(z->f(vcat(p[1:(j-1)], z, p[(j+1):end])), p[j])
-end
-
-# The matern correlation for smoothness 1/2 and 3/2:
-mtn_cor_12(x) = exp(-abs(x))
-mtn_cor_32(x) = exp(-abs(x))*(1.0 + abs(x))
-function mtn_cor(nu, x)
- abs(x)<1.0e-8 && return 1.0
- out = (sqrt(2.0*nu)*x)^nu
- out *= besselk(nu, sqrt(2.0*nu)*x)
- out /= gamma(nu)*2.0^(nu-1.0)
- out
-end
-
-# The butterworth-like function. Several tweaks for numerical reasons to avoid
-# breaking AD.
-function butter(w,p)
- iszero(w) && return p[1] # tweak one: a branch for w=0.
- arg = abs2(w/p[2])^p[3]
- arg > 1.0e18 && return 0.0 # tweak two: a branch for the arg being too large.
- p[1]/(1.0 + arg)
-end
-
-function _logistic(x, p_0, p_1, shape, center)
- alpha = 1.0/(1.0 + exp(-shape*(x-center)))
- (1-alpha)*p_0 + alpha*p_1
-end
-
-function _logistic_multiple(x, shape, center, a::NTuple, b::NTuple)
- map(z->_logistic(x, z[1], z[2], shape, center), zip(a,b))
-end
-
-# Marginal spectral density:
-function Sf(w,x,p)
- (_i, _v) = _logistic_multiple(x, p[16], p[15], (exp(p[5]), p[6]), (exp(p[7]), p[8]))
- extra_parm = _logistic(x, exp(p[17]), exp(p[18]), p[16], p[15])
- extra = 1.0+butter(sinpi(w), (extra_parm, exp(p[20]), exp(p[19])))
- scale_x = x <= p[15] ? 1.0 : butter(x-p[15], (1.0, exp(p[21]), exp(p[22])))
- scale_x*extra*(_i*sinpi(w)^2 + 1)^(-_v-1/2)
-end
-
-function Cw(w,x,y,p)
- x == y && return 1.0
- (gp1x, gp1y) = (_logistic(z, exp(p[9]), exp(p[12]), p[16], p[15]) for z in (x,y))
- (gp2x, gp2y) = (_logistic(z, exp(p[10]), exp(p[13]), p[16], p[15]) for z in (x,y))
- (gp3x, gp3y) = (_logistic(z, exp(p[11]), exp(p[14]), p[16], p[15]) for z in (x,y))
- nux = x <= p[15] ? 0.5 : 0.75
- nuy = y <= p[15] ? 0.5 : 0.75
- gamx = butter(sinpi(w), (gp1x, gp2x, gp3x))
- gamy = butter(sinpi(w), (gp1y, gp2y, gp3y))
- gamxy = sqrt(0.5*(gamx+gamy))
- gamx_y = (gamx^0.25)*(gamy^0.25)
- return (gamx_y/gamxy)*mtn_cor(0.5*(nux+nuy), abs(x-y)/gamxy)
-end
-
-# Phase term:
-phase(w,x,y,p) = exp(2.0*pi*im*p[23]*sinpi(w)*(x-y))
-function phased1(w,x,y,p)
- 2.0*pi*im*sinpi(w)*(x-y)*phase(w,x,y,p)
-end
-
-# Whole integrand, and manual derivative with respect to phase parameter:
-integrand(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phase(w,x,y,p)
-integrand_dp1(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phased1(w,x,y,p)
-
-# Integrand with no phase for autodiff:
-integrand_np(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)
-
-# Derivatives of the kernel:
-dfuns =convert(Vector{Function},
- [(w,x,y,p)->coord_deriv(z->integrand_np(w,x,y,z), p, j)*phase(w,x,y,p)
- for j in 5:22])
-
-# do the phase functions derivative by hand because ForwardDiff doesn't
-# support complex numbers....
-push!(dfuns, integrand_dp1)
-
-# The nugget function, added as a post-call. This derivative function gets
-# added to the list of kernel derivatives in the space-time domain.
-nugfn(x,y,p) = Float64(x==y)*p[24]^2 + Float64(x[1]==y[1])*p[25]^2
-dnugfn1(x,y,p) = Float64(x==y)*2.0*p[24]
-dnugfn2(x,y,p) = Float64(x[1]==y[1])*2.0*p[25]
-
D example/paperscripts/model_untransformed.jl => example/paperscripts/model_untransformed.jl +0 -88
@@ 1,88 0,0 @@
-
-using SpecialFunctions, ForwardDiff, DelimitedFiles
-
-# For the derivatives of the kernel:
-function coord_deriv(f, p, j)
- ForwardDiff.derivative(z->f(vcat(p[1:(j-1)], z, p[(j+1):end])), p[j])
-end
-
-# The matern correlation for smoothness 1/2 and 3/2:
-mtn_cor_12(x) = exp(-abs(x))
-mtn_cor_32(x) = exp(-abs(x))*(1.0 + abs(x))
-function mtn_cor(nu, x)
- abs(x)<1.0e-8 && return 1.0
- out = (sqrt(2.0*nu)*x)^nu
- out *= besselk(nu, sqrt(2.0*nu)*x)
- out /= gamma(nu)*2.0^(nu-1.0)
- out
-end
-
-# The butterworth-like function. Several tweaks for numerical reasons to avoid
-# breaking AD.
-function butter(w,p)
- iszero(w) && return p[1] # tweak one: a branch for w=0.
- arg = abs2(w/p[2])^p[3]
- arg > 1.0e18 && return 0.0 # tweak two: a branch for the arg being too large.
- p[1]/(1.0 + arg)
-end
-
-function _logistic(x, p_0, p_1, shape, center)
- alpha = 1.0/(1.0 + exp(-shape*(x-center)))
- (1-alpha)*p_0 + alpha*p_1
-end
-
-function _logistic_multiple(x, shape, center, a::NTuple, b::NTuple)
- map(z->_logistic(x, z[1], z[2], shape, center), zip(a,b))
-end
-
-# Marginal spectral density:
-function Sf(w,x,p)
- (_i, _v) = _logistic_multiple(x, p[16], p[15], (exp(p[5]), p[6]), (exp(p[7]), p[8]))
- extra_parm = _logistic(x, p[17], p[18], p[16], p[15])
- extra = 1.0+butter(sinpi(w), (extra_parm, p[20], p[19]))
- scale_x = x <= p[15] ? 1.0 : butter(x-p[15], (1.0, p[21], p[22]))
- scale_x*extra*(_i*sinpi(w)^2 + 1)^(-_v-1/2)
-end
-
-function Cw(w,x,y,p)
- x == y && return 1.0
- (gp1x, gp1y) = (_logistic(z, exp(p[9]), exp(p[12]), p[16], p[15]) for z in (x,y))
- (gp2x, gp2y) = (_logistic(z, p[10], p[13], p[16], p[15]) for z in (x,y))
- (gp3x, gp3y) = (_logistic(z, p[11], p[14], p[16], p[15]) for z in (x,y))
- nux = x <= p[15] ? 0.5 : 0.75
- nuy = y <= p[15] ? 0.5 : 0.75
- gamx = butter(sinpi(w), (gp1x, gp2x, gp3x))
- gamy = butter(sinpi(w), (gp1y, gp2y, gp3y))
- gamxy = sqrt(0.5*(gamx+gamy))
- gamx_y = (gamx^0.25)*(gamy^0.25)
- return (gamx_y/gamxy)*mtn_cor(0.5*(nux+nuy), abs(x-y)/gamxy)
-end
-
-# Phase term:
-phase(w,x,y,p) = exp(2.0*pi*im*p[23]*sinpi(w)*(x-y))
-function phased1(w,x,y,p)
- 2.0*pi*im*sinpi(w)*(x-y)*phase(w,x,y,p)
-end
-
-# Whole integrand, and manual derivative with respect to phase parameter:
-integrand(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phase(w,x,y,p)
-integrand_dp1(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phased1(w,x,y,p)
-
-# Integrand with no phase for autodiff:
-integrand_np(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)
-
-# Derivatives of the kernel:
-dfuns =convert(Vector{Function},
- [(w,x,y,p)->coord_deriv(z->integrand_np(w,x,y,z), p, j)*phase(w,x,y,p)
- for j in 5:22])
-
-# do the phase functions derivative by hand because ForwardDiff doesn't
-# support complex numbers....
-push!(dfuns, integrand_dp1)
-
-# The nugget function, added as a post-call. This derivative function gets
-# added to the list of kernel derivatives in the space-time domain.
-nugfn(x,y,p) = Float64(x==y)*p[24]^2 + Float64(x[1]==y[1])*p[25]^2
-dnugfn1(x,y,p) = Float64(x==y)*2.0*p[24]
-dnugfn2(x,y,p) = Float64(x[1]==y[1])*2.0*p[25]
-
D example/paperscripts/nsscale.jl => example/paperscripts/nsscale.jl +0 -35
@@ 1,35 0,0 @@
-
-using LinearAlgebra
-
-mutable struct NonstationaryScale{T,B,F,N} <: Function
- locs::NTuple{N,T}
- ix::NTuple{N,Int64}
-end
-
-function nwts(pt, locs::NTuple{N,T})::NTuple{N, Float64} where{N,T}
- tmp = map(z->exp(-abs(pt-z)/50), locs)
- tmp ./ sum(tmp)
-end
-
-# Nonstationary scale function evaluate at space-time location
-function (NS::NonstationaryScale{T,false,0,N})(x, y, p) where{T,N}
- out = 1.0
- for _x in (x,y)
- (time, space) = _x
- wts = nwts(time, NS.locs)
- out *= sum(z->p[z[1]]*z[2], zip(NS.ix, wts))
- end
- out
-end
-
-# Efficient derivatives---but need the chain rule for two applications!
-function (NS::NonstationaryScale{T,true,J,N})(x, y, p) where{T,J,N}
- (wx, wy) = map(pt->nwts(pt[1], NS.locs), (x,y))
- # Evaluations for x and y:
- (nx, ny) = map(w->sum(z->p[z[1]]*z[2], zip(NS.ix, w)), (wx,wy))
- # Derivatives for x and y:
- (dx, dy) = map(w->w[J]*2.0*p[NS.ix[J]], (wx,wy))
- # Return chain rule output:
- nx*dy + dx*ny
-end
-
D example/paperscripts/untransformer.jl => example/paperscripts/untransformer.jl +0 -29
@@ 1,29 0,0 @@
-
-function transf(v)
- [v[1],
- v[2],
- v[3],
- v[4],
- v[5],
- v[6],
- v[7],
- v[8],
- v[9],
- exp(v[10]),
- exp(v[11]),
- v[12],
- exp(v[13]),
- exp(v[14]),
- v[15],
- v[16],
- exp(v[17]),
- exp(v[18]),
- exp(v[19]),
- exp(v[20]),
- exp(v[21]),
- exp(v[22]),
- v[23],
- v[24],
- v[25]]
-end
-
D example/paperscripts/writer.jl => example/paperscripts/writer.jl +0 -31
@@ 1,31 0,0 @@
-using DelimitedFiles
-
-function tablesave!(name, M::Matrix; rlabels=nothing, clabels=nothing)
- out = string.(M)
- if !isnothing(rlabels)
- out = hcat(string.(rlabels), out)
- end
- if !isnothing(clabels)
- is_row_nothing = isnothing(rlabels)
- rw1 = string.(clabels)
- if !is_row_nothing
- pushfirst!(rw1, "")
- end
- rw1[end] = rw1[end] * "\\\\"
- writedlm("header_"*name, reshape(rw1, 1, length(rw1)), '&')
- end
- out[:,end] .= out[:,end] .* repeat(["\\\\"], size(out,1))
- writedlm(name, out, '&')
-end
-
-function gnuplot_save_matrix!(name, M::Matrix{Float64}, row_pts, col_pts)
- out = zeros(size(M,1)+1, size(M,2)+1)
- out[1,2:end] .= col_pts
- out[2:end,1] .= row_pts
- out[2:end, 2:end] .= M
- writedlm(name, out, ',')
-end
-
-function gnuplot_save_vector!(name, M, pts)
- writedlm(name, hcat(pts, M), ',')
-end
A example/withvecchia.jl => example/withvecchia.jl +57 -0
@@ 0,0 1,57 @@
+
+using Vecchia, ForwardDiff
+include("basic_example.jl")
+
+# This is a struct that is itself callable. But the point here is that sometimes
+# the parameters will be of different types, so this wraps a Dict{Type,
+# HSKernel}. It makes it possible to squeeze out some redundant computations for
+# AD gradients/Hessians.
+const KWrapper = HSKernelADWrapper{Int64}(Dict{Type,HSKernel}())
+
+# See the README for Vecchia.jl for more details, but this function is the only
+# place in the Vecchia.jl codebase where the kernel is evaluated to assemble
+# matrix blocks. So if you write your own method for this here, you can use your
+# own exotic types for the kernel.
+#
+# Here, you see that there's some extra logic about making sure that we have a
+# kernel that accepts parameters of type eltype(params).
+function Vecchia.updatebuf!(buf, pts1, pts2,
+ kfn::HSKernelADWrapper,
+ params; skipltri=false)
+ # Check if the wrapper has a kernel of the right type:
+ if !haskey(kfn.kernels, eltype(params))
+ newkernel = HSKernel(params, xpts, integrand, timelen)
+ kfn.kernels[eltype(params)] = newkernel
+ kernel = newkernel
+ else
+ kernel = kfn.kernels[eltype(params)]
+ end
+ # Now check if the parameters agree. If they don't, rebuild the kernel:
+ if !HalfSpectral.check_hskernel_params(kernel, params)
+ newkernel = HSKernel(params, xpts, integrand, timelen)
+ kfn.kernels[eltype(params)] = newkernel
+ kernel = newkernel
+ end
+ # From here, we know that the kernel is the right type and that the parameters
+ # agree. So we can really blaze through the rest of this:
+ for k in 1:length(pts2)
+ ptk = pts2[k]
+ for j in 1:length(pts1)
+ buf[j,k] = kernel(pts1[j], ptk)
+ end
+ end
+ nothing
+end
+
+# The point of the fancy wrapper struct and the new method is that now the
+# kernel composes nicely with the approximate Vecchia likelihood, which gives
+# O(n) complexity approximations that are in general pretty accurate. Here's an
+# nll, gradient, and Hessian that you can now plug in to your favorite
+# optimizer. An expressive kernel, a linear-cost approximation, second order
+# optimization....what's not to love?
+const rdata = randn(length(pts))
+const cfg = Vecchia.kdtreeconfig(rdata, pts, 20, 3, KWrapper)
+nll(p) = Vecchia.nll(cfg, p)
+grad(p) = ForwardDiff.gradient(smart_nll, p)
+hess(p) = ForwardDiff.hessian(smart_nll, p)
+
M src/HalfSpectral.jl => src/HalfSpectral.jl +7 -114
@@ 1,125 1,18 @@
module HalfSpectral
- export timegrid, tftkernel
+ using LinearAlgebra
- using LinearAlgebra, Statistics, FFTW
+ export HSKernel, HSKernelADWrapper
- abstract type TimeGridResult end
+ include("structstypes.jl")
- struct ZeroFunction <: Function
- end
- (zrf::ZeroFunction)(x,y,p) = 0.0
+ include("radix2.jl")
- struct IdFunction <: Function
- end
- (idf::IdFunction)(x,p) = 1.0
+ include("utils.jl")
- struct TimeGrid <: TimeGridResult
- Tv :: Vector{Float64}
- Xt :: Vector{Int64}
- dt :: Float64
- tol :: Float64
- end
+ include("methods.jl")
- struct TimeGridFailure <: TimeGridResult
- Tv ::Vector{Float64}
- tol ::Float64
- end
-
- mutable struct TimeFFTKernel{T,F<:Function,N,G<:Function,H<:Function,M} <: Function
- V :: TimeGrid
- Fn :: F
- D :: Dict{Tuple{Int64, T, T}, Float64}
- p :: NTuple{N, Float64}
- X :: AbstractVector{T}
- addfn :: G
- mulfn :: H
- nobuild_ix :: NTuple{M, Int64}
- padsz::Int64
- longmemory::Bool
- end
-
- function timegrid(X; tol=0.075, warn_tol=2*tol)::TimeGridResult
- dt, dtmax = extrema(abs.(diff(X)))
- d_reldiff = rem(dtmax-dt, dt)/dt
- if d_reldiff <= warn_tol
- d_reldiff > tol && (@info "Only the acceptable tol was acheived.")
- Xd = (X .- X[1])./dt
- Xd_int = Int64.(round.(Xd))
- minimum(diff(Xd_int)) > 0 || throw(error("time grid construction failed."))
- return TimeGrid(X, Xd_int, dt, d_reldiff)
- end
- return TimeGridFailure(X, d_reldiff)
- end
-
- function tftkernel(integrand::F, tg::TimeGrid, X::AbstractVector{T}, p;
- addfn::G=ZeroFunction(),
- mulfn::H=IdFunction(),
- nob=NTuple{0, Int64}(),
- padsz::Int64=7,
- longmemory=false) where{T,F,G,H}
- dict = Dict{Tuple{Int64, T, T}, Float64}()
- lp = length(p)
- nnb = length(nob)
- out = TimeFFTKernel{T,F,lp,G,H,nnb}(tg, integrand, dict, tuple(p...),
- X, addfn, mulfn, nob, padsz, longmemory)
- build!(out, p)
- return out
- end
-
- function build!(K::TimeFFTKernel{T,F,N,G,H,M}, p) where{T,F,N,G,H,M}
- ftn = nextprod([2,3,5,7], K.padsz*K.V.Xt[end])
- (K.longmemory && iseven(ftn)) && (ftn += 1)
- wgd = collect(range(-0.5, 0.5, length=ftn+1)[1:ftn])
- if K.longmemory
- @assert !(in(0.0, wgd)) "debug, this should be impossible."
- end
- ix = sort(unique(map(z->abs(z[2]-z[1]), Iterators.product(K.V.Xt, K.V.Xt))))
- plan = plan_ifft(wgd)
- for (x,y) in Iterators.product(K.X, K.X)
- covxy = real(plan*(fftshift([K.Fn(w,x,y,p) for w in wgd])))[ix.+1]
- newv = Dict(zip(zip(ix, fill(x, length(ix)), fill(y, length(ix))), covxy))
- merge!(K.D, newv)
- end
- K.p = convert(NTuple{N, Float64}, Tuple(p))
- return nothing
- end
-
- # Avoid the runtime cost of checking the parameters AND the post-function:
- function (K::TimeFFTKernel{T,F,N,ZeroFunction,IdFunction,0})(x, y) where{T,F,N}
- opt1 = (abs(x[1]-y[1]), x[2], y[2])
- haskey(K.D, opt1) && return K.D[opt1]
- throw(error("No computed covariance for these two points: $x, $y"))
- end
-
- # Avoid the runtime cost of checking the parameters, but call the post
- # computation function:
- function (K::TimeFFTKernel{T,F,N,G,H,M})(x, y) where{T,F,N,G,H,M}
- opt1 = (abs(x[1]-y[1]), x[2], y[2])
- haskey(K.D, opt1) && return K.mulfn(x,y,K.p)*K.D[opt1] + K.addfn(x,y,K.p)
- throw(error("No computed covariance for these two points: $x, $y"))
- end
-
- # Safer checked version:
- function (K::TimeFFTKernel{T,F,N,M})(x, y, p) where{T,F,N,M}
- update_p_flag = false
- for (j, (Kpj, pj)) in enumerate(zip(K.p, p))
- if Kpj != pj
- if !in(j, K.nobuild_ix)
- build!(K, p)
- update_p_flag = false
- break
- else
- update_p_flag = true
- end
- end
- end
- if update_p_flag
- K.p = tuple(p...)
- end
- K(x, y)
- end
+ include("interp.jl")
end
-
A src/methods.jl => src/methods.jl +45 -0
@@ 0,0 1,45 @@
+
+# TODO (cg 2022/02/24 11:50): Make sure that this doesn't allocate. Also, get
+# rid of this if thing about x[2:end] being a singleton or not.
+function (K::HSKernel{T,Q})(x, y) where{T,Q}
+ (tx, ty) = (x[1], y[1])
+ if length(x) == 2
+ (_x, _y) = (x[2], y[2])
+ else
+ (_x, _y) = (x[2:end], y[2:end])
+ end
+ key = (_x, _y)
+ haskey(K.store, key) ||throw(error("No key for location pair $key."))
+ ix = round(Int64, abs(tx-ty)+1)
+ K.store[key][ix]
+end
+
+function (K::HSKernel{T,Q})(tdif::Int64, x, y) where{T,Q}
+ key = (x, y)
+ haskey(K.store, key) ||throw(error("No key for location pair $key."))
+ K.store[key][tdif+1]
+end
+
+function (K::HSKernel{T,Q})(tdif::Float64, x, y) where{T,Q}
+ key = (x, y)
+ haskey(K.store, key) ||throw(error("No key for location pair $key."))
+ K.store[key][round(Int64, tdif)+1]
+end
+
+# Hopefully in this format the compiler will understand that it should hoist the
+# check outside of a hot loop, for example.
+function (K::HSKernel{T,Q})(x, y, p) where{T,Q}
+ check_hskernel_params(K, p)
+ K(x,y)
+end
+
+function Base.display(K::HSKernel)
+ println("Half-spectral kernel with:")
+ println(" - parameters $(round.(K.params, digits=3))")
+end
+
+#function updatewrapper(W::HSKernelWrapper{Q}, hsk::HSKernel{T,Q}) where{T,Q}
+# W.kernels[T] = hsk
+#end
+
+
A src/radix2.jl => src/radix2.jl +83 -0
@@ 0,0 1,83 @@
+
+# Very midly changed derived product from:
+# https://github.com/JuliaApproximation/FastTransforms.jl.
+#
+# copyright Richard Mikael Slevinsky and other contributors
+# originally distributed under MIT "Expat" license
+
+function generic_fft_pow2!(x::Vector{T}) where{T}
+ n,big2=length(x),2one(T)
+ nn,j=n÷2,1
+ for i=1:2:n-1
+ if j>i
+ x[j], x[i] = x[i], x[j]
+ x[j+1], x[i+1] = x[i+1], x[j+1]
+ end
+ m = nn
+ while m ≥ 2 && j > m
+ j -= m
+ m = m÷2
+ end
+ j += m
+ end
+ logn = 2
+ while logn < n
+ θ=-big2/logn
+ wtemp = sinpi(θ/2)
+ wpr, wpi = -2wtemp^2, sinpi(θ)
+ wr, wi = one(T), zero(T)
+ for m=1:2:logn-1
+ for i=m:2logn:n
+ j=i+logn
+ mixr, mixi = wr*x[j]-wi*x[j+1], wr*x[j+1]+wi*x[j]
+ x[j], x[j+1] = x[i]-mixr, x[i+1]-mixi
+ x[i], x[i+1] = x[i]+mixr, x[i+1]+mixi
+ end
+ wr = (wtemp=wr)*wpr-wi*wpi+wr
+ wi = wi*wpr+wtemp*wpi+wi
+ end
+ logn = logn << 1
+ end
+ return x
+end
+
+function rradix2(x)
+ @assert isinteger(log2(length(x))) "This only applies to vectors of length 2^k."
+ _x = zeros(eltype(x), 2*length(x))
+ ix = 1
+ @inbounds @simd for j in 1:2:length(_x)
+ _x[j] = x[ix]
+ ix += 1
+ end
+ generic_fft_pow2!(_x)
+ _x[1:2:end]
+end
+
+function irradix2(x)
+ @assert isinteger(log2(length(x))) "This only applies to vectors of length 2^k."
+ _x = zeros(eltype(x), 2*length(x))
+ ix = 1
+ @inbounds @simd for j in 1:2:length(_x)
+ _x[j] = x[ix]
+ ix += 1
+ end
+ generic_fft_pow2!(_x)
+ _x[1:2:end]./length(x)
+end
+
+# General purpose.
+function radix2(x, iflag)
+ @assert isinteger(log2(length(x))) "This only applies to vectors of length 2^k."
+ _x = zeros(eltype(x), 2*length(x))
+ ix = 1
+ @inbounds @simd for j in 1:2:length(_x)
+ _x[j] = x[ix]
+ ix += 1
+ end
+ generic_fft_pow2!(_x)
+ if iflag
+ return complex.(_x[1:2:end], _x[2:2:end])./length(x)
+ else
+ return complex.(_x[1:2:end], _x[2:2:end])
+ end
+end
A src/structstypes.jl => src/structstypes.jl +35 -0
@@ 0,0 1,35 @@
+
+# (H)alf-(S)pectral Kernel. Defaults to long memory-compatible gridding.
+# Unlike the original implementation, I'm going to leave the nugget and
+# nonstationary scale outside this kernel, or if somebody wants they'll have to
+# put it into the integrand function.
+#
+# As another simplification, I'm just going to assume that the data put in is
+# equally spaced and with unit spacing and that dt=1. So all of that rescaling
+# is just going to go into the parameters.
+#
+# In general, this version is built with NOT updating in mind, instead just
+# building a new kernel.
+struct HSKernel{T,Q}
+ params::Vector{T}
+ store::Dict{Tuple{Q, Q}, Vector{T}}
+end
+
+# Integrand needs to be of the signature (f, x, y, params) ->
+# C_f(x,y)*sqrt(S_x(f)*S_y(f))*... and stuff.
+function HSKernel(params, locs, integrand, timelen; padsize=7, use_sinpi=false)
+ ftsz = nextprod([2], max(512, padsize*timelen))
+ sfgd = _fftshift(collect(range(-0.5, 0.5, length=ftsz+1)[1:ftsz]))
+ lprs = Iterators.product(locs, locs)
+ Dv = map(xy->_subdict(params, xy, sfgd, integrand, timelen, use_sinpi), lprs)
+ HSKernel(params, Dict(Dv))
+end
+
+struct HSKernelADWrapper{Q}
+ kernels::Dict{Type, HSKernel}
+end
+
+#function HSKernelADWrapper(Q::Type)
+# HSKernelWrapper{Q}(Dict{Type, HSKernel}())
+#end
+
A src/utils.jl => src/utils.jl +26 -0
@@ 0,0 1,26 @@
+
+_fftshift(v) = circshift(v, div(length(v), 2))
+_ifftshift(v) = circshift(v, div(length(v)+1, 2))
+
+function _subdict(params, xy, sfgd, integrand, timelen, use_sinpi)
+ (x,y) = xy
+ if use_sinpi
+ integrand_grid = [integrand(sinpi(f),x,y,params) for f in sfgd]
+ else
+ integrand_grid = [integrand(f,x,y,params) for f in sfgd]
+ end
+ covxy = irradix2(integrand_grid)[1:(timelen+1)]
+ xy => covxy
+end
+
+# Put this at the top of an assembly block and then you don't have to check the
+# parameters with every call.
+function check_hskernel_params(K::HSKernel{T,Q}, p) where{T,Q}
+ eltype(p) == T || return false
+ length(p) == length(K.params) || return false
+ @inbounds for j in eachindex(K.params)
+ K.params[j] == p[j] || return false
+ end
+ true
+end
+
A test/Manifest.toml => test/Manifest.toml +280 -0
@@ 0,0 1,280 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.1"
+manifest_format = "2.0"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "9489214b993cd42d17f44c36e359bf6a7c919abf"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.15.0"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.3"
+
+[[deps.CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
+[[deps.Compat]]
+deps = ["Dates", "LinearAlgebra", "UUIDs"]
+git-tree-sha1 = "924cdca592bc16f14d2f7006754a621735280b74"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.1.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.13"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DiffResults]]
+deps = ["StaticArrays"]
+git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.0.3"
+
+[[deps.DiffRules]]
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.11.0"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.FiniteDifferences]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Printf", "Random", "Richardson", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "0ee1275eb003b6fc7325cb14301665d1072abda1"
+uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
+version = "0.12.24"
+
+[[deps.ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "2f18915445b248731ec5db4e4a17e451020bf21e"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.30"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "b3364212fb5d870f724876ffcd34dd8ec6d98918"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.7"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.15"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.9"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NaNMath]]
+git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.0"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.4.2"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.Richardson]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "e03ca566bec93f8a3aeb059c8ef102f268a38949"
+uuid = "708f8203-808e-40c0-ba2d-98a6953ed40d"
+version = "1.4.0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "a9e798cae4867e3a41cae2dd9eb60c047f1212db"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.1.6"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "383a578bdf6e6721f480e749d503ebc8405a0b22"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.4.6"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
A test/Project.toml => test/Project.toml +7 -0
@@ 0,0 1,7 @@
+[deps]
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
A test/runtests.jl => test/runtests.jl +43 -0
@@ 0,0 1,43 @@
+
+using Test, HalfSpectral, StaticArrays, LinearAlgebra, QuadGK
+using FiniteDifferences, ForwardDiff
+
+# Simple SDF and a simple nontrivial coherence, which make an integrand:
+sdf(f,x,p) = p[1]*(p[2]^2 + f^2)^(-p[3] - 1/2)
+coh(f,x,y,p) = exp(-p[4]*(1+(p[5]*f)^2)*abs(x-y))
+integrand(f,x,y,p) = coh(f,x,y,p)*sqrt(sdf(f,x,p)*sdf(f,y,p))
+
+const pts = reduce(vcat, map(x->[@SVector [t, x] for t in 1:100], 1:5))
+const params = ones(5)
+const kernel = HalfSpectral.HSKernel(params, 1:5, integrand, 100)
+const x = randn(length(pts))
+
+function qform(p)
+ Kf = cholesky!(Symmetric([kernel(x,y) for x in pts, y in pts]))
+ sum(_x->_x^2, Kf.U'\x)
+end
+
+# test 1: integral accuracy.
+@testset "accuracy" begin
+ for (x1,x2) in ((1,1), (1,2), (3,5))
+ for k in 1:5:30
+ qk = quadgk(f->cos(2*pi*(k-1)*f)*integrand(f,x1,x2,params), -1/2, 1/2)[1]
+ @test isapprox(qk, kernel.store[(x1,x2)][k], rtol=1e-2)
+ end
+ end
+end
+
+# test 2: positive definite.
+@testset "positive definite" begin
+ K = Symmetric([kernel(x,y) for x in pts, y in pts])
+ @test isposdef(K)
+end
+
+# test3: Accurate AD:
+@testset "autodiff" begin
+ fdd = central_fdm(10,1)
+ fdg = FiniteDifferences.grad(fdd, qform, params)[1]
+ adg = ForwardDiff.gradient(qform, params)
+ @test isapprox(fdg, adg, atol=1e-8)
+end
+