[Pkg-octave-commit] [octave-divand] 01/05: Imported Upstream version 1.1.2
Rafael Laboissière
rlaboiss-guest at moszumanska.debian.org
Fri Sep 2 14:35:01 UTC 2016
This is an automated email from the git hooks/post-receive script.
rlaboiss-guest pushed a commit to branch master
in repository octave-divand.
commit 301f4a1ef443155e3591ba94b0ff91a4fe6ebc95
Author: Rafael Laboissiere <rafael at debian.org>
Date: Sat Aug 20 16:19:08 2016 -0300
Imported Upstream version 1.1.2
---
.hg_archival.txt | 5 +
CITATION | 22 +++
COPYING | 337 ++++++++++++++++++++++++++++++++++++
DESCRIPTION | 10 ++
INDEX | 124 +++++++++++++
NEWS | 21 +++
doc/gmd-7-225-2014.pdf | Bin 0 -> 1329711 bytes
inst/@CatBlockMat/CatBlockMat.m | 49 ++++++
inst/@CatBlockMat/append.m | 35 ++++
inst/@CatBlockMat/ctranspose.m | 33 ++++
inst/@CatBlockMat/double.m | 29 ++++
inst/@CatBlockMat/end.m | 24 +++
inst/@CatBlockMat/full.m | 24 +++
inst/@CatBlockMat/mldivide.m | 30 ++++
inst/@CatBlockMat/mtimes.m | 66 +++++++
inst/@CatBlockMat/size.m | 31 ++++
inst/@CatBlockMat/subsref.m | 41 +++++
inst/@CovarFun/CovarFun.m | 60 +++++++
inst/@CovarFun/mldivide.m | 48 +++++
inst/@CovarIS/CovarIS.m | 31 ++++
inst/@CovarIS/diag.m | 45 +++++
inst/@CovarIS/factorize.m | 41 +++++
inst/@CovarIS/full.m | 24 +++
inst/@CovarIS/inv.m | 23 +++
inst/@CovarIS/mldivide.m | 25 +++
inst/@CovarIS/mtimes.m | 37 ++++
inst/@CovarIS/subsref.m | 67 +++++++
inst/@CovarLanczos/CovarLanczos.m | 28 +++
inst/@CovarLanczos/diag.m | 30 ++++
inst/@CovarLanczos/mtimes.m | 25 +++
inst/@CovarParam/CovarParam.m | 69 ++++++++
inst/@CovarParam/diag.m | 24 +++
inst/@CovarParam/full.m | 25 +++
inst/@CovarParam/isscalar.m | 24 +++
inst/@CovarParam/isvector.m | 25 +++
inst/@CovarParam/mldivide.m | 47 +++++
inst/@CovarParam/mtimes.m | 39 +++++
inst/@CovarParam/mtimescorr.m | 40 +++++
inst/@CovarParam/size.m | 31 ++++
inst/@CovarParam/subsref.m | 35 ++++
inst/@CovarSMW/CovarSMW.m | 46 +++++
inst/@CovarSMW/diag.m | 24 +++
inst/@CovarSMW/full.m | 24 +++
inst/@CovarSMW/mldivide.m | 25 +++
inst/@CovarSMW/mtimes.m | 24 +++
inst/@CovarSMW/size.m | 31 ++++
inst/@DiagBlockMat/DiagBlockMat.m | 44 +++++
inst/@DiagBlockMat/append.m | 32 ++++
inst/@DiagBlockMat/diag.m | 30 ++++
inst/@DiagBlockMat/full.m | 24 +++
inst/@DiagBlockMat/inv.m | 28 +++
inst/@DiagBlockMat/mldivide.m | 43 +++++
inst/@DiagBlockMat/mtimes.m | 33 ++++
inst/@DiagBlockMat/size.m | 31 ++++
inst/@DiagBlockMat/sum.m | 16 ++
inst/@MatFun/MatFun.m | 39 +++++
inst/@MatFun/ctranspose.m | 25 +++
inst/@MatFun/double.m | 24 +++
inst/@MatFun/end.m | 23 +++
inst/@MatFun/full.m | 24 +++
inst/@MatFun/mtimes.m | 40 +++++
inst/@MatFun/size.m | 31 ++++
inst/cat_cell_array.m | 56 ++++++
inst/colordisp.m | 27 +++
inst/colormsg.m | 71 ++++++++
inst/conjugategradient.m | 198 +++++++++++++++++++++
inst/divand.m | 290 +++++++++++++++++++++++++++++++
inst/divand_addc.m | 41 +++++
inst/divand_background.m | 275 +++++++++++++++++++++++++++++
inst/divand_background_components.m | 80 +++++++++
inst/divand_constr_advec.m | 67 +++++++
inst/divand_diagnose.m | 47 +++++
inst/divand_eof_contraint.m | 68 ++++++++
inst/divand_error.m | 88 ++++++++++
inst/divand_factorize.m | 82 +++++++++
inst/divand_kernel.m | 84 +++++++++
inst/divand_laplacian.m | 119 +++++++++++++
inst/divand_metric.m | 56 ++++++
inst/divand_obs.m | 83 +++++++++
inst/divand_operators.m | 112 ++++++++++++
inst/divand_orthogonalize_modes.m | 48 +++++
inst/divand_pc_michol.m | 50 ++++++
inst/divand_pc_none.m | 30 ++++
inst/divand_pc_sqrtiB.m | 31 ++++
inst/divand_realistic_example.m | 99 +++++++++++
inst/divand_rms.m | 37 ++++
inst/divand_simple_example.m | 68 ++++++++
inst/divand_solve.m | 87 ++++++++++
inst/divand_solve_eof.m | 73 ++++++++
inst/interp_regular.m | 39 +++++
inst/localize_regular_grid.m | 67 +++++++
inst/localize_separable_grid.m | 90 ++++++++++
inst/mtimescorr.m | 46 +++++
inst/private/geomean.m | 19 ++
inst/sparse_diag.m | 26 +++
inst/sparse_diff.m | 81 +++++++++
inst/sparse_gradient.m | 59 +++++++
inst/sparse_interp.m | 127 ++++++++++++++
inst/sparse_pack.m | 35 ++++
inst/sparse_shift.m | 75 ++++++++
inst/sparse_stagger.m | 78 +++++++++
inst/sparse_trim.m | 61 +++++++
inst/statevector_init.m | 74 ++++++++
inst/statevector_pack.m | 75 ++++++++
inst/statevector_unpack.m | 64 +++++++
inst/test_1dvar.m | 30 ++++
inst/test_2dvar.m | 56 ++++++
inst/test_2dvar_adv.m | 42 +++++
inst/test_2dvar_check.m | 64 +++++++
inst/test_2dvar_check_correrr.m | 70 ++++++++
inst/test_2dvar_constrain.m | 84 +++++++++
inst/test_2dvar_cyclic.m | 97 +++++++++++
inst/test_2dvar_eof_check.m | 97 +++++++++++
inst/test_2dvar_lenxy.m | 54 ++++++
inst/test_2dvar_pcg.m | 65 +++++++
inst/test_2dvar_rellen.m | 47 +++++
inst/test_3dvar.m | 58 +++++++
inst/test_3dvar_large_stacked.m | 56 ++++++
inst/test_4dvar.m | 53 ++++++
inst/test_divand.m | 55 ++++++
inst/test_interp_1d.m | 37 ++++
inst/test_interp_2d.m | 39 +++++
inst/test_interp_regular.m | 83 +++++++++
inst/test_sparse_diff.m | 81 +++++++++
124 files changed, 6936 insertions(+)
diff --git a/.hg_archival.txt b/.hg_archival.txt
new file mode 100644
index 0000000..7a6a6ec
--- /dev/null
+++ b/.hg_archival.txt
@@ -0,0 +1,5 @@
+repo: 3a718b90cae30db215ea7c893d3677ec9de03f5f
+node: 3f827f9ac45150bed477e569746202c0a32d492b
+branch: default
+latesttag: null
+latesttagdistance: 5
diff --git a/CITATION b/CITATION
new file mode 100644
index 0000000..46cb256
--- /dev/null
+++ b/CITATION
@@ -0,0 +1,22 @@
+Please cite divand as follows if you use this software in a publication:
+
+ Barth, A., Beckers, J.-M., Troupin, C., Alvera-Azcárate, A., and
+ Vandenbulcke, L.: divand-1.0: n-dimensional variational data analysis
+ for ocean observations, Geosci. Model Dev., 7, 225-241,
+ doi:10.5194/gmd-7-225-2014, 2014.
+ URL http://www.geosci-model-dev.net/7/225/2014/
+
+A BibTeX entry for LaTeX users is:
+
+ @Article{gmd-7-225-2014,
+ AUTHOR = {Barth, A. and Beckers, J.-M. and Troupin, C. and Alvera-Azc\'arate, A. and Vandenbulcke, L.},
+ TITLE = {divand-1.0: n-dimensional variational data analysis for ocean observations},
+ JOURNAL = {Geoscientific Model Development},
+ VOLUME = {7},
+ YEAR = {2014},
+ NUMBER = {1},
+ PAGES = {225--241},
+ URL = {http://www.geosci-model-dev.net/7/225/2014/},
+ DOI = {10.5194/gmd-7-225-2014}
+ }
+
diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000..bdb226b
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,337 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The 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 Library 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, see <http://www.gnu.org/licenses/>.
+
+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 Library General
+Public License instead of this License.
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..11b8e52
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,10 @@
+Name: divand
+Version: 1.1.2
+Date: 2014-06-18
+Author: Alexander Barth <barth.alexander at gmail.com>
+Maintainer: Alexander Barth <barth.alexander at gmail.com>
+Title: divand
+Description: divand performs an n-dimensional variational analysis (interpolation) of arbitrarily located observations.
+Depends: octave (>= 3.4.0)
+License: GPLv2+
+Url: http://modb.oce.ulg.ac.be/mediawiki/index.php/divand, http://www.geosci-model-dev.net/7/225/2014/gmd-7-225-2014.html
diff --git a/INDEX b/INDEX
new file mode 100644
index 0000000..b1815bf
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,124 @@
+divand >> divand
+High-level functions
+ divand
+ divand_metric
+ divand_kernel
+Functions for advanced usage
+ divand_addc
+ divand_background_components
+ divand_background
+ divand_constr_advec
+ divand_diagnose
+ divand_eof_contraint
+ divand_error
+ divand_factorize
+ divand_laplacian
+ divand_obs
+ divand_operators
+ divand_orthogonalize_modes
+ divand_pc_none
+ divand_pc_michol
+ divand_pc_sqrtiB
+ divand_rms
+ divand_solve_eof
+ divand_solve
+Utility functions
+ cat_cell_array
+ colordisp
+ colormsg
+ conjugategradient
+ interp_regular
+ localize_regular_grid
+ localize_separable_grid
+ mtimescorr
+ statevector_init
+ statevector_pack
+ statevector_unpack
+Special matrices
+ @CatBlockMat/CatBlockMat
+ @CatBlockMat/append
+ @CatBlockMat/ctranspose
+ @CatBlockMat/double
+ @CatBlockMat/end
+ @CatBlockMat/full
+ @CatBlockMat/mldivide
+ @CatBlockMat/mtimes
+ @CatBlockMat/size
+ @CatBlockMat/subsref
+ @CovarFun/CovarFun
+ @CovarFun/mldivide
+ @CovarIS/CovarIS
+ @CovarIS/diag
+ @CovarIS/factorize
+ @CovarIS/full
+ @CovarIS/inv
+ @CovarIS/mldivide
+ @CovarIS/mtimes
+ @CovarIS/subsref
+ @CovarLanczos/CovarLanczos
+ @CovarLanczos/diag
+ @CovarLanczos/mtimes
+ @CovarParam/CovarParam
+ @CovarParam/diag
+ @CovarParam/full
+ @CovarParam/isscalar
+ @CovarParam/isvector
+ @CovarParam/mldivide
+ @CovarParam/mtimescorr
+ @CovarParam/mtimes
+ @CovarParam/size
+ @CovarParam/subsref
+ @CovarSMW/CovarSMW
+ @CovarSMW/diag
+ @CovarSMW/full
+ @CovarSMW/mldivide
+ @CovarSMW/mtimes
+ @CovarSMW/size
+ @DiagBlockMat/append
+ @DiagBlockMat/DiagBlockMat
+ @DiagBlockMat/diag
+ @DiagBlockMat/full
+ @DiagBlockMat/inv
+ @DiagBlockMat/mldivide
+ @DiagBlockMat/mtimes
+ @DiagBlockMat/size
+ @DiagBlockMat/sum
+ @MatFun/ctranspose
+ @MatFun/double
+ @MatFun/end
+ @MatFun/full
+ @MatFun/MatFun
+ @MatFun/mtimes
+ @MatFun/size
+Sparse matrix operations
+ sparse_diag
+ sparse_diff
+ sparse_gradient
+ sparse_interp
+ sparse_pack
+ sparse_shift
+ sparse_stagger
+ sparse_trim
+Test script
+ test_1dvar
+ test_2dvar_adv
+ test_2dvar_check_correrr
+ test_2dvar_check
+ test_2dvar_constrain
+ test_2dvar_cyclic
+ test_2dvar_eof_check
+ test_2dvar_lenxy
+ test_2dvar
+ test_2dvar_pcg
+ test_2dvar_rellen
+ test_3dvar_large_stacked
+ test_3dvar
+ test_4dvar
+ test_divand
+ test_interp_1d
+ test_interp_2d
+ test_interp_regular
+ test_sparse_diff
+Example
+ divand_simple_example
+ divand_realistic_example
\ No newline at end of file
diff --git a/NEWS b/NEWS
new file mode 100644
index 0000000..005aac0
--- /dev/null
+++ b/NEWS
@@ -0,0 +1,21 @@
+Summary of important user-visible changes for divand 1.1.2:
+------------------------------------------------------------
+
+ ** Include missing coefficient 2 in divand_metric.
+ ** Remove dependency to mapping package (which is currently unmaintained)
+
+Summary of important user-visible changes for divand 1.1.1:
+------------------------------------------------------------
+
+ ** Improved documentation
+
+Summary of important user-visible changes for divand 1.1:
+------------------------------------------------------------
+
+ ** new preconditioner: divand_pc_michol and divand_pc_sqrtiB
+ for iterative method
+
+Summary of important user-visible changes for divand 1.0:
+------------------------------------------------------------
+
+ ** Initial release
diff --git a/doc/gmd-7-225-2014.pdf b/doc/gmd-7-225-2014.pdf
new file mode 100644
index 0000000..3f536fc
Binary files /dev/null and b/doc/gmd-7-225-2014.pdf differ
diff --git a/inst/@CatBlockMat/CatBlockMat.m b/inst/@CatBlockMat/CatBlockMat.m
new file mode 100644
index 0000000..56cfe7b
--- /dev/null
+++ b/inst/@CatBlockMat/CatBlockMat.m
@@ -0,0 +1,49 @@
+% Create a matrix which represents the concatenation of a series of matrices.
+%
+% C = CatBlockMat(dim,X1,X2,...)
+%
+% Create a matrix C which represents the concatenation of a series of matrices
+% X1,X2,... over dimension dim. The matrices X1,X2,... can be special matrix
+% objects. The resulting matrix C behaves as a regular matrix.
+
+function retval = CatBlockMat(dim,varargin)
+
+self.d = 1:2; % only for matrices now
+self.dim = dim;
+self.B = varargin;
+self.N = length(varargin);
+self.i = zeros(self.N+1,1);
+
+self.sz = size(self.B{1});
+self.sz(dim) = 0;
+
+
+for i=1:self.N
+ sz = size(self.B{i});
+
+ if any(self.sz(self.d ~= dim) ~= sz(self.d ~= dim))
+ error('arguments dimensions are not consistent.');
+ end
+
+ self.i(i+1) = self.i(i) + sz(dim);
+end
+
+self.sz(dim) = self.i(end);
+
+retval = class(self,'CatBlockMat');
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/append.m b/inst/@CatBlockMat/append.m
new file mode 100644
index 0000000..c5d7658
--- /dev/null
+++ b/inst/@CatBlockMat/append.m
@@ -0,0 +1,35 @@
+% Append a matrix to a CatBlockMat matrix.
+%
+% A = append(B,M)
+%
+% Append the matrix M to a CatBlockMat matrix B.
+
+
+function self = append(self,M)
+
+sz = size(M);
+
+if any(self.sz(self.d ~= self.dim) ~= sz(self.d ~= self.dim))
+ error('arguments dimensions are not consistent.');
+end
+
+self.B{end+1} = M;
+self.N = self.N+1;
+
+
+self.i(end+1) = self.i(end) + sz(self.dim);
+self.sz(self.dim) = self.i(end);
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/ctranspose.m b/inst/@CatBlockMat/ctranspose.m
new file mode 100644
index 0000000..f29bed1
--- /dev/null
+++ b/inst/@CatBlockMat/ctranspose.m
@@ -0,0 +1,33 @@
+% Conjugate transpose of a CatBlockMat matrix.
+%
+% T = ctranspose(C)
+%
+% Returns the conjugate transpose of matrix C.
+
+
+function T = ctranspose(self)
+
+B = cell(1,self.N);
+
+for l=1:self.N
+ B{l} = self.B{l}';
+end
+
+d = 3-self.dim;
+T = CatBlockMat(d,B{:});
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/double.m b/inst/@CatBlockMat/double.m
new file mode 100644
index 0000000..2ad2bfc
--- /dev/null
+++ b/inst/@CatBlockMat/double.m
@@ -0,0 +1,29 @@
+% Convert a CatBlockMat matrix to a regular matrix.
+%
+% F = double(C)
+%
+% Convert the CatBlockMat C matrix to a regular matrix F.
+
+function F = double(self)
+
+b = {};
+for l=1:self.N
+ b{l} = double(self.B{l});
+end
+
+F = cat(self.dim,b{:});
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/end.m b/inst/@CatBlockMat/end.m
new file mode 100644
index 0000000..f21d55a
--- /dev/null
+++ b/inst/@CatBlockMat/end.m
@@ -0,0 +1,24 @@
+% Last index of a CatBlockMat matrix.
+%
+% e = end(C,k)
+%
+% Return the last index of the CatBlockMat C along dimension k.
+
+
+function e = end(self,k,n)
+
+e = self.sz(k);
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/full.m b/inst/@CatBlockMat/full.m
new file mode 100644
index 0000000..0246f55
--- /dev/null
+++ b/inst/@CatBlockMat/full.m
@@ -0,0 +1,24 @@
+% Convert a CatBlockMat matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the CatBlockMat C matrix to a regular matrix F.
+
+function F = full(self)
+
+F = cat(self.dim,self.B{:});
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/mldivide.m b/inst/@CatBlockMat/mldivide.m
new file mode 100644
index 0000000..4344c5b
--- /dev/null
+++ b/inst/@CatBlockMat/mldivide.m
@@ -0,0 +1,30 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of CatBlockMat A and another
+% matrix B.
+
+
+function p = mldivide(self,b)
+
+if self.sz(2) ~= size(b,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+p = full(self) \ b;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/mtimes.m b/inst/@CatBlockMat/mtimes.m
new file mode 100644
index 0000000..3d3b8cd
--- /dev/null
+++ b/inst/@CatBlockMat/mtimes.m
@@ -0,0 +1,66 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the CatBlockMat A and another matrix B.
+
+function p = mtimes(A,B)
+
+if size(A,2) ~= size(B,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+if isa(A,'CatBlockMat')
+
+ if A.dim == 1
+ p = zeros(size(A,1),size(B,2));
+
+ for l=1:A.N
+ i = A.i(l)+1:A.i(l+1);
+ p(i,:) = A.B{l} * B;
+ end
+ else
+ for l=1:A.N
+ i = (A.i(l)+1) :A.i(l+1);
+
+ % must use subsref explicetly
+ % subsref(B,S) = B(i,:)
+ % http://www.mathworks.com/support/solutions/en/data/1-18J9B/
+
+ S.type = '()';
+ S.subs = {i, ':'};
+ tmp = A.B{l} * subsref(B,S);
+
+ if l==1
+ p = tmp;
+ else
+ p = p + tmp;
+ end
+ end
+ end
+
+elseif isa(B,'CatBlockMat') && B.dim == 2
+ p = zeros(size(A,1),size(B,2));
+
+ for l=1:B.N
+ j = B.i(l)+1:B.i(l+1);
+ p(:,j) = A * B.B{l};
+ end
+else
+ warning('divand:fullmat','use full matrices for product');
+ p = full(A) * full(B);
+end
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/size.m b/inst/@CatBlockMat/size.m
new file mode 100644
index 0000000..341a061
--- /dev/null
+++ b/inst/@CatBlockMat/size.m
@@ -0,0 +1,31 @@
+% Size of a CatBlockMat matrix.
+%
+% sz = size(C)
+% sz = size(C,dim)
+%
+% Return the size of the CatBlockMat C matrix. If dim is specified, then only
+% the size along dimension dim is returned.
+
+function sz = size(self,dim);
+
+sz = self.sz;
+
+if nargin == 2
+ sz = sz(dim);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CatBlockMat/subsref.m b/inst/@CatBlockMat/subsref.m
new file mode 100644
index 0000000..0d0343b
--- /dev/null
+++ b/inst/@CatBlockMat/subsref.m
@@ -0,0 +1,41 @@
+% Subreference of a CatBlockMat matrix.
+%
+% B = subsref(C,idx)
+%
+% Perform the subscripted element selection operation according to
+% the subscript specified by idx.
+%
+% See also: subsref.
+
+function B = subsref(self,S)
+
+if strcmp(S.type,'()')
+ i = S.subs{self.dim};
+ l = find(self.i+1 == i(1));
+
+ itest = self.i(l)+1:self.i(l+1);
+
+ if any(i ~= itest) || ~all(strcmp(':',S.subs(self.d ~= self.dim)))
+ error('only whole blocks can be referenced');
+ end
+
+ B = self.B{l};
+else
+ error('Attempt to reference field of non-structure array');
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarFun/CovarFun.m b/inst/@CovarFun/CovarFun.m
new file mode 100644
index 0000000..34b1ae4
--- /dev/null
+++ b/inst/@CovarFun/CovarFun.m
@@ -0,0 +1,60 @@
+% Create a covariance matrix based on a function handler.
+%
+% CF = = CovarFun(n,fun)
+% CF = = CovarFun(n,fun,...)
+%
+% Create the covariance matrix CF of size n-by-n based on the function handler
+% fun. Optional key word parameters are:
+% 'isvec': 1 if function is vectorized and 0 (default) if not.
+% 'pc': function hander of the preconditioner for pcg
+% 'M1','M2','matit','tol': parameters for pcg.
+%
+% See also: pcg
+
+
+function retval = CovarFun(n,fun,varargin)
+
+isvec = 0;
+self.pc = [];
+self.M1 = [];
+self.M2 = [];
+self.tol = 1e-6;
+self.maxit = 100;
+
+for i=1:2:length(varargin)
+ if strcmp(varargin{i},'isvec')
+ isvec = varargin{i+1};
+ elseif strcmp(varargin{i},'pc')
+ self.pc = varargin{i+1};
+ elseif strcmp(varargin{i},'M1')
+ self.M1 = varargin{i+1};
+ elseif strcmp(varargin{i},'M2')
+ self.M2 = varargin{i+1};
+ elseif strcmp(varargin{i},'maxit')
+ self.maxit = varargin{i+1};
+ elseif strcmp(varargin{i},'tol')
+ self.tol = varargin{i+1};
+ else
+ error('CovarFun:param','unknown param');
+ end
+end
+
+self.fun = fun;
+
+retval = class(self,'CovarFun',MatFun([n,n],fun,fun,isvec));
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarFun/mldivide.m b/inst/@CovarFun/mldivide.m
new file mode 100644
index 0000000..1b38e6e
--- /dev/null
+++ b/inst/@CovarFun/mldivide.m
@@ -0,0 +1,48 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of CovarParam A and another
+% matrix B using pcg.
+
+function p = mldivide(self,b)
+
+if size(self,2) ~= size(b,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+if isempty(self.M1)
+ args = {self.pc};
+else
+ args = {self.M1,self.M2};
+end
+
+for i= 1:size(b,2)
+
+ [p(:,i),flag,relres,iter] = pcg(self.fun, b(:,i), self.tol,...
+ self.maxit,args{:});
+
+ %whos b
+ self.pc
+ fprintf('iter %d %g \n',iter,relres);
+ if (flag ~= 0)
+ error('CovarFun:pcg', ['Preconditioned conjugate gradients method'...
+ ' did not converge %d %g %g'],flag,relres,iter);
+ end
+
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/CovarIS.m b/inst/@CovarIS/CovarIS.m
new file mode 100644
index 0000000..d715d05
--- /dev/null
+++ b/inst/@CovarIS/CovarIS.m
@@ -0,0 +1,31 @@
+% Covariance matrix with a sparse inverse matrix.
+%
+% C = CovarIS(IS)
+%
+% Create the covariance matrix C whose inverse is the sparse matrix IS.
+% To accelerate matrix product C*x, IS can be factorized.
+
+function retval = CovarIS(IS);
+
+self.IS = IS;
+self.f = 0;
+self.RP = [];
+self.q = [];
+self.PP = [];
+
+retval = class(self,'CovarIS');
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/diag.m b/inst/@CovarIS/diag.m
new file mode 100644
index 0000000..0eacb6c
--- /dev/null
+++ b/inst/@CovarIS/diag.m
@@ -0,0 +1,45 @@
+% Diagonal elements.
+%
+% d = diag(C)
+%
+% Return the diagonal elements of the CovarIS matrix C.
+
+function d = diag(self)
+
+N = size(self.IS,1);
+d = zeros(N,1);
+
+
+if self.f
+ t0 = cputime;
+
+ for i=1:N
+ % show progress every 5 seconds
+ if (t0 + 5 < cputime)
+ t0 = cputime;
+ fprintf('%d/%d\n',i,N);
+ end
+
+ d(i) = self.PP(i,:) * (self.RP \ (self.RP' \ self.PP(i,:)'));
+ end
+else
+ disp('factorize CovarIS');
+ self_fact = factorize(self);
+ d = diag(self_fact);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/factorize.m b/inst/@CovarIS/factorize.m
new file mode 100644
index 0000000..36b01f9
--- /dev/null
+++ b/inst/@CovarIS/factorize.m
@@ -0,0 +1,41 @@
+% Factorize matrix.
+%
+% self = factorize(self)
+%
+% Factorize matrix self using the Cholesky decomposition.
+
+function self = factorize(self)
+
+
+if issparse(self.IS)
+ [self.RP, self.q, self.PP] = chol (self.IS);
+
+ if self.q ~= 0
+ error('factoziation failed (matrix is not positive definite)');
+ end
+
+
+ %length(find(self.PP))/numel(self.PP)
+ %length(find(self.RP))/numel(self.RP)
+
+else
+ self.RP = chol(self.IS);
+ self.PP = speye(size(self.IS));
+end
+
+self.f = 1;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/full.m b/inst/@CovarIS/full.m
new file mode 100644
index 0000000..d020374
--- /dev/null
+++ b/inst/@CovarIS/full.m
@@ -0,0 +1,24 @@
+% Convert a CovarIS matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the CovarIS C matrix to a regular matrix F.
+
+function M = full(self)
+
+M = inv(self.IS);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/inv.m b/inst/@CovarIS/inv.m
new file mode 100644
index 0000000..4a11f37
--- /dev/null
+++ b/inst/@CovarIS/inv.m
@@ -0,0 +1,23 @@
+% Inverse of a CovarIS matrix.
+%
+% M = inv(C)
+%
+% Return the inverse of the CovarIS matrix C.
+
+function M = inv(self)
+
+M = self.IS;
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/mldivide.m b/inst/@CovarIS/mldivide.m
new file mode 100644
index 0000000..2a75e81
--- /dev/null
+++ b/inst/@CovarIS/mldivide.m
@@ -0,0 +1,25 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of CovarIS A and another
+% matrix B.
+
+function q = mldivide(self,b)
+
+q = self.IS * b;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/mtimes.m b/inst/@CovarIS/mtimes.m
new file mode 100644
index 0000000..fc5a9d2
--- /dev/null
+++ b/inst/@CovarIS/mtimes.m
@@ -0,0 +1,37 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the CovarIS matrix A and another matrix B.
+
+function p = mtimes(self,b)
+
+%keyboard
+if self.f
+ %'fact'
+ p = self.PP * (self.RP \ (self.RP' \ (self.PP' * b)));
+else
+% 'direct'
+% tic
+% A = (self.IS + self.IS')/2;
+% p = A \ b;
+% toc
+% tic
+ p = self.IS \ b;
+% toc
+% 'end'
+end
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarIS/subsref.m b/inst/@CovarIS/subsref.m
new file mode 100644
index 0000000..336c6db
--- /dev/null
+++ b/inst/@CovarIS/subsref.m
@@ -0,0 +1,67 @@
+% Subreference of a CovarIS matrix.
+%
+% B = subsref(C,idx)
+%
+% Perform the subscripted element selection operation according to
+% the subscript specified by idx.
+%
+% See also: subsref.
+
+function x = subsref(self,idx)
+
+assert(length(idx) == 1)
+assert(strcmp(idx.type,'()'))
+N = size(self.IS,1);
+
+i = idx.subs{1};
+
+if strcmp(i,':')
+ i = 1:N;
+end
+
+j = idx.subs{2};
+
+if strcmp(j,':')
+ j = 1:N;
+end
+
+x = zeros(length(i),length(j));
+
+if self.f
+ t0 = cputime;
+
+ k = 0;
+ for jj=1:length(j)
+ for ii=1:length(i)
+
+ % show progress every 5 seconds
+ if (t0 + 5 < cputime)
+ t0 = cputime;
+ fprintf('%d/%d\n',k,numel(x));
+ end
+
+ x(ii,jj) = self.PP(i(ii),:) * (self.RP \ (self.RP' \ self.PP(j(jj),:)'));
+ k = k+1;
+ end
+ end
+else
+ disp('factorize CovarIS');
+ self_fact = factorize(self);
+ x = subsref(self_fact,idx);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarLanczos/CovarLanczos.m b/inst/@CovarLanczos/CovarLanczos.m
new file mode 100644
index 0000000..db5e7b9
--- /dev/null
+++ b/inst/@CovarLanczos/CovarLanczos.m
@@ -0,0 +1,28 @@
+% Covariance matrix reconstructed from Lanczos vectors.
+%
+% A = CovarLanczos(Q,T)
+%
+% Return a covariance matrix, factorized into A = Q * T * Q'
+% where Q'*Q = I and T is a tridiagonal matrix.
+
+function retval = CovarLanczos(Q,T)
+
+self.Q = Q;
+self.T = T;
+
+retval = class(self,'CovarLanczos');
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarLanczos/diag.m b/inst/@CovarLanczos/diag.m
new file mode 100644
index 0000000..b44bf50
--- /dev/null
+++ b/inst/@CovarLanczos/diag.m
@@ -0,0 +1,30 @@
+% Diagonal elements.
+%
+% d = diag(C)
+%
+% Return the diagonal elements of the CovarLanczos matrix C.
+
+function d = diag(self)
+
+N = size(self.Q,1);
+d = zeros(N,1);
+
+for i=1:N
+ d(i) = self.Q(i,:) * (self.T \ self.Q(i,:)');
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarLanczos/mtimes.m b/inst/@CovarLanczos/mtimes.m
new file mode 100644
index 0000000..32a0188
--- /dev/null
+++ b/inst/@CovarLanczos/mtimes.m
@@ -0,0 +1,25 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the CovarLanczos matrix A and another matrix
+% B.
+
+function p = mtimes(self,b)
+
+p = self.Q * (self.T \ (self.Q' * b));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/CovarParam.m b/inst/@CovarParam/CovarParam.m
new file mode 100644
index 0000000..b2f1d50
--- /dev/null
+++ b/inst/@CovarParam/CovarParam.m
@@ -0,0 +1,69 @@
+% Parametric covariance matrix with a guassian kernel.
+%
+% C = CovarParam({x,y,...},variance,len)
+%
+% Creates a parametric covariance matrix C with a guassian
+% kernel, a correlation length of len and a variance of
+% variance. The first parameter is a cell array with the coordinates.
+
+function retval = CovarParam(xi,variance,len,varargin)
+
+xi = squeeze(cat_cell_array(xi));
+m = size(xi,1);
+
+
+% fun = @(xi1,xi2) exp(- sum( (xi1 - xi2).^2 , 2) / len^2);
+%
+% tic
+% for j=1:m
+% for i=1:m
+% R2(i,j) = fun(xi(i,:),xi(j,:));
+% end
+% end
+% toc
+%
+% maxdiff(R,R2);
+
+self.tol = 1e-6;
+self.maxit = 100;
+self.var_add = 1;
+
+for i=1:2:length(varargin)
+ if strcmp(varargin{i},'tolerance')
+ self.tol = varargin{i+1};
+ elseif strcmp(varargin{i},'maxit')
+ self.maxit = varargin{i+1};
+ elseif strcmp(varargin{i},'fraccorr')
+ self.var_add = varargin{i+1};
+ end
+end
+
+self.m = m;
+self.xi = xi;
+self.len = len;
+
+if isscalar(variance)
+ variance = ones(self.m,1) * variance;
+end
+
+self.variance = variance;
+self.std = sqrt(variance);
+self.S = sparse_diag(self.std);
+
+retval = class(self,'CovarParam');
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/diag.m b/inst/@CovarParam/diag.m
new file mode 100644
index 0000000..df115d8
--- /dev/null
+++ b/inst/@CovarParam/diag.m
@@ -0,0 +1,24 @@
+% Diagonal elements.
+%
+% d = diag(C)
+%
+% Return the diagonal elements of the CovarParam matrix C.
+
+function d = diag(self);
+
+d = self.variance;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/full.m b/inst/@CovarParam/full.m
new file mode 100644
index 0000000..c6eda6d
--- /dev/null
+++ b/inst/@CovarParam/full.m
@@ -0,0 +1,25 @@
+% Convert a CovarParam matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the CovarParam C matrix to a regular matrix F.
+
+
+function F = full(self)
+
+F = self * eye(size(self));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/isscalar.m b/inst/@CovarParam/isscalar.m
new file mode 100644
index 0000000..78407bc
--- /dev/null
+++ b/inst/@CovarParam/isscalar.m
@@ -0,0 +1,24 @@
+% Return true if argument is a scalar.
+%
+% b = isscalar(X)
+%
+% Return true if X is a scalar.
+
+function b = isscalar(self)
+
+b = isequal(size(self),[1 1]);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/isvector.m b/inst/@CovarParam/isvector.m
new file mode 100644
index 0000000..34fb7df
--- /dev/null
+++ b/inst/@CovarParam/isvector.m
@@ -0,0 +1,25 @@
+% Return true if argument is a vector.
+%
+% b = isscalar(X)
+%
+% Return true if X is a vector. A vector is a N-by-1 or 1-by-N matrix.
+
+function b = isvector(self)
+
+sz = size(self);
+b = length(sz) == 2 && (sz(1) == 1 || sz(2) == 1);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/mldivide.m b/inst/@CovarParam/mldivide.m
new file mode 100644
index 0000000..67917e1
--- /dev/null
+++ b/inst/@CovarParam/mldivide.m
@@ -0,0 +1,47 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of CovarParam matrix A and
+% another matrix B.
+
+function q = mldivide(self,b)
+
+q = zeros(size(b));
+b = self.S \ b;
+
+for i=1:size(b,2)
+ M = [];
+ %M = @(x) ifft(x);
+
+ if 1
+ [q(:,i),flag,relres,iter] = pcg(@(x) mtimescorr(self,x), b(:,i),self.tol,self.maxit,M);
+ else
+
+ [q(:,i),flag,relres,iter] = pcg(@(x) fft(mtimescorr(self,ifft(x))), fft(b(:,i)),self.tol,self.maxit,M);
+ q(:,i) = real(ifft(q(:,i)));
+ end
+
+
+ if flag > 0
+ relres,iter
+ warning('divand:pcgFlag','pcg flag %d',flag);
+ end
+end
+
+q = self.S \ q;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/mtimes.m b/inst/@CovarParam/mtimes.m
new file mode 100644
index 0000000..1b0d9a7
--- /dev/null
+++ b/inst/@CovarParam/mtimes.m
@@ -0,0 +1,39 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the CovarParam matrix A and another matrix
+% B.
+
+function p = mtimes(self,b)
+
+m = self.m;
+p = zeros(size(b));
+
+% for j=1:m
+% for i=1:m
+% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 );
+% p(i) = p(i) + exp(-d2 / self.len^2) * b(j);
+% end
+% end
+%
+
+b = self.S*b;
+p = mtimescorr(self,b);
+p = self.S*p;
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/mtimescorr.m b/inst/@CovarParam/mtimescorr.m
new file mode 100644
index 0000000..9e227bc
--- /dev/null
+++ b/inst/@CovarParam/mtimescorr.m
@@ -0,0 +1,40 @@
+% Matrix product with correlation matrix.
+%
+% p = mtimescorr(A,B)
+%
+% Return the matrix product between the CovarParam matrix A (with unit variance)
+% and another matrix B.
+
+
+function p = mtimescorr(self,b)
+
+m = self.m;
+p = zeros(size(b));
+
+% for j=1:m
+% for i=1:m
+% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 );
+% p(i) = p(i) + exp(-d2 / self.len^2) * b(j);
+% end
+% end
+%
+
+p = mtimescorr(self.xi,self.len,b);
+%p = mtimescorr_f(self.xi,self.len,b);
+p = self.var_add * p + (1-self.var_add) * b;
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/size.m b/inst/@CovarParam/size.m
new file mode 100644
index 0000000..b960ce1
--- /dev/null
+++ b/inst/@CovarParam/size.m
@@ -0,0 +1,31 @@
+% Size of a CovarParam matrix.
+%
+% sz = size(C)
+% sz = size(C,dim)
+%
+% Return the size of the CovarParam C matrix. If dim is specified, then only
+% the size along dimension dim is returned.
+
+function sz = size(self,dim);
+
+sz = [self.m, self.m];
+
+if nargin == 2
+ sz = sz(dim);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarParam/subsref.m b/inst/@CovarParam/subsref.m
new file mode 100644
index 0000000..96c8f46
--- /dev/null
+++ b/inst/@CovarParam/subsref.m
@@ -0,0 +1,35 @@
+% Subreference of a CovarParam matrix.
+%
+% B = subsref(C,idx)
+%
+% Perform the subscripted element selection operation according to
+% the subscript specified by idx.
+%
+% See also: subsref.
+
+function C = subsref(self,idx)
+
+assert(strcmp(idx.type,'()'))
+
+i = idx.subs{1};
+j = idx.subs{2};
+
+assert(isequal(i,j));
+
+
+C = CovarParam(self.xi(i,:),self.variance(i),self.len);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/CovarSMW.m b/inst/@CovarSMW/CovarSMW.m
new file mode 100644
index 0000000..5584d2c
--- /dev/null
+++ b/inst/@CovarSMW/CovarSMW.m
@@ -0,0 +1,46 @@
+% Covariance matrix than can be inverted using the Sherman-Morrison-Woodbury
+% forumla.
+%
+% CSMW = CovarSMW(C,B)
+%
+% Create the covariance matrix CSMW = C + B*B' than can be inverted efficiently
+% using the Sherman-Morrison-Woodbury formula, if size(B,2) is much smaller than
+% size(C,1):
+%
+% inv(C + B*B') = inv(C) - inv(C)*B*inv(B'*inv(C)*B + I)*B'*inv(C)
+%
+% The symetric matrix C should implement the methods: size, diag, mtimes,
+% mldivde and the matrix B should implement the methods: size, tranpose, mtimes,
+% sum.
+
+function retval = CovarSMW(C,B)
+
+if size(C,1) ~= size(C,2)
+ error('C should be square matrix');
+end
+
+if size(C,1) ~= size(B,1)
+ error('size of C and B are not compatible');
+end
+
+self.C = C;
+self.B = B;
+self.D = inv(B' * (C \ B) + eye(size(B,2)));
+
+retval = class(self,'CovarSMW');
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/diag.m b/inst/@CovarSMW/diag.m
new file mode 100644
index 0000000..a58a148
--- /dev/null
+++ b/inst/@CovarSMW/diag.m
@@ -0,0 +1,24 @@
+% Diagonal elements.
+%
+% d = diag(C)
+%
+% Return the diagonal elements of the CovarSWM matrix C.
+
+function d = diag(self);
+
+d = diag(self.C) + sum(self.B.^2,2);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/full.m b/inst/@CovarSMW/full.m
new file mode 100644
index 0000000..2083c27
--- /dev/null
+++ b/inst/@CovarSMW/full.m
@@ -0,0 +1,24 @@
+% Convert a CovarSWM matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the CovarSWM C matrix to a regular matrix F.
+
+function F = full(self)
+
+F = self.C + self.B * self.B';
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/mldivide.m b/inst/@CovarSMW/mldivide.m
new file mode 100644
index 0000000..c95c537
--- /dev/null
+++ b/inst/@CovarSMW/mldivide.m
@@ -0,0 +1,25 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of CovarSWM A and another
+% matrix B.
+
+function Q = mldivide(self,A)
+
+Q = self.C \ A - self.C \ (self.B * (self.D * (self.B' * (self.C \ A))));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/mtimes.m b/inst/@CovarSMW/mtimes.m
new file mode 100644
index 0000000..c87a349
--- /dev/null
+++ b/inst/@CovarSMW/mtimes.m
@@ -0,0 +1,24 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the CovarSWM matrix A and another matrix B.
+
+function P = mtimes(self,A)
+
+P = self.C*A + self.B * (self.B'*A);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@CovarSMW/size.m b/inst/@CovarSMW/size.m
new file mode 100644
index 0000000..c40a239
--- /dev/null
+++ b/inst/@CovarSMW/size.m
@@ -0,0 +1,31 @@
+% Size of a CovarSWM matrix.
+%
+% sz = size(C)
+% sz = size(C,dim)
+%
+% Return the size of the CovarSWM C matrix. If dim is specified, then only
+% the size along dimension dim is returned.
+
+function sz = size(self,dim);
+
+sz = size(self.C);
+
+if nargin == 2
+ sz = sz(dim);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/DiagBlockMat.m b/inst/@DiagBlockMat/DiagBlockMat.m
new file mode 100644
index 0000000..74df161
--- /dev/null
+++ b/inst/@DiagBlockMat/DiagBlockMat.m
@@ -0,0 +1,44 @@
+% Block diagonal matrix.
+%
+% DB = DiagBlockMat(A1,A2,...,An)
+%
+% Create the block diagonal matrix DB from the matrices A1, A2,...,An
+%
+% / A1 0 0 ... 0 \
+% | 0 A2 0 : |
+% DB = | 0 0 A3 |
+% | : . 0 |
+% \ 0 ... 0 An /
+
+function retval = DiagBlockMat(varargin)
+
+self.B = varargin;
+self.N = length(varargin);
+self.i = zeros(self.N+1,1);
+self.j = zeros(self.N+1,1);
+
+for i=1:self.N
+ sz = size(self.B{i});
+ self.i(i+1) = self.i(i) + sz(1);
+ self.j(i+1) = self.j(i) + sz(2);
+end
+
+self.sz = [self.i(end), self.j(end)];
+
+retval = class(self,'DiagBlockMat');
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/append.m b/inst/@DiagBlockMat/append.m
new file mode 100644
index 0000000..cf8c661
--- /dev/null
+++ b/inst/@DiagBlockMat/append.m
@@ -0,0 +1,32 @@
+% Append a matrix to a DiagBlockMat matrix.
+%
+% A = append(B,M)
+%
+% Append the matrix M to a DiagBlockMat matrix B along its diagonal.
+
+function self = append(self,M)
+
+self.B{end+1} = M;
+self.N = self.N+1;
+
+sz = size(M);
+
+self.i(end+1) = self.i(end) + sz(1);
+self.j(end+1) = self.j(end) + sz(2);
+
+self.sz = [self.i(end), self.j(end)];
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/diag.m b/inst/@DiagBlockMat/diag.m
new file mode 100644
index 0000000..e87f3d0
--- /dev/null
+++ b/inst/@DiagBlockMat/diag.m
@@ -0,0 +1,30 @@
+% Diagonal elements.
+%
+% d = diag(C)
+%
+% Return the diagonal elements of the matrix C.
+
+function d = diag(self)
+
+t = cell(1,self.N);
+
+for l=1:self.N
+ t{l} = diag(self.B{l});
+end
+
+d = cat(1,t{:});
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/full.m b/inst/@DiagBlockMat/full.m
new file mode 100644
index 0000000..8b40c16
--- /dev/null
+++ b/inst/@DiagBlockMat/full.m
@@ -0,0 +1,24 @@
+% Convert a DiagBlockMat matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the DiagBlockMat C matrix to a regular matrix F.
+
+function F = full(self)
+
+F = blkdiag(self.B{:});
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/inv.m b/inst/@DiagBlockMat/inv.m
new file mode 100644
index 0000000..9e1bd01
--- /dev/null
+++ b/inst/@DiagBlockMat/inv.m
@@ -0,0 +1,28 @@
+% Inverse of a DiagBlockMat matrix.
+%
+% M = inv(C)
+%
+% Return the inverse of the DiagBlockMat matrix C.
+
+function iA = inv(A)
+
+iA = A;
+
+for l=1:A.N
+ iA.B{l} = inv(A.B{l});
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/mldivide.m b/inst/@DiagBlockMat/mldivide.m
new file mode 100644
index 0000000..5c9f97c
--- /dev/null
+++ b/inst/@DiagBlockMat/mldivide.m
@@ -0,0 +1,43 @@
+% Matrix left division A \ B.
+%
+% p = mldivide(A,B)
+%
+% Return the matrix product between the inverse of DiagBlockMat A and another
+% matrix B
+
+function p = mldivide(self,b)
+
+if size(self,2) ~= size(b,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+
+if isa(b,'CatBlockMat')
+ B = {};
+ for l=1:self.N
+ B{l} = self.B{l} \ b(self.j(l)+1:self.j(l+1),:);
+ end
+
+ p = CatBlockMat(1,B{:});
+else
+
+ p = zeros(size(self,1),size(b,2));
+
+ for l=1:self.N
+ p(self.i(l)+1:self.i(l+1),:) = self.B{l} \ b(self.j(l)+1:self.j(l+1),:);
+ end
+end
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/mtimes.m b/inst/@DiagBlockMat/mtimes.m
new file mode 100644
index 0000000..b50ecfc
--- /dev/null
+++ b/inst/@DiagBlockMat/mtimes.m
@@ -0,0 +1,33 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the DiagBlockMat matrix A and another matrix
+% B.
+
+function p = mtimes(self,b)
+
+if size(self,2) ~= size(b,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+p = zeros(size(self,1),size(b,2));
+
+for l=1:self.N
+ p(self.i(l)+1:self.i(l+1),:) = self.B{l} * b(self.j(l)+1:self.j(l+1),:);
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/size.m b/inst/@DiagBlockMat/size.m
new file mode 100644
index 0000000..37b3417
--- /dev/null
+++ b/inst/@DiagBlockMat/size.m
@@ -0,0 +1,31 @@
+% Size of a DiagBlockMat matrix.
+%
+% sz = size(C)
+% sz = size(C,dim)
+%
+% Return the size of the DiagBlockMat C matrix. If dim is specified, then only
+% the size along dimension dim is returned.
+
+function sz = size(self,dim);
+
+sz = self.sz;
+
+if nargin == 2
+ sz = sz(dim);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@DiagBlockMat/sum.m b/inst/@DiagBlockMat/sum.m
new file mode 100644
index 0000000..1165ba9
--- /dev/null
+++ b/inst/@DiagBlockMat/sum.m
@@ -0,0 +1,16 @@
+% Compute the sum.
+%
+% S = sum (X, DIM)
+%
+% Compute the sum along dimension DIM.
+
+
+function d = sum(self,i)
+
+t = cell(1,self.N);
+
+for l=1:self.N
+ t{l} = sum(self.B{l},i);
+end
+
+d = cat(3-i,t{:});
\ No newline at end of file
diff --git a/inst/@MatFun/MatFun.m b/inst/@MatFun/MatFun.m
new file mode 100644
index 0000000..ce77744
--- /dev/null
+++ b/inst/@MatFun/MatFun.m
@@ -0,0 +1,39 @@
+% Matrix operator object based on a function handel.
+%
+% MF = MatFun(sz,fun,funt,isvec)
+%
+% Create a operator of size sz based on the linear function fun and its adjoint
+% funt. For a matrix X of approriate size:
+%
+% MF*X is fun(X) and MF'*X is funt(X)
+%
+% If isvec is 0 (default), then the function must be applied on the columns of X
+% individually, otherwise fun and funt are assumed to be "vectorized".
+
+function retval = MatFun(sz,fun,funt,isvec)
+
+if nargin == 3
+ isvec = 0;
+end
+
+self.sz = sz;
+self.fun = fun;
+self.funt = funt;
+self.isvec = isvec;
+
+retval = class(self,'MatFun');
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/ctranspose.m b/inst/@MatFun/ctranspose.m
new file mode 100644
index 0000000..1506d4a
--- /dev/null
+++ b/inst/@MatFun/ctranspose.m
@@ -0,0 +1,25 @@
+% Conjugate transpose of a MatFun matrix.
+%
+% T = ctranspose(C)
+%
+% Returns the conjugate transpose of matrix C.
+
+function T = ctranspose(self)
+
+T = MatFun(self.sz([2 1]),self.funt,self.fun);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/double.m b/inst/@MatFun/double.m
new file mode 100644
index 0000000..dae75e7
--- /dev/null
+++ b/inst/@MatFun/double.m
@@ -0,0 +1,24 @@
+% Convert a MatFun matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the MatFun C matrix to a regular matrix F.
+
+function d = double(self)
+
+d = full(self);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/end.m b/inst/@MatFun/end.m
new file mode 100644
index 0000000..6a43e58
--- /dev/null
+++ b/inst/@MatFun/end.m
@@ -0,0 +1,23 @@
+% Last index of a MatFun matrix.
+%
+% e = end(C,k)
+%
+% Return the last index of the MatFun matrix C along dimension k.
+
+function e = end(self,k,n)
+
+e = self.sz(k);
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/full.m b/inst/@MatFun/full.m
new file mode 100644
index 0000000..99357c8
--- /dev/null
+++ b/inst/@MatFun/full.m
@@ -0,0 +1,24 @@
+% Convert a MatFun matrix to a regular matrix.
+%
+% F = full(C)
+%
+% Convert the MatFun C matrix to a regular matrix F.
+
+function M = full(self);
+
+M = self * eye(self.sz(2));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/mtimes.m b/inst/@MatFun/mtimes.m
new file mode 100644
index 0000000..8333c10
--- /dev/null
+++ b/inst/@MatFun/mtimes.m
@@ -0,0 +1,40 @@
+% Matrix product A*B.
+%
+% p = mtimes(A,B)
+%
+% Return the matrix product between the MatFun matrix A and another matrix B.
+
+function p = mtimes(A,B)
+
+if size(A,2) ~= size(B,1)
+ error('Inner matrix dimensions must agree.');
+end
+
+if isa(A,'MatFun')
+ if A.isvec
+ p = A.fun(B);
+ else
+ p = zeros(size(A,1),size(B,2));
+
+ for l=1:size(B,2)
+ p(:,l) = A.fun(B(:,l));
+ end
+ end
+else
+ error('not implemented');
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/@MatFun/size.m b/inst/@MatFun/size.m
new file mode 100644
index 0000000..c5a4750
--- /dev/null
+++ b/inst/@MatFun/size.m
@@ -0,0 +1,31 @@
+% Size of a MatFun matrix.
+%
+% sz = size(C)
+% sz = size(C,dim)
+%
+% Return the size of the MatFun matrix C. If dim is specified, then only
+% the size along dimension dim is returned.
+
+function sz = size(self,dim);
+
+sz = self.sz;
+
+if nargin == 2
+ sz = sz(dim);
+end
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/cat_cell_array.m b/inst/cat_cell_array.m
new file mode 100644
index 0000000..471ff79
--- /dev/null
+++ b/inst/cat_cell_array.m
@@ -0,0 +1,56 @@
+% Concatenate a cell array.
+%
+% x = cat_cell_array(c)
+%
+% Concatenate all arrays of the cell array c along the dimension length(c)+1
+
+
+function x = cat_cell_array(x)
+
+do_grid = 0;
+
+if iscell(x)
+ n = length(x);
+
+ if do_grid
+ all_vectors = 1;
+
+ for i=1:n
+ all_vectors = all_vectors & isvector(x{i});
+ end
+
+ if all_vectors
+ [x{:}] = ndgrid(x{:});
+ end
+ end
+
+ sz = size(x{1});
+ for i=2:n
+ if ~isequal(sz,size(x{i}))
+ %sz
+ %size(x{i})
+ error('all arguments in cell array should have the same size');
+ end
+ end
+
+ x = cat(n+1,x{:});
+ end
+end
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/colordisp.m b/inst/colordisp.m
new file mode 100644
index 0000000..b5a1e2d
--- /dev/null
+++ b/inst/colordisp.m
@@ -0,0 +1,27 @@
+% Display a message in a color (followed by a newline).
+%
+% colordisp(msg,color)
+%
+% The same as disp but in color.
+%
+% See also
+% colormsg
+
+function colordisp(msg,color)
+colormsg (msg,color)
+fprintf(1,'\n');
+
+% Copyright (C) 2004 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/colormsg.m b/inst/colormsg.m
new file mode 100644
index 0000000..44ed938
--- /dev/null
+++ b/inst/colormsg.m
@@ -0,0 +1,71 @@
+% Display a message in a color (without a newline).
+%
+% colormsg(msg,color)
+%
+% The message msg in printed on the terminal the specified color. Acceptable
+% values color color are:
+% black
+% red
+% green
+% yellow
+% blue
+% magenta
+% cyan
+% white
+%
+% In Matlab the message is displayed in black (due to the lack of ANSI escape
+% code support).
+%
+% See also
+% colordisp
+
+function colormsg (msg,color)
+
+%getenv('TERM')
+%if strcmp(getenv('TERM'),'xterm') % && exist('puts') > 1
+% only in octave
+if exist('puts') > 1
+ esc = char(27);
+
+ % ANSI escape codes
+ colors.black = [esc, '[40m'];
+ colors.red = [esc, '[41m'];
+ colors.green = [esc, '[42m'];
+ colors.yellow = [esc, '[43m'];
+ colors.blue = [esc, '[44m'];
+ colors.magenta = [esc, '[45m'];
+ colors.cyan = [esc, '[46m'];
+ colors.white = [esc, '[47m'];
+
+ reset = [esc, '[0m'];
+
+ c = getfield(colors,color);
+
+ %oldpso = page_screen_output (0);
+
+ try
+ fprintf([c, msg, reset]);
+ %puts (reset); % reset terminal
+ catch
+ %page_screen_output (oldpso);
+ end
+else
+ fprintf(msg);
+end
+
+
+
+% Copyright (C) 2004 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/conjugategradient.m b/inst/conjugategradient.m
new file mode 100644
index 0000000..5ec2c9d
--- /dev/null
+++ b/inst/conjugategradient.m
@@ -0,0 +1,198 @@
+% Solve a linear system with the preconditioned conjugated-gradient method.
+%
+% [x] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
+% [x,Q,T] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
+%
+% solves for x
+% A x = b
+% using the preconditioned conjugated-gradient method.
+% It provides also an approximation of A:
+% A \sim Q*T*Q'
+%
+% J(x) = 1/2 (x' b - x' A x)
+% ∇ J = b - A x
+% A x = b - ∇ J
+% b = ∇ J(0)
+%
+% the columns of Q are the Lanczos vectors
+
+% Alexander Barth 2010,2012-2013
+
+function [x,Q,T,diag] = conjugategradient(fun,b,varargin)
+
+n = length(b);
+
+% default parameters
+maxit = min(n,20);
+minit = 0;
+tol = 1e-6;
+pc = [];
+x0 = [];
+renorm = 0;
+
+prop = varargin;
+
+for i=1:2:length(prop)
+ if strcmp(prop{i},'minit')
+ minit = prop{i+1};
+ elseif strcmp(prop{i},'maxit')
+ maxit = prop{i+1};
+ elseif strcmp(prop{i},'tol')
+ tol = prop{i+1};
+ elseif strcmp(prop{i},'x0')
+ x0 = prop{i+1};
+ elseif strcmp(prop{i},'pc')
+ pc = prop{i+1};
+ elseif strcmp(prop{i},'renorm')
+ renorm = prop{i+1};
+ else
+
+ end
+end
+
+
+if isempty(x0)
+ % random initial vector
+ x0 = randn(n,1);
+end
+
+
+if isempty(pc)
+ % no pre-conditioning
+ pc = @(x) x;
+elseif isnumeric(pc)
+ invM = pc;
+ pc = @(x)( invM *x);
+end
+
+
+tol2 = tol^2;
+
+delta = [];
+gamma = [];
+q = NaN*zeros(n,1);
+
+%M = inv(invM);
+%E = chol(M);
+
+
+% initial guess
+x = x0;
+
+% gradient at initial guess
+r = b - fun(x);
+
+% quick exit
+if r'*r < tol2
+ return
+end
+
+% apply preconditioner
+z = pc(r);
+
+% first search direction == gradient
+p = z;
+
+% compute: r' * inv(M) * z (we will need this product at several
+% occasions)
+
+zr_old = r'*z;
+
+% r_old: residual at previous iteration
+r_old = r;
+
+for k=1:maxit
+ if k <= n && nargout > 1
+ % keep at most n vectors
+ Q(:,k) = r/sqrt(zr_old);
+ end
+
+ % compute A*p
+ Ap = fun(p);
+ %maxdiff(A*p,Ap)
+
+ % how far do we need to go in direction p?
+ % alpha is determined by linesearch
+
+ % alpha z'*r / (p' * A * p)
+ alpha(k) = zr_old / ( p' * Ap);
+
+ % get new estimate of x
+ x = x + alpha(k)*p;
+
+ % recompute gradient at new x. Could be done by
+ % r = b-fun(x);
+ % but this does require an new call to fun
+ r = r - alpha(k)*Ap;
+
+ if renorm
+ r = r - Q(:,1:k) * Q(:,1:k)' * r ;
+ end
+ %maxdiff(r,b-fun(x))
+
+ % apply pre-conditionner
+ z = pc(r);
+
+
+ zr_new = r'*z;
+
+ if r'*r < tol2 && k >= minit
+ break
+ end
+
+ %Fletcher-Reeves
+ beta(k+1) = zr_new / zr_old;
+ %Polak-Ribiere
+ %beta(k+1) = r'*(r-r_old) / zr_old;
+ %Hestenes-Stiefel
+ %beta(k+1) = r'*(r-r_old) / (p'*(r-r_old));
+ %beta(k+1) = r'*(r-r_old) / (r_old'*r_old);
+
+
+ % norm(p)
+ p = z + beta(k+1)*p;
+ zr_old = zr_new;
+ r_old = r;
+end
+
+
+%disp('alpha and beta')
+%figure,plot(beta(2:end))
+%rg(alpha)
+%rg(beta(2:end))
+
+if nargout > 1
+ kmax = size(Q,2);
+
+ delta(1) = 1/alpha(1);
+
+ %delta(1) - Q(:,1)'*invM*A*invM*Q(:,1)
+
+
+ for k=1:kmax-1
+ delta(k+1) = 1/alpha(k+1) + beta(k+1)/alpha(k);
+ gamma(k) = -sqrt(beta(k+1))/alpha(k);
+ end
+
+ T = sparse([1:kmax 1:kmax-1 2:kmax ],...
+ [1:kmax 2:kmax 1:kmax-1],...
+ [delta gamma gamma]);
+
+ diag.iter = k;
+ diag.relres = sqrt(r'*r);
+end
+
+% Copyright (C) 2004 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand.m b/inst/divand.m
new file mode 100644
index 0000000..1b52df2
--- /dev/null
+++ b/inst/divand.m
@@ -0,0 +1,290 @@
+% Compute a variational analysis of arbitrarily located observations.
+%
+% [fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
+% [fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...);
+% [fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
+%
+% Perform an n-dimensional variational analysis of the observations f located at
+% the coordinates x. The array fi represent the interpolated field at the grid
+% defined by the coordinates xi and the scales factors pmn.
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+%
+% pmn: scale factor of the grid. pmn is a cell array with n elements. Every
+% element represents the scale factor of the corresponding dimension. Its
+% inverse is the local resolution of the grid in a particular dimension.
+%
+% xi: cell array with n elements. Every element represents a coordinate
+% of the final grid on which the observations are interpolated
+% x: cell array with n elements. Every element represents a coordinate of
+% the observations
+% f: value of the observations *minus* the background estimate (m-by-1 array).
+% (see note)
+% len: correlation length
+% lambda: signal-to-noise ratio of observations (if lambda is a scalar).
+% The larger this value is, the closer is the field fi to the
+% observation.
+% if lambda is a scalar:
+% R = 1/lambda I, where R is the observation error covariance
+% matrix),
+% if lambda is a vector
+% a vector (R = diag(lambda)) or a matrix or a CovarParam object (R =
+% lambda).
+%
+%
+% Optional input arguments specified as pairs of keyword and values:
+% 'velocity', vel: velocity of advection constraint. The default is
+% no-advection constraint
+%
+% 'alpha': alpha is vector of coefficients multiplying various terms in the
+% cost function. The first element multiplies the norm.
+% The other i-th element of alpha multiplies the (i+1)-th derivative.
+% Per default, the highest derivative is m = ceil(1+n/2) where n is the
+% dimension of the problem.
+%
+% The values of alpha is the (m+1)th row of the Pascal triangle:
+% m=0 1
+% m=1 1 1
+% m=1 1 2 1 (n=1,2)
+% m=2 1 3 3 1 (n=3,4)
+% ...
+%
+% 'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or
+% off (0, default). If on, they will be returned as the last output
+% argument
+%
+% 'EOF', EOF: sub-space constraint. Orthogonal (EOF' WE^2 EOF = I) (units of
+% EOF: m^(-n/2))
+%
+% 'EOF_scaling', EOF_scaling: (dimensional)
+%
+% 'constraint': a structure with user specified constrain
+%
+% 'moddim': modulo for cyclic dimension (vector with n elements).
+% Zero is used for non-cyclic dimensions. Halo points should
+% not be included for cyclic dimensions. For example if the first dimension
+% is cyclic, then the grid point corresponding to mask(1,j) should be
+% between mask(end,1) (left neighbor) and mask(2,j) (right neighbor)
+%
+% 'fracdim': fractional indices (n-by-m array). If this array is specified,
+% then x and xi are not used.
+%
+% 'inversion': direct solver ('chol' for Cholesky factorization) or a
+% interative solver ('pcg' for preconditioned conjugate gradient) can be
+% used.
+%
+% 'compPC': function that returns a preconditioner for the primal formulation
+% if inversion is set to 'pcg'. The function has the following arguments:
+%
+% [M1,M2] = compPC(iB,H,R)
+%
+% where iB is the inverse background error covariance, H the observation
+% operator and R the error covariance of the observation. The used
+% preconditioner M will be M = M1 * M2 (or just M = M1 if M2 is empty).
+% Per default a modified incomplete Cholesky factorization will be used a
+% preconditioner.
+%
+% Note: 'velocity' and 'constraint' may appear multiple times
+%
+% Output:
+% fi: the analysed field
+% err: error variance of the analysis field relative to the error variance of
+% the background
+% s: structure
+% s.iB: adimensional
+% s.E: scaled EOF (dimensional)
+%
+% Note:
+% If zero is not a valid first guess for your variable (as it is the case for
+% e.g. ocean temperature), you have to subtract the first guess from the
+% observations before calling divand and then add the first guess back in.
+%
+% Example:
+% see divand_simple_example.m
+%
+
+
+function varargout = divand(mask,pmn,xi,x,f,len,lambda,varargin)
+
+% default values
+
+velocity = {};
+EOF = [];
+diagnostics = 0;
+EOF_lambda = 0;
+primal = 1;
+factorize = 1;
+tol = 1e-6;
+maxit = 100;
+minit = 10;
+constraints = {};
+inversion = 'chol';
+moddim = [];
+fracindex = [];
+alpha = [];
+keepLanczosVectors = 0;
+compPC = @divand_pc_sqrtiB;
+
+prop = varargin;
+for i=1:length(prop)/2
+ if strcmp(prop{2*i-1},'velocity')
+ velocity{end+1} = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'EOF')
+ EOF = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'EOF_scaling')
+ EOF_lambda = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'primal')
+ primal = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'factorize')
+ factorize = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'inversion')
+ inversion = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'tol')
+ tol = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'maxit')
+ maxit = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'minit')
+ minit = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'diagnostics')
+ diagnostics = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'constraint')
+ constraints{end+1} = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'fracindex')
+ fracindex = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'moddim')
+ moddim = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'alpha')
+ alpha = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'keepLanczosVectors')
+ keepLanczosVectors = prop{2*i};
+ elseif strcmpi(prop{2*i-1},'compPC')
+ compPC = prop{2*i};
+ else
+ warning('divand:unknownProperty','unknown property "%s" is irgnored',prop{2*i-1});
+ end
+end
+
+% check inputs
+
+if ~any(mask(:))
+ error('no sea points in mask');
+end
+
+s = divand_background(mask,pmn,len,alpha,moddim);
+s.betap = 0;
+s.EOF_lambda = EOF_lambda;
+s.primal = primal;
+s.factorize = factorize;
+s.tol = tol;
+s.maxit = maxit;
+s.minit = minit;
+s.inversion = inversion;
+s.keepLanczosVectors = keepLanczosVectors;
+s.compPC = compPC;
+
+% remove non-finite elements from observations
+f = f(:);
+valid = isfinite(f);
+x = cat_cell_array(x);
+
+if ~all(valid)
+ fprintf(1,'remove %d (out of %d) non-finite elements from observation vector\n',sum(~valid),numel(f));
+ x = reshape(x,[length(f) s.n]);
+ f = f(valid);
+ x = reshape(x(repmat(valid,[1 s.n])),[length(f) s.n]);
+
+ if ~isempty(fracindex)
+ fracindex = fracindex(:,valid);
+ end
+
+ if isscalar(lambda)
+ % do nothing
+ elseif isvector(lambda)
+ lambda = lambda(valid);
+ elseif ismatrix(lambda)
+ lambda = lambda(valid,valid);
+ end
+end
+
+apply_EOF_contraint = ~(isempty(EOF) | all(EOF_lambda == 0));
+
+s.mode = 1;
+
+if ~apply_EOF_contraint
+ s.betap = 0;
+else
+ if s.mode==0
+ s.betap = max(EOF_lambda)/s.coeff; % units m^(-n)
+ elseif s.mode==1
+ s.betap = max(max(EOF_lambda)-1,0)/s.coeff;
+ end
+end
+
+%assert(s.betap,0,1e-8)
+% increase contraint on total enegery to ensure system is still positive defined
+%s.betap
+s.iB = s.iB + s.betap * s.WE'*s.WE;
+
+% add observation constrain to cost function
+s = divand_addc(s,divand_obs(s,xi,x,f,lambda,fracindex));
+
+% add advection constraint to cost function
+for i=1:length(velocity)
+ %s = divand_advection(s,velocity);
+ s = divand_addc(s,divand_constr_advec(s,velocity{i}));
+end
+
+
+% add all additional constrains
+for i=1:length(constraints)
+ s = divand_addc(s,constraints{i});
+end
+
+
+if apply_EOF_contraint
+ s = divand_eof_contraint(s,EOF_lambda,EOF);
+end
+
+% factorize a posteori error covariance matrix
+% or compute preconditioner
+s = divand_factorize(s);
+
+if ~apply_EOF_contraint
+ [fi,s] = divand_solve(s,f);
+else
+ [fi,s] = divand_solve_eof(s,f);
+end
+
+varargout{1} = fi;
+
+if nargout-diagnostics >= 2
+ err = divand_error(s);
+ varargout{2} = err;
+end
+
+if diagnostics
+ s.B = CovarIS(s.iB);
+ [s.Jb,s.Jo,s.Jeof,s.J] = divand_diagnose(s,fi,f);
+ s.valid = valid;
+ varargout{nargout} = s;
+end
+
+
+% Copyright (C) 2008-2014 Alexander Barth <barth.alexander at gmail.com>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
+% LocalWords: fi divand pmn len diag CovarParam vel ceil moddim fracdim
\ No newline at end of file
diff --git a/inst/divand_addc.m b/inst/divand_addc.m
new file mode 100644
index 0000000..68cdb9a
--- /dev/null
+++ b/inst/divand_addc.m
@@ -0,0 +1,41 @@
+% Add a constraint to the cost function.
+%
+% s = divand_addc(s,constrain)
+%
+% Include in the structure s the specified constrain.
+%
+% Input:
+% s: structure created by divand_background
+% constrain: The parameter constrain has the following fields: R (a covariance
+% matrix), H (extraction operator) and yo (specified value for the
+% constrain).
+%
+% Output:
+% s: structure to be used by divand_factorize
+
+function s = divand_addc(s,constrain)
+
+if isfield(s,'R')
+ s.R = append(s.R,constrain.R);
+ s.H = append(s.H,constrain.H);
+ s.yo = cat(1,s.yo,constrain.yo);
+else
+ s.H = CatBlockMat(1,constrain.H);
+ s.R = DiagBlockMat(constrain.R);
+ s.yo = constrain.yo;
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_background.m b/inst/divand_background.m
new file mode 100644
index 0000000..1afdc1e
--- /dev/null
+++ b/inst/divand_background.m
@@ -0,0 +1,275 @@
+% Form the inverse of the background error covariance matrix.
+%
+% s = divand_background(mask,pmn,Labs,alpha,moddim)
+%
+% Form the inverse of the background error covariance matrix with
+% finite-difference operators on a curvilinear grid
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+%
+% pmn: scale factor of the grid.
+%
+% Labs: correlation length
+%
+% alpha: a dimensional coefficients for norm, gradient, laplacian,...
+% alpha is usually [1 2 1].
+%
+%
+% Output:
+% s: stucture containing
+% s.iB: inverse of the background error covariance
+% s.L: spatial average correlation length
+% s.n: number of dimenions
+% s.coeff: scaling coefficient such that the background variance
+% diag(inv(iB)) is one far away from the boundary.
+
+function s = divand_background(mask,pmn,Labs,alpha,moddim,mapindex)
+
+sz = size(mask);
+% number of dimensions
+pmn = cat_cell_array(pmn);
+n = size(pmn,ndims(pmn));
+
+if isempty(moddim)
+ moddim = zeros(1,n);
+end
+
+if nargin == 5
+ mapindex = [];
+end
+
+iscyclic = moddim > 0;
+
+
+if isscalar(Labs)
+ Labs = Labs * ones([size(mask) n]);
+elseif iscell(Labs)
+ for i=1:n
+ if isscalar(Labs{i})
+ Labs{i} = Labs{i} * ones(size(mask));
+ else
+ if ~isequal(size(mask),size(Labs{i}))
+ error('mask (%s) and correlation length (%s) have incompatible size',...
+ formatsize(size(mask)),formatsize(size(Labs{i})));
+ end
+ end
+ end
+
+ Labs = cat_cell_array(Labs);
+else
+ if ~isequal(size(mask),size(Labs))
+ error('mask (%s) and correlation length (%s) have incompatible size',...
+ formatsize(size(mask)),formatsize(size(Labs)));
+ end
+
+ Labs = repmat(Labs,[ones(1,ndims(mask)) n]);
+end
+
+if isempty(alpha)
+ % kernel should has be continuous derivative
+
+ % highest derivative in cost function
+ m = ceil(1+n/2);
+
+ % alpha is the (m+1)th row of the Pascal triangle:
+ % m=0 1
+ % m=1 1 1
+ % m=1 1 2 1
+ % m=2 1 3 3 1
+ % ...
+
+ alpha = arrayfun (@(k) nchoosek(m,k), 0:m);
+end
+
+%if ~isequal([size(mask) n],size(pmn))
+% error('mask (%s) and metric (%s) have incompatible size',formatsize(size(mask)),formatsize(size(pmn)));
+%end
+
+s = divand_operators(mask,pmn,Labs.^2,iscyclic,mapindex);
+D = s.D; % laplacian (a dimensional, since nu = Labs.^2)
+sv = s.sv;
+n = s.n;
+
+
+% Labsp: 1st index represents the dimensions
+Labsp = permute(Labs,[n+1 1:n]);
+pmnp = permute(pmn,[n+1 1:n]);
+
+% must handle the case when Labs is zero in some dimension
+% thus reducing the effective dimension
+
+% mean correlation length in every dimension
+Ld = mean(reshape(Labsp,[n sv.numels_all]),2);
+neff = sum(Ld > 0);
+%L = mean(Ld(Ld > 0));
+
+% geometric mean
+L = geomean(Ld(Ld > 0));
+
+
+% norm taking only dimension into account with non-zero correlation
+% WE: units length^(neff/2)
+d = prod(pmnp(find(Ld > 0),:),1)';
+WE = sparse_diag(statevector_pack(sv,1./sqrt(d)));
+
+
+% scale iB such that the diagonal of inv(iB) is 1 far from
+% the boundary
+% we use the effective dimension neff to take into account that the
+% correlation length-scale might be zero in some directions
+
+Ln = prod(Ld(Ld > 0));
+
+if any(Ld <= 0)
+ pmnd = mean(reshape(pmnp,[n sv.numels_all]),2);
+ %Ln = Ln * prod(pmnd(Ld <= 0));
+end
+
+
+try
+ coeff = divand_kernel(neff,alpha) * Ln; % units length^n
+catch ME
+ if strcmp(ME.identifier, 'divand:divand_kernel:unsupported_alpha')
+ warning('divand:divand_background:noscaling','no scaling');
+ coeff = 1;
+ else
+ rethrow(ME);
+ end
+end
+
+pmnv = reshape(pmn,numel(mask),n);
+pmnv(:,Ld == 0) = 1;
+
+% staggered version of norm
+
+for i=1:n
+ S = sparse_stagger(sz,i,iscyclic(i));
+ ma = (S * mask(:) == 1);
+ d = sparse_pack(ma) * prod(S * pmnv,2);
+ d = 1./d;
+ s.WEs{i} = sparse_diag(sqrt(d));
+end
+
+% staggered version of norm scaled by length-scale
+
+s.WEss = cell(n,1);
+s.Dxs = cell(n,1);
+for i=1:n
+ Li2 = Labsp(i,:).^2;
+
+ S = sparse_stagger(sz,i,iscyclic(i));
+ % mask for staggered variable
+ m = (S * mask(:) == 1);
+
+ tmp = sparse_pack(m) * sqrt(S*Li2(:));
+
+ s.WEss{i} = sparse_diag(tmp) * s.WEs{i};
+% s.Dxs{i} = sparse_diag(sqrt(tmp)) * s.Dx{i};
+end
+
+% adjust weight of halo points
+if ~isempty(mapindex)
+ % ignore halo points at the center of the cell
+
+ WE = sparse_diag(s.isinterior) * WE;
+
+ % divide weight be two at the edged of halo-interior cell
+ % weight of the grid points between halo and interior points
+ % are 1/2 (as there are two) and interior points are 1
+ for i=1:n
+ s.WEs{i} = sparse_diag(sqrt(s.isinterior_stag{i})) * s.WEs{i};
+ end
+end
+
+s.WE = WE;
+s.coeff = coeff;
+% number of dimensions
+s.n = n;
+
+[iB_,iB] = divand_background_components(s,alpha);
+
+if 0
+iB_ = cell(length(alpha),1);
+
+% constrain of total norm
+
+iB_{1} = (1/coeff) * (WE'*WE);
+
+
+% loop over all derivatives
+
+for j=2:length(alpha)
+ % exponent of laplacian
+ k = floor((j-2)/2);
+
+ if mod(j,2) == 0
+ % constrain of derivative with uneven order (j-1)
+ % (gradient, gradient*laplacian,...)
+ % normalized by surface
+
+ iB_{j} = sparse(size(D,1),size(D,1));
+
+ for i=1:n
+ Dx = s.WEss{i} * s.Dx{i} * D^k;
+
+ iB_{j} = iB_{j} + Dx'*Dx;
+ end
+ else
+ % constrain of derivative with even order (j-1)
+ % (laplacian, biharmonic,...)
+
+ % normalize by surface of each cell
+ % such that inner produces (i.e. WE'*WE)
+ % become integrals
+ % WD: units length^(n/2)
+
+ WD = WE * D^(k+1);
+ iB_{j} = WD'*WD;
+ end
+
+ iB_{j} = iB_{j}/coeff;
+end
+
+
+% sum all terms of iB
+% iB is adimentional
+iB = alpha(1) * iB_{1};
+for j=2:length(alpha)
+ iB = iB + alpha(j) * iB_{j};
+end
+
+s.coeff = coeff;
+end
+
+% inverse of background covariance matrix
+s.iB = iB;
+
+% individual terms of background covariance matrix (corresponding alpha)
+s.iB_ = iB_;
+
+s.L = L;
+s.Ln = Ln;
+
+s.moddim = moddim;
+s.iscyclic = iscyclic;
+
+s.alpha = alpha;
+s.neff = neff;
+s.WE = WE; % units length^(n/2)
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_background_components.m b/inst/divand_background_components.m
new file mode 100644
index 0000000..e02a7c4
--- /dev/null
+++ b/inst/divand_background_components.m
@@ -0,0 +1,80 @@
+% Form the different components of the background error covariance matrix.
+%
+% [iB_,iB] = divand_background_components(s,alpha)
+%
+% Compute the components of the background error covariance matrix iB_ and
+% their sum based on alpha (the a-dimensional coefficients for norm, gradient,
+% laplacian,...).
+
+function [iB_,iB] = divand_background_components(s,alpha)
+
+WE = s.WE;
+D = s.D;
+coeff = s.coeff;
+n = s.n;
+
+iB_ = cell(length(alpha),1);
+
+% constrain of total norm
+
+iB_{1} = (1/coeff) * (WE'*WE);
+
+
+% loop over all derivatives
+
+for j=2:length(alpha)
+ % exponent of laplacian
+ k = floor((j-2)/2);
+
+ if mod(j,2) == 0
+ % constrain of derivative with uneven order (j-1)
+ % (gradient, gradient*laplacian,...)
+ % normalized by surface
+
+ iB_{j} = sparse(size(D,1),size(D,1));
+
+ for i=1:n
+ Dx = s.WEss{i} * s.Dx{i} * D^k;
+ %Dx = s.WEs{i} * s.Dxs{i} * D^k;
+ iB_{j} = iB_{j} + Dx'*Dx;
+ end
+ else
+ % constrain of derivative with even order (j-1)
+ % (laplacian, biharmonic,...)
+
+ % normalize by surface of each cell
+ % such that inner produces (i.e. WE'*WE)
+ % become integrals
+ % WD: units length^(n/2)
+
+ WD = WE * D^(k+1);
+ iB_{j} = WD'*WD;
+ end
+
+ iB_{j} = iB_{j}/coeff;
+end
+
+
+% sum all terms of iB
+% iB is adimentional
+iB = alpha(1) * iB_{1};
+for j=2:length(alpha)
+ iB = iB + alpha(j) * iB_{j};
+end
+
+% LocalWords: iB divand
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_constr_advec.m b/inst/divand_constr_advec.m
new file mode 100644
index 0000000..b3c3757
--- /dev/null
+++ b/inst/divand_constr_advec.m
@@ -0,0 +1,67 @@
+% Create the advection constrain.
+%
+% c = divand_constr_advec(s,velocity)
+%
+% Create the advection constrain using the specified velocity.
+%
+% Input:
+% s: structure created by divand_background
+% velocity: cell array of velocity vectors
+%
+% Output:
+% c: structure to be used by divand_addc with the following fields: R (a
+% covariance matrix), H (extraction operator) and yo (specified value for
+% the constrain).
+
+function c = divand_constr_advec(s,velocity)
+
+velocity = cat_cell_array(velocity);
+
+% check for NaNs
+nancount = sum(isnan(velocity(:)));
+if nancount > 0
+ error('%d velocity values are equal to NaN',nancount);
+end
+
+mask = s.mask;
+n = s.n;
+iscyclic = s.iscyclic;
+
+sz = size(mask);
+
+A = sparse(0);
+
+vel = reshape(velocity,numel(mask),n);
+
+
+
+for i=1:n
+ S = sparse_stagger(sz,i,iscyclic(i));
+ m = (S * mask(:) == 1);
+
+ d = vel(:,i);
+
+ A = A + sparse_diag(d(mask==1)) * sparse_pack(mask) * S' * sparse_pack(m)' * s.Dx{i};
+end
+
+l = size(A,1);
+
+c.H = A;
+c.yo = sparse(l,1);
+c.R = speye(l);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_diagnose.m b/inst/divand_diagnose.m
new file mode 100644
index 0000000..2f7a515
--- /dev/null
+++ b/inst/divand_diagnose.m
@@ -0,0 +1,47 @@
+% Computes diagnostics for the analysis.
+%
+% [Jb,Jo,Jeof,J] = divand_diagnose(s,fi,d)
+%
+% Computes the value of different terms of the cost-function
+%
+% Input:
+% s: structure created by divand_solve
+% fi: analysis
+% d: observations
+%
+% Output:
+% Jb: background constrain
+% Jeof: observation constrain
+% Jeof: EOF constrain
+
+function [Jb,Jo,Jeof,J] = divand_diagnose(s,fi,d)
+
+d = s.yo;
+xa = statevector_pack(s.sv,fi);
+
+Jb = xa' * s.iB * xa - s.betap * (s.WE'*xa)' * (s.WE'*xa);
+Jo = (s.H*xa - d)' * (s.R \ (s.H*xa - d));
+
+if all(s.EOF_lambda == 0)
+ Jeof = 0;
+else
+ Jeof = s.betap * ((s.WE'*xa)' * (s.WE'*xa)) - (s.E'*xa)' * (s.E'*xa);
+end
+
+Jo = full(Jo);
+
+J = Jb + Jo + Jeof;
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_eof_contraint.m b/inst/divand_eof_contraint.m
new file mode 100644
index 0000000..3964317
--- /dev/null
+++ b/inst/divand_eof_contraint.m
@@ -0,0 +1,68 @@
+% Include constraint from EOFs.
+%
+% s = divand_eof_contraint(s,EOF_lambda,EOF)
+%
+% Include the constraint from the EOF and their eigenvalues (EOF_lambda)
+% in the cost function described by the structure s.
+%
+
+function s = divand_eof_contraint(s,EOF_lambda,EOF)
+
+iB = s.iB;
+sv = s.sv;
+WE = s.WE;
+
+coeff = s.coeff;
+
+% remove land points
+E = statevector_pack(sv,EOF);
+% units m^(-n/2)
+
+% normalize by surface
+% for coeff see divand_background.m
+%E = WE * E;
+%EOF_lambda
+
+E = WE^2 * E * diag(sqrt(EOF_lambda / coeff));
+%% units: m^(n) m^(-n/2) sqrt(m^(-n)) = 1
+
+if 0
+
+BE = iB \ E;
+
+
+
+A = BE * inv(E'*BE - eye(size(E,2)));
+Beof_var = - sum(A .* BE,2);
+
+% scaling factor to enshure that:
+% scaling * inv(iB - E E') has a variance of 1 far from boundaries
+
+scaling = 1/(1 + mean(Beof_var))
+scaling = 1
+
+%disp('apply scaling')
+E = 1/sqrt(scaling) * E;
+iB = iB/scaling;
+
+s.scaling = scaling;
+
+end
+
+s.E = E;
+s.iB = iB;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_error.m b/inst/divand_error.m
new file mode 100644
index 0000000..854ab4d
--- /dev/null
+++ b/inst/divand_error.m
@@ -0,0 +1,88 @@
+% Compute the expected error variance of the analysis.
+%
+% err = divand_error(s)
+%
+% Compute the expected error variance of the analysis based on the
+% background error covariance, observation error covariance and data
+% distribution.
+%
+% Input:
+% s: structure created by divand_factorize
+%
+% Output:
+% err: expected error variance of the analysis
+
+
+function err = divand_error(s)
+
+H = s.H;
+sv = s.sv;
+N = sum(s.mask(:)==1);
+n = s.n;
+
+errp = zeros(N,1);
+
+if s.primal
+ errp = diag(s.P);
+
+ %global P_CL
+ %errp = diag(P_CL);
+else
+ %dual
+ B = s.B;
+ R = s.R;
+
+ % fun(x) computes (H B H' + R)*x
+ fun = @(x) H * (B * (H'*x)) + R*x;
+ t0 = cputime;
+
+ for i=1:N
+ % show progress every 5 seconds
+ if (t0 + 5 < cputime)
+ t0 = cputime;
+ fprintf('%d/%d\n',i,N);
+ end
+
+ e = zeros(N,1);
+ e(i) = 1;
+ % analyzed covariance
+ % P = B - B * H' * inv(H * B * H' + R) * H * B
+ % analyzed variance
+ % P = e' B e - e' * B * H' * inv(H * B * H' + R) * H * B * e
+
+ Be = B * e;
+ errp(i) = e' * Be;
+ HBe = H*Be;
+ %keyboard
+
+ % tmp is inv(H * B * H' + R) * H * B * e
+ %tic
+ [tmp,flag,relres,iter] = pcg(fun, HBe,s.tol,s.maxit,s.funPC);
+ %toc
+
+ if (flag ~= 0)
+ error('divand:pcg', ['Preconditioned conjugate gradients method'...
+ ' did not converge %d %g %g'],flag,relres,iter);
+ end
+
+ errp(i) = errp(i) - HBe' * tmp;
+ end
+end
+
+err = statevector_unpack(sv,errp);
+err(~s.mask) = NaN;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_factorize.m b/inst/divand_factorize.m
new file mode 100644
index 0000000..93cad16
--- /dev/null
+++ b/inst/divand_factorize.m
@@ -0,0 +1,82 @@
+% Factorize some matrices to increase performance.
+%
+% s = divand_factorize(s)
+%
+% Create the inverse of the posteriori covariance matrix
+% and factorize
+%
+% Input:
+% s: structure created by divand_background
+%
+% Output:
+% s.P: factorized P
+
+
+function s = divand_factorize(s)
+
+R = s.R;
+iB = s.iB;
+H = s.H;
+
+if s.primal
+ if strcmp(s.inversion,'chol')
+ iP = iB + H'*(R\H);
+ P = CovarIS(iP);
+
+ % Cholesky factor of the inverse of a posteriori
+ % error covariance iP
+ if s.factorize
+ P = factorize(P);
+ end
+
+ s.P = P;
+ else
+ tic
+ [s.M1,s.M2] = s.compPC(iB,H,R);
+ s.pc_time = toc();
+ end
+else % dual
+ %C = H * (iB \ H') + R;
+ s.B = CovarIS(iB);
+ % pre-conditioning function for conjugate gradient
+ s.funPC = [];
+
+ if s.factorize
+ s.B = factorize(s.B);
+
+ % iM = inv(H * B * H' + sparse_diag(diag(R)))
+ % iM * (H * B * H' + R) is identify matrix if R is diagonal
+
+ M = H * (s.B * H') + sparse_diag(diag(R));
+ %M = H * (s.B * H') + sparse_diag(sum(R,1));
+ iM = CovarIS(M);
+ iM = factorize(iM);
+
+ % pre-conditioning function for conjugate gradient
+ s.funPC = @(x) iM*x;
+
+ if 0
+ C = H * (s.B * H') + sparse_diag(diag(R));
+ [RC, q, QC] = chol (C);
+ s.funPC = @(x) RC' \ (QC*x);
+
+ end
+ end
+end
+
+% LocalWords: divand
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_kernel.m b/inst/divand_kernel.m
new file mode 100644
index 0000000..f0ad197
--- /dev/null
+++ b/inst/divand_kernel.m
@@ -0,0 +1,84 @@
+% Return the analytical kernel and normalization factor.
+%
+% [mu,K] = divand_kernel(n,alpha)
+% [mu,K] = divand_kernel(n,alpha,r)
+%
+% Analytical (normalized) kernels for infinite domain in dimension n and for
+% coefficients alpha
+% Input
+% n: number of dimensions
+% alpha: coefficients
+% r (optional): distance from origin
+% Output:
+% K: kernel function evaluate at the values of r if present or a function handle
+% mu: normalization factor
+
+function [mu,K,rh] = divand_kernel(n,alpha,r)
+
+
+% remove trailling zeros
+ind = max(find(alpha ~= 0));
+alpha = alpha(1:ind);
+
+m = length(alpha)-1;
+K = [];
+if isequal(arrayfun(@(k) nchoosek(m,k),0:m),alpha)
+ % alpha are binomial coefficients
+
+ [mu,K] = divand_kernel_binom(n,m);
+else
+ error('divand:divand_kernel:unsupported_alpha',...
+ 'unsupported sequence of alpha')
+end
+
+if nargin == 3
+ % evaluate the kernel for the given values of r
+ K = K(r);
+end
+
+if nargout == 3
+ % determine at which distance r K(r) = 1/2
+ rh = abs(fzero(@(r) K(abs(r))-.5,1));
+end
+
+end
+
+function [mu,K] = divand_kernel_binom(n,m)
+
+
+% %if isequal(alpha,1)
+
+nu = m-n/2;
+mu = (4*pi)^(n/2) * gamma(m) / gamma(nu);
+K = @(x) divand_rbesselk(nu,x);
+
+if nu <= 0
+ warning('divand:nonorm','No normalization possible. Extend parameter alpha.');
+ mu = 1;
+end
+
+end
+
+
+function [K] = divand_rbesselk(nu,r)
+r = abs(r);
+K = 2/gamma(nu) * ((r/2).^nu .* besselk(nu,r));
+K(r == 0) = 1;
+end
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_laplacian.m b/inst/divand_laplacian.m
new file mode 100644
index 0000000..a7fe537
--- /dev/null
+++ b/inst/divand_laplacian.m
@@ -0,0 +1,119 @@
+% Create the laplacian operator.
+%
+% Lap = divand_laplacian(mask,pmn,nu,iscyclic)
+%
+% Form a Laplacian using finite differences
+% assumes that gradient is zero at "coastline"
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+% pmn: scale factor of the grid.
+% nu: diffusion coefficient of the Laplacian
+% field of the size mask or cell arrays of fields
+%
+% Output:
+% Lap: sparce matrix represeting a Laplaciant
+%
+%
+function Lap = divand_laplacian(mask,pmn,nu,iscyclic)
+
+pmn = cat_cell_array(pmn);
+
+% number of dimensions
+n = size(pmn,ndims(pmn));
+
+sz = size(mask);
+
+if nargin == 2
+ nu = ones(sz);
+end
+
+if isequal(sz,size(nu))
+ % assume diffusion is the same along all dimensions
+ nu = repmat(nu,[ones(1,n) n]);
+end
+
+% nu: 1st index represents the dimensions
+nu = permute(nu,[n+1 1:n]);
+
+% extraction operator of sea points
+H = sparse_pack(mask==1);
+
+
+sz = size(mask);
+
+%DD = sparse(0);
+DD = sparse(prod(sz),prod(sz));
+
+Pmn = reshape(pmn,[prod(sz) n]);
+for i=1:n
+ % operator for staggering in dimension i
+ S = sparse_stagger(sz,i,iscyclic(i));
+
+ % d = 1 for interfaces surounded by water, zero otherwise
+ d = S * mask(:) == 1;
+
+ % metric
+ for j = 1:n
+ tmp = S * Pmn(:,j);
+
+ if j==i
+ d = d .* tmp;
+ else
+ d = d ./ tmp;
+ end
+ end
+
+ % tmp: "diffusion coefficient" along dimension i
+ tmp = nu(i,:);
+
+ %whos tmp d S nu
+ %keyboard
+ d = d .* (S * tmp(:));
+
+ % Flux operators D
+ % zero at "coastline"
+
+ D = sparse_diag(d) * sparse_diff(sz,i,iscyclic(i));
+ szt = sz;
+
+ if ~iscyclic(i)
+ % extx: extension operator
+ szt(i) = szt(i)+1;
+ extx = sparse_trim(szt,i)';
+
+ D = extx * D;
+ else
+ % shift back one grid cell along dimension i
+ D = sparse_shift(sz,i,iscyclic(i))' * D;
+ end
+
+ % add laplacian along dimension i
+ DD = DD + sparse_diff(szt,i,iscyclic(i)) * D;
+end
+
+ivol = prod(pmn,n+1);
+
+% Laplacian on regualar grid DD
+DD = sparse_diag(ivol) * DD;
+
+% Laplacian on grid with on sea points Lap
+Lap = H * DD * H';
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_metric.m b/inst/divand_metric.m
new file mode 100644
index 0000000..006a716
--- /dev/null
+++ b/inst/divand_metric.m
@@ -0,0 +1,56 @@
+% Compute metric scale factors.
+%
+% [pm,pn] = divand_metric(lon,lat)
+%
+% Compute metric scale factors pm and pn based on
+% longitude lon and latitude lat. The variables pm and pn
+% represent the inverse of the local resolution in meters using
+% the mean Earth radius.
+
+
+
+function [pm,pn] = divand_metric(lon,lat)
+
+sz = size(lon);
+i = 2:sz(1)-1;
+j = 2:sz(2)-1;
+
+
+dx = distance(lat(i-1,:),lon(i-1,:),lat(i+1,:),lon(i+1,:))/2;
+dx = cat(1,dx(1,:),dx,dx(end,:));
+
+dy = distance(lat(:,j-1),lon(:,j-1),lat(:,j+1),lon(:,j+1))/2;
+dy = cat(2,dy(:,1),dy,dy(:,end));
+
+dx = real(dx);
+dy = real(dy);
+
+dx = deg2m(dx);
+dy = deg2m(dy);
+
+
+pm = 1./dx;
+pn = 1./dy;
+
+end
+
+function dy = deg2m(dlat)
+% Mean radius (http://en.wikipedia.org/wiki/Earth_radius)
+R = 6371.009e3;
+
+dy = dlat*(2*pi*R)/360;
+end
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_obs.m b/inst/divand_obs.m
new file mode 100644
index 0000000..5d6ea9c
--- /dev/null
+++ b/inst/divand_obs.m
@@ -0,0 +1,83 @@
+% Include the constrain from the observations.
+%
+% s = divand_obs(s,xi,x,lambda,I)
+%
+% Set observations of variational problem.
+% It is assumed that the each coordinate depends only on one
+% index. If this is not the case, then matrix I must be provided.
+%
+% Input:
+% s: structure created by divand_background
+% xi: coordinates of observations*
+% x: coordinates of grid*
+% lambda: signal-to-noise ratio of observations
+% I (optional): fractional indexes of location of observation
+% within the grid
+%
+% Output:
+% s: structure to be used by divand_factorize
+%
+% Note:
+% *these parameters can either be specified as a cell
+% array of all dimenions:
+% xi = {Xi,Yi,Zi}
+% or as n+1 dimensional array
+
+function c = divand_obs(s,xi,x,yo,lambda,I)
+
+if nargin == 5
+ I = [];
+end
+
+xi = cat_cell_array(xi);
+x = cat_cell_array(x);
+
+mask = s.mask;
+iscyclic = s.iscyclic;
+moddim = s.moddim;
+
+
+if isempty(I)
+ I = localize_separable_grid(x,mask,xi);
+end
+
+[H,out] = sparse_interp(mask,I,iscyclic);
+
+nout = sum(out);
+if nout ~= 0
+ fprintf(1,'Observations out of domain: %d\n',nout);
+end
+
+
+H = H * sparse_pack(mask)';
+
+% iB is scaled such that diag(inv(iB)) is 1 far from the
+% boundary
+
+if isscalar(lambda)
+ R = 1/lambda * speye(size(H,1));
+elseif isvector(lambda)
+ R = sparse_diag(lambda);
+else
+ R = lambda;
+end
+
+c.H = H;
+s.out = out;
+c.R = R;
+c.yo = yo;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_operators.m b/inst/divand_operators.m
new file mode 100644
index 0000000..04d1148
--- /dev/null
+++ b/inst/divand_operators.m
@@ -0,0 +1,112 @@
+% Generate the gradient and Laplacian operators.
+%
+% s = divand_operators(mask,pmn,nu,iscyclic)
+%
+% Form sparse matrixes representing the gradient and Laplacian using
+% finite differences
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+% pmn: scale factor of the grid.
+% nu: diffusion coefficient of the Laplacian
+%
+% Output:
+% s: stucture containing
+% s.Dx: cell array of the gradient
+% s.D: Laplaciant
+% s.sv: structure describing the state vector
+% s.mask: land-sea mask
+% s.WE: diagonal matrix where each element is the surface
+% of a grid cell
+
+
+function s = divand_operators(mask,pmn,nu,iscyclic,mapindex)
+
+pmn = cat_cell_array(pmn);
+
+% number of dimensions
+n = size(pmn,ndims(pmn));
+
+sz = size(mask);
+
+sv = statevector_init(mask);
+
+if ~isempty(mapindex)
+ % mapindex is unpacked and referers to unpacked indices
+
+ % range of packed indices
+ % land point map to 1, but those points are remove by statevector_pack
+ i2 = statevector_unpack(sv,(1:sv.n)',1);
+ mapindex_packed = statevector_pack(sv,i2(mapindex));
+
+ % applybc*x applies the boundary conditions to x
+ i = [1:sv.n]';
+ applybc = sparse(i,mapindex_packed(i),ones(1,sv.n),sv.n,sv.n);
+
+ % a halo point is points which maps to a (different) interior point
+ % a interior point maps to itself
+ s.isinterior = (i == mapindex_packed(i));
+ s.isinterior_unpacked = statevector_unpack(sv,s.isinterior);
+
+ s.mapindex_packed = mapindex_packed;
+end
+
+D = divand_laplacian(mask,pmn,nu,iscyclic);
+
+% XXX remove this WE
+d = statevector_pack(sv,1./(prod(pmn,n+1)));
+WE = sparse_diag(sqrt(d));
+
+for i=1:n
+ S = sparse_stagger(sz,i,iscyclic(i));
+
+ % mask on staggered grid
+ ma = (S * mask(:) == 1);
+ s.mask_stag{i} = ma;
+
+ d = sparse_pack(ma) * prod(S * reshape(pmn,numel(mask),n),2);
+ d = 1./d;
+ s.WEs{i} = sparse_diag(sqrt(d));
+end
+
+s.Dx = cell(n,1);
+[s.Dx{:}] = sparse_gradient(mask,pmn,iscyclic);
+
+
+if ~isempty(mapindex)
+ D = applybc * D * applybc;
+ WE = sparse_diag(s.isinterior) * WE;
+
+ for i=1:n
+ S = sparse_stagger(sz,i,iscyclic(i));
+ s.isinterior_stag{i} = sparse_pack(s.mask_stag{i}) * S * s.isinterior_unpacked(:);
+
+ % the results of 's.Dx{i} * field' satisfies the bc is field does
+ % there is need to reapply the bc on the result
+ s.Dx{i} = s.Dx{i} * applybc;
+ end
+
+ s.applybc = applybc;
+end
+
+s.D = D;
+s.sv = sv;
+s.mask = mask;
+s.WE = WE;
+s.n = n;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_orthogonalize_modes.m b/inst/divand_orthogonalize_modes.m
new file mode 100644
index 0000000..cfecfdb
--- /dev/null
+++ b/inst/divand_orthogonalize_modes.m
@@ -0,0 +1,48 @@
+% Orthogonalize EOF modes.
+%
+% [U3,WE] = divand_orthogonalize_mode(mask,pmn,U)
+%
+% The modes U are orthogonalized using the mask and metric (pmn). WE is square
+% root of norm. Diagonal of WE.^2 is the surface of corresponding cell (units:
+% m). The output U3 represent the normalized EOFs (units: m-1)
+
+function [U3,WE] = divand_orthogonalize_modes(mask,pmn,U);
+
+pmn = cat_cell_array(pmn);
+% number of dimensions
+n = size(pmn,ndims(pmn));
+
+sv = statevector_init(mask);
+d = statevector_pack(sv,1./(prod(pmn,n+1)));
+
+% spare root of the norm
+WE = sparse_diag(sqrt(d));
+
+U2 = statevector_pack(sv,U);
+U2 = WE * U2;
+
+r = size(U2,2);
+[U3,S3,V3] = svds(U2,r);
+%U3 = reshape(U3,[size(mask) r]);
+
+U3 = inv(WE) * U3;
+
+%assert(U3' * WE^2 * U3, eye(r),1e-8)
+
+U3 = statevector_unpack(sv,U3);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_pc_michol.m b/inst/divand_pc_michol.m
new file mode 100644
index 0000000..acbf36b
--- /dev/null
+++ b/inst/divand_pc_michol.m
@@ -0,0 +1,50 @@
+% Compute a preconditioner using a modified incomplete Cholesky decomposition.
+%
+% [M1,M2] = divand_pc_michol(iB,H,R,icholparam)
+%
+% Compute preconditioner matrices M1 and M2 based on
+% the inverse background error covariance iB, observation operator
+% H and observation R covariance R. icholparam is a structure will parameters
+% for ichol. The default value is struct('michol','on'). A modified incomplete
+% Cholesky factorization of the matrix iP = iB + H'*(R\H) is computed per
+% default.
+% M2 is the transpose of M1 for this preconditioner.
+%
+% The function ichol is necessary.
+%
+% See also:
+% ichol, divand_pc_sqrtiB
+
+function [M1,M2] = divand_pc_michol(iB,H,R,icholparam)
+
+if nargin == 3
+ icholparam = struct('michol','on');
+end
+
+iP = iB + H'*(R\H);
+
+if which('ichol')
+ M1 = ichol(iP,icholparam);
+else
+ warning('divand-noichol','ichol is not available')
+ M1 = [];
+end
+
+M2 = M1';
+
+% LocalWords: preconditioner Cholesky diavnd pc michol iB icholparam ichol struct iP divand sqrtiB
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_pc_none.m b/inst/divand_pc_none.m
new file mode 100644
index 0000000..abb3739
--- /dev/null
+++ b/inst/divand_pc_none.m
@@ -0,0 +1,30 @@
+% No preconditioner is used.
+%
+% [M1,M2] = divand_pc_none(iB,H,R)
+%
+% Dummy function for requiring that no preconditioner is used in divand.
+%
+% See also:
+% diavnd_pc_michol, diavnd_pc_sqrtiB
+
+function [M1,M2] = divand_pc_none(iB,H,R)
+
+M1 = [];
+M2 = [];
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_pc_sqrtiB.m b/inst/divand_pc_sqrtiB.m
new file mode 100644
index 0000000..a5a0852
--- /dev/null
+++ b/inst/divand_pc_sqrtiB.m
@@ -0,0 +1,31 @@
+% Compute a preconditioner using the Cholesky decomposition.
+%
+% [M1,M2] = divand_pc_michol(iB,H,R)
+%
+% Compute preconditioner matrices M1 and M2 based on
+% the Cholesky decomposition of iB. The matrices H and R are not used.
+% M2 is the transpose of M1 for this preconditioner.
+%
+% See also:
+% chol, divand_pc_michol
+
+function [M1,M2] = divand_pc_sqrtiB(iB,H,R)
+
+M1 = chol(iB);
+M2 = M1';
+% LocalWords: preconditioner diavnd pc michol iB Cholesky chol divand sqrtiB
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_realistic_example.m b/inst/divand_realistic_example.m
new file mode 100644
index 0000000..aea762a
--- /dev/null
+++ b/inst/divand_realistic_example.m
@@ -0,0 +1,99 @@
+% A realistic example of divand in 2 dimensions
+% with salinity observations in the Mediterranean Sea at 30m.
+%
+% Data is available at
+% http://modb.oce.ulg.ac.be/mediawiki/index.php/Divand
+
+% load observations
+
+load data.dat
+
+% extract columns of data
+lonobs = data(:,1);
+latobs = data(:,2);
+S = data(:,3);
+
+% load bathymetry
+% bath is negative in water and positive in air
+bat = ncread('diva_bath.nc','bat');
+lon = ncread('diva_bath.nc','lon');
+lat = ncread('diva_bath.nc','lat');
+
+% compute grid metric
+% pm and pn are in meters^-1
+[lon,lat] = ndgrid(lon,lat);
+[pm,pn] = divand_metric(lon,lat);
+
+% compute mean and anomalies
+Smean = mean(S);
+Sanom = S - mean(S);
+
+% correlation length (in meters)
+len = 200e3;
+len = 100e3;
+% signal-to-noise ratio
+lambda = 50.5;
+
+% land-sea mask
+% mask everything below 30 m
+mask = bat < -30;
+
+% mask for the analysis
+maska = mask;
+% remove Atlantic
+maska(1:60,130:end) = 0;
+% remove Black Sea
+maska(323:end,100:end) = 0;
+
+% make analysis
+Sa = divand(maska,{pm,pn},{lon,lat},{lonobs,latobs},Sanom,len,lambda);
+
+% add mean back
+Sa2 = Sa + Smean;
+
+% create plots
+
+% color axis
+ca = [36.2520 39.4480];
+
+% aspect ratio
+ar = [1 cos(mean(lat(:)) * pi/180) 1];
+
+
+subplot(2,1,1);
+scatter(lonobs,latobs,10,S)
+caxis(ca);
+colorbar
+hold on, contour(lon,lat,mask,0.5,'k')
+xlim(lon([1 end]))
+ylim(lat([1 end]))
+title('Observations');
+set(gca,'DataAspectRatio',ar,'Layer','top')
+
+subplot(2,1,2);
+pcolor(lon,lat,Sa2), shading flat,colorbar
+hold on
+plot(lonobs,latobs,'k.','MarkerSize',1)
+caxis(ca)
+contour(lon,lat,mask,0.5,'k')
+xlim(lon([1 end]))
+ylim(lat([1 end]))
+title('Analysis');
+set(gca,'DataAspectRatio',ar,'Layer','top')
+
+set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 15 15 ]);
+print -dpng divand_realistic_example.png
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_rms.m b/inst/divand_rms.m
new file mode 100644
index 0000000..a902e48
--- /dev/null
+++ b/inst/divand_rms.m
@@ -0,0 +1,37 @@
+% The RMS error between two variables.
+%
+% r = divand_rms(x,y,norm)
+%
+% Returns rms between x and y (taking norm into account if present)
+
+% Alexander Barth
+function r = divand_rms(x,y,norm)
+
+d = x-y;
+
+if nargin == 3
+ d = d .* sqrt(norm);
+end
+
+m = ~isnan(d);
+r = mean(d(m).^2);
+
+if nargin == 3
+ r = r/mean(norm(m));
+end
+
+r = sqrt(r);
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_simple_example.m b/inst/divand_simple_example.m
new file mode 100644
index 0000000..70cad0b
--- /dev/null
+++ b/inst/divand_simple_example.m
@@ -0,0 +1,68 @@
+% A simple example of divand in 2 dimensions
+% with observations from an analytical function.
+
+
+% observations
+x = rand(75,1);
+y = rand(75,1);
+f = sin(x*6) .* cos(y*6);
+
+% final grid
+[xi,yi] = ndgrid(linspace(0,1,30));
+
+% reference field
+fref = sin(xi*6) .* cos(yi*6);
+
+% all points are valid points
+mask = ones(size(xi));
+
+% this problem has a simple cartesian metric
+% pm is the inverse of the resolution along the 1st dimension
+% pn is the inverse of the resolution along the 2nd dimension
+
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+% correlation length
+len = 0.1;
+
+% signal-to-noise ratio
+lambda = 20;
+
+% fi is the interpolated field
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,len,lambda);
+
+% plotting of results
+subplot(1,2,1);
+pcolor(xi,yi,fref);
+shading flat,colorbar
+ca = caxis;
+set(findobj(gcf,'FontSize',10),'FontSize',16)
+title('Reference field and observation loc.');
+hold on
+plot(x,y,'k.');
+
+subplot(1,2,2);
+hold on
+pcolor(xi,yi,fi);
+shading flat,colorbar
+caxis(ca);
+set(findobj(gcf,'FontSize',10),'FontSize',16)
+title('Interpolated field');
+
+set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 25 12 ]);
+print -dpng divand_simple_example.png
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_solve.m b/inst/divand_solve.m
new file mode 100644
index 0000000..3c6c2de
--- /dev/null
+++ b/inst/divand_solve.m
@@ -0,0 +1,87 @@
+% Solve the variational problem.
+%
+% fi = divand_solve(s,yo)
+%
+% Derive the analysis based on all contraints included in s and using the
+% observations yo
+%
+% Input:
+% s: structure created by divand_factorize
+% yo: value of the observations
+%
+% Output:
+% fi: analyzed field
+
+function [fi,s] = divand_solve(s,yo)
+
+H = s.H;
+sv = s.sv;
+R = s.R;
+% fix me: ignore argument
+yo = s.yo;
+
+
+if s.primal
+ if strcmp(s.inversion,'chol')
+ P = s.P;
+ fpi = P * (H'* (R \ yo(:)));
+ else
+ %H = double(H);
+ HiRyo = H'* (R \ yo(:));
+
+ fun = @(x) s.iB*x + H'*(R\ (H * x));
+
+ if ~s.keepLanczosVectors
+ [fpi,s.flag,s.relres,s.iter] = pcg(fun,HiRyo,s.tol,s.maxit,s.M1,s.M2);
+
+ if s.flag ~= 0
+ warning('divand-noconvergence','pcg did not converge');
+ end
+
+ else
+ %pc = @(x) s.M2 \ (s.M1 \ x);
+ pc = [];
+ x0 = zeros(size(s.iB,1),1);
+ [fpi,Q,T,diag] = conjugategradient(fun,HiRyo,'tol',s.tol,...
+ 'maxit',s.maxit,...
+ 'minit',s.minit,...
+ 'x0',x0,'renorm',1,'pc',pc);
+ s.iter = diag.iter;
+ s.relres = diag.relres;
+ s.P = CovarLanczos(Q,T);
+ end
+ end
+else % dual
+ B = s.B;
+ % fun(x) computes (H B H' + R)*x
+ fun = @(x) H * (B * (H'*x)) + R*x;
+
+ [tmp,flag,relres,iter] = pcg(fun, yo,s.tol,s.maxit,s.funPC);
+
+ if (flag ~= 0)
+ error('divand:pcg', ['Preconditioned conjugate gradients method'...
+ ' did not converge %d %g %g'],flag,relres,iter);
+ end
+
+ fpi = B * (H'*tmp);
+end
+
+
+[fi] = statevector_unpack(sv,fpi);
+
+fi(~s.mask) = NaN;
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/divand_solve_eof.m b/inst/divand_solve_eof.m
new file mode 100644
index 0000000..f0c628a
--- /dev/null
+++ b/inst/divand_solve_eof.m
@@ -0,0 +1,73 @@
+% Solve the variational problem with the contraints from EOFs.
+%
+% [fi,s] = divand_solve_eof(s,yo,EOF_lambda,EOF)
+%
+% Derive the analysis based on all contraints included in s and using the
+% observations yo and EOFs
+%
+% Input:
+% s: structure created by divand_factorize
+% yo: value of the observations
+% EOF_lambda: weight of each EOF (adimentional)
+% EOF: matrix containing the EOFs (units m^(-n/2))
+%
+% Output:
+% fi: analyzed field
+% s.E: scaled EOFs
+
+
+function [fi,s] = divand_solve_eof(s,yo,EOF_lambda,EOF)
+
+H = s.H;
+sv = s.sv;
+R = s.R;
+P = s.P;
+E = s.E;
+
+
+r = size(E,2);
+
+%E = zeros(s.sv.n,r);
+%for i=1:r
+% E(:,i) = statevector_pack(sv,EOF(:,i));
+%end
+
+yo2 = (H'* (R \ yo(:)));
+
+% "classical analysis"
+xa = P * yo2;
+
+% apply Pa to all EOFs
+PaE = P * E;
+
+%keyboard
+
+M = E'*PaE - eye(r);
+
+% analysis with EOFs
+xa = xa - PaE * pinv(M) * (E'*xa);
+
+[fi] = statevector_unpack(sv,xa);
+
+fi(~s.mask) = NaN;
+
+
+% debugging information
+s.cond = cond(M);
+%lam = eig(M);
+%s.lambda_min = min(lam);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/interp_regular.m b/inst/interp_regular.m
new file mode 100644
index 0000000..0e3cf19
--- /dev/null
+++ b/inst/interp_regular.m
@@ -0,0 +1,39 @@
+% Interpolation matrix for a n-dimensional interpolation.
+%
+% [H,out] = interp_regular(mask,x,xi)
+%
+% Creates sparse interpolation matrix H for n-dimensional interpolation problem.
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+% x: cell array with n elements. Every element represents a coordinate of the
+% observations
+% xi: cell array with n elements. Every element represents a coordinate of the
+% final grid on which the observations are interpolated.
+%
+% Output:
+% H: sparse matrix representing the linear interpolation
+% out: 1 for each observation out of the grid, 0 otherwise
+
+function [H,out] = interp_regular(mask,x,xi)
+
+%I = localize_regular_grid(xi,mask,x);
+I = localize_separable_grid(xi,mask,x);
+[H,out] = sparse_interp(mask,I);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/localize_regular_grid.m b/inst/localize_regular_grid.m
new file mode 100644
index 0000000..81efd79
--- /dev/null
+++ b/inst/localize_regular_grid.m
@@ -0,0 +1,67 @@
+% Derive fractional indices on a regular grid.
+%
+% I = localize_regular_grid(xi,mask,x)
+%
+% Derive fractional indices where xi are the points to localize in the
+% regular grid x (constant increment in all dimension).
+% The output I is an n-by-m array where n number of dimensions and m number of
+% observations
+
+function I = localize_regular_grid(xi,mask,x)
+
+x = cat_cell_array(x);
+xi = cat_cell_array(xi);
+
+% m is the number of arbitrarily distributed observations
+tmp = size(xi);
+mi = prod(tmp(1:end-1));
+
+% sz is the size of the grid
+tmp = size(x);
+sz = tmp(1:end-1);
+m = prod(sz);
+
+% n dimension of the problem
+n = tmp(end);
+
+xi = reshape(xi,mi,n);
+
+scale = [1 cumprod(sz(1:end-1))];
+
+A = zeros(n,n);
+
+x0 = x([0:n-1]' * m +1)';
+
+E = eye(n);
+
+I = zeros(n,size(xi,1));
+
+for i=1:n
+ I(i,:) = (xi(:,i) - x0(i));
+
+ % linear index of element (2,1,1,...,:),
+ % (1,2,1,...,:), (1,1,2,...,:) ...
+
+ l = scale*E(:,i) + [0:n-1]' * m +1;
+ A(i,:) = x(l)- x0(:);
+end
+
+
+I = A \ I + 1;
+
+% LocalWords: indices sz
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/localize_separable_grid.m b/inst/localize_separable_grid.m
new file mode 100644
index 0000000..daaf776
--- /dev/null
+++ b/inst/localize_separable_grid.m
@@ -0,0 +1,90 @@
+% Derive fractional indices on a separable grid.
+%
+% I = localize_separable_grid(xi,mask,x)
+%
+% Derive fractional indices where xi are the points to localize in the
+% separable grid x (every dimension in independent on other dimension).
+% The output I is an n-by-m array where n number of dimensions and m number of
+% observations
+
+function I = localize_separable_grid(xi,mask,x)
+
+x = cat_cell_array(x);
+xi = cat_cell_array(xi);
+
+% m is the number of arbitrarily distributed observations
+tmp = size(xi);
+mi = prod(tmp(1:end-1));
+
+% sz is the size of the grid
+tmp = size(x);
+sz = tmp(1:end-1);
+
+if isvector(x)
+ n = 1;
+ mi = length(xi);
+ I = zeros(1,mi);
+ I(1,:) = interp1(x,1:length(x),xi);
+else
+ % n dimension of the problem
+ n = tmp(end);
+ m = prod(sz);
+
+ xi = reshape(xi,mi,n);
+ x = reshape(x,[m n]);
+
+ IJ = cell(n,1);
+ vi = cell(n,1);
+ X = cell(n,1);
+ XI = cell(n,1);
+ I = zeros(n,mi);
+
+ for i=1:n
+ vi{i} = 1:sz(i);
+ X{i} = reshape(x(:,i),sz);
+ XI{i} = xi(:,i);
+ end
+
+ [IJ{:}] = ndgrid(vi{:});
+
+
+ for i=1:n
+ I(i,:) = interpn(X{:},IJ{i},XI{:});
+ end
+
+end
+
+
+% handle rounding errors
+% snap to domain bounding box if difference does not exceeds tol
+tol = 50*eps;
+
+for i=1:n
+ % upper bound
+ ind = sz(i) < I(i,:) & I(i,:) <= sz(i) + tol;
+ I(i,ind) = sz(i);
+
+ % lower bound
+ ind = 1 < I(i,:) & I(i,:) <= 1 + tol;
+ I(i,ind) = 1;
+end
+
+
+
+
+% LocalWords: indices sz tol
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/mtimescorr.m b/inst/mtimescorr.m
new file mode 100644
index 0000000..5eb7b21
--- /dev/null
+++ b/inst/mtimescorr.m
@@ -0,0 +1,46 @@
+% Product between a Gaussian covariance matrix and a vector.
+%
+% p = mtimescorr(xi,len,b)
+%
+% Compute the product between a Gaussian covariance with a length scale len and
+% grid points located in xi and the vector b. The covariance matrix has a unit
+% variance.
+
+function p = mtimescorr(xi,len,b)
+
+%'mat'
+m = size(xi,1);
+p = zeros(size(b));
+
+% for j=1:m
+% for i=1:m
+% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 );
+% p(i) = p(i) + exp(-d2 / self.len^2) * b(j);
+% end
+% end
+%
+
+
+for i=1:m
+ X = repmat(xi(i,:),[m 1]);
+ d2 = sum( (X - xi).^2, 2) ;
+ ed = exp(-d2 / len^2);
+ p(i,:) = ed' * b;
+end
+
+% LocalWords: mtimescorr len
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/private/geomean.m b/inst/private/geomean.m
new file mode 100644
index 0000000..bbc59f7
--- /dev/null
+++ b/inst/private/geomean.m
@@ -0,0 +1,19 @@
+% for compatibility with octave
+
+function m = geomean(x,varargin)
+
+m = prod(x,varargin{:})^(1/length(x));
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/sparse_diag.m b/inst/sparse_diag.m
new file mode 100644
index 0000000..8cd1c09
--- /dev/null
+++ b/inst/sparse_diag.m
@@ -0,0 +1,26 @@
+% Create diagonal sparse matrix.
+%
+% s = sparse_diag(d)
+%
+% Create diagonal sparse matrix s based on the diagonal elements d. The
+% parameter d is transformed into a vector.
+
+function s = sparse_diag(d)
+
+s = spdiags(d(:),0,numel(d),numel(d));
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/sparse_diff.m b/inst/sparse_diff.m
new file mode 100644
index 0000000..b09cd6d
--- /dev/null
+++ b/inst/sparse_diff.m
@@ -0,0 +1,81 @@
+% Sparse operator for differentiation.
+%
+% diffx = sparse_diff(sz1,m,cyclic)
+%
+% Sparse operator for differentiation along dimension m for "collapsed" matrix
+% of the size sz1.
+%
+% Input:
+% sz1: size of rhs
+% m: dimension to differentiate
+% cyclic: true if domain is cyclic along dimension m. False is the
+% default value
+
+function diffx = sparse_diff(sz1,m,cyclic)
+
+if nargin == 2
+ cyclic = false;
+end
+
+n1 = prod(sz1);
+sz2 = sz1;
+
+if ~cyclic
+ sz2(m) = sz2(m)-1;
+end
+
+n2 = prod(sz2);
+
+n = length(sz1);
+
+vi = cell(1,n);
+for i=1:n
+ vi{i} = 1:sz2(i);
+end
+
+IJ = cell(n,1);
+
+[IJ{:}] = ndgrid(vi{:});
+
+for i=1:n
+ IJ{i} = IJ{i}(:);
+end
+
+L1 = [1:n2]';
+
+L2 = sub2ind(sz1,IJ{:});
+one = ones(size(L1));
+
+IJ{m} = IJ{m} + 1;
+
+if cyclic
+ IJ{m} = mod(IJ{m}-1,sz1(m))+1;
+end
+
+
+L2o = sub2ind(sz1,IJ{:});
+
+diffx = sparse( ...
+ [L1; L1; ], ...
+ [L2; L2o; ], ...
+ [-one; one ], n2 , n1 );
+
+
+
+% Copyright (C) 2009,2012 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
+
+% LocalWords: diffx sz
diff --git a/inst/sparse_gradient.m b/inst/sparse_gradient.m
new file mode 100644
index 0000000..ff6412a
--- /dev/null
+++ b/inst/sparse_gradient.m
@@ -0,0 +1,59 @@
+% Sparse operator for a gradient.
+%
+% [Dx1,Dx2,...,Dxn] = sparse_gradient(mask,pmn)
+%
+% Form the gradient using finite differences in all n-dimensions
+%
+% Input:
+% mask: binary mask delimiting the domain. 1 is inside and 0 outside.
+% For oceanographic application, this is the land-sea mask.
+%
+% pmn: scale factor of the grid.
+%
+% Output:
+% Dx1,Dx2,...,Dxn: sparce matrix represeting a gradient along
+% different dimensions
+
+function varargout = sparse_gradient(mask,pmn,iscyclic)
+
+H = sparse_pack(mask);
+
+sz = size(mask);
+n = size(pmn,ndims(pmn));
+Pmn = reshape(pmn,[prod(sz) n]);
+
+if nargin == 2
+ iscyclic = zeros(n,1);
+end
+
+for i=1:n
+ % staggering operator
+ S = sparse_stagger(sz,i,iscyclic(i));
+
+ % mask for staggered variable
+ m = (S * mask(:) == 1);
+
+ d = m .* (S * Pmn(:,i));
+
+ varargout{i} = sparse_pack(m) * sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)) * H';
+% varargout{i} = sparse_diag(d) * sparse_diff(sz,i) * H';
+% size(varargout{i})
+end
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/sparse_interp.m b/inst/sparse_interp.m
new file mode 100644
index 0000000..ffd8ebd
--- /dev/null
+++ b/inst/sparse_interp.m
@@ -0,0 +1,127 @@
+% Create a sparse interpolation matrix.
+%
+% [H,out] = sparse_interp(mask,I)
+%
+% Create interpolation matrix from mask and fractional indexes I
+%
+% Input:
+% mask: 0 invalid and 1 valid points (n-dimensional array)
+% I: fractional indexes (2-dim array n by mi, where mi is the number of points to interpolate)
+% Ouput:
+% H: sparse matrix with interpolation coefficients
+% out: true if value outside of grid
+
+function [H,out] = sparse_interp(mask,I,iscyclic)
+
+if ndims(mask) ~= size(I,1) && ~(isvector(mask) && size(I,1) == 1)
+ error('sparse_interp: inconsistent arguments')
+end
+
+if isvector(mask)
+ sz = numel(mask);
+ mask = reshape(mask,[sz 1]);
+else
+ sz = size(mask);
+end
+
+% n dimension of the problem
+n = size(I,1);
+m = prod(sz);
+
+if nargin == 2
+ iscyclic = zeros(1,n);
+end
+
+
+% mi is the number of arbitrarly distributed observations
+mi = size(I,2);
+
+% handle cyclic dimensions
+for i = 1:n
+ if iscyclic(i)
+ % bring I(i,:) inside the interval [1 sz(i)+1[
+ % since the i-th dimension is cyclic
+
+ I(i,:) = mod(I(i,:)-1,sz(i))+1;
+ end
+end
+
+
+scale = [1 cumprod(sz(1:end-1))];
+
+% integer index
+ind = floor(I);
+inside = true(1,mi);
+
+for i = 1:n
+ if ~iscyclic(i)
+ % make a range check only for non-cyclic dimension
+
+ % handle border cases
+ p = find(I(i,:) == sz(i));
+ ind(i,p) = sz(i)-1;
+
+ inside = inside & 1 <= ind(i,:) & ind(i,:) < sz(i);
+ end
+end
+
+%whos ind
+
+% interpolation coefficient
+coeff(:,:,2) = I - ind;
+coeff(:,:,1) = 1 - coeff(:,:,2);
+
+% consider now only points inside
+iind = find(inside);
+coeff2 = coeff(:,iind,:);
+ind2 = ind(:,iind);
+mip = length(iind); % number of obs. inside
+
+si = repmat(1:mip,[2^n 1]);
+sj = ones(2^n, mip);
+ss = ones(2^n, mip);
+
+% loop over all corner of hypercube
+for i=1:2^n
+
+ % loop over all dimensions
+ for j=1:n
+ bit = bitget(i,j);
+
+ % the j-index of the i-th corner has the index ip
+ % this index ip is zero-based
+ ip = (ind2(j,:) + bit - 1);
+
+ % ip must be [0 and sz(j)-1] (zero-based)
+ % we know aleady that the point is inside the domain
+ % so, if it is outside this range then it is because of periodicity
+ ip = mod(ip,sz(j));
+
+ sj(i,:) = sj(i,:) + scale(j) * ip;
+ ss(i,:) = ss(i,:) .* coeff2(j,:,bit + 1);
+ end
+end
+
+% sj must refer to a valid point or its interpolation coefficient
+% must be zero
+
+inside(iind) = all(mask(sj) | ss == 0,1);
+H = sparse(iind(si(:)),sj(:),ss(:),mi,m);
+
+out = ~inside;
+% LocalWords: interp
+
+% Copyright (C) 2004 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/sparse_pack.m b/inst/sparse_pack.m
new file mode 100644
index 0000000..e8d1888
--- /dev/null
+++ b/inst/sparse_pack.m
@@ -0,0 +1,35 @@
+% Create a sparse matrix which packs an array under the control of a mask.
+%
+% H = sparse_pack(mask)
+%
+% Create a sparse matrix H which packs an array. It selects the elements where
+% mask is true.
+
+function H = sparse_pack(mask)
+
+
+j = find(mask)';
+m = length(j);
+i = 1:m;
+s = ones(1,m);
+
+n = numel(mask);
+
+%whos i j s
+H = sparse(i,j,s,m,n);
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/sparse_shift.m b/inst/sparse_shift.m
new file mode 100644
index 0000000..b5e9e3c
--- /dev/null
+++ b/inst/sparse_shift.m
@@ -0,0 +1,75 @@
+% Sparse operator shifting a field in a given dimension.
+%
+% function S = sparse_shift(sz1,m,cyclic)
+%
+% Sparse operator shifting a field in the dimension m. The field is a
+% "collapsed" matrix of the size sz1.
+%
+% Input:
+% sz1: size of rhs
+% m: dimension to shift
+% cyclic: true if domain is cyclic along dimension m. False is the
+% default value
+
+function S = sparse_shift(sz1,m,cyclic)
+
+if nargin == 2
+ cyclic = false;
+end
+
+n1 = prod(sz1);
+
+sz2 = sz1;
+
+if ~cyclic
+ sz2(m) = sz2(m)-1;
+end
+
+n2 = prod(sz2);
+
+n = length(sz1);
+
+vi = cell(1,n);
+for i=1:n
+ vi{i} = 1:sz2(i);
+end
+
+IJ = cell(n,1);
+
+[IJ{:}] = ndgrid(vi{:});
+
+for i=1:n
+ IJ{i} = IJ{i}(:);
+end
+
+L1 = [1:n2]';
+one = ones(size(L1));
+
+IJ{m} = IJ{m} + 1;
+
+if cyclic
+ IJ{m} = mod(IJ{m}-1,sz1(m))+1;
+end
+
+
+L2 = sub2ind(sz1,IJ{:});
+
+S = sparse(L1,L2,one,n2,n1);
+
+
+
+% Copyright (C) 2012 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/sparse_stagger.m b/inst/sparse_stagger.m
new file mode 100644
index 0000000..6a25473
--- /dev/null
+++ b/inst/sparse_stagger.m
@@ -0,0 +1,78 @@
+% Sparse operator for staggering.
+%
+% S = sparse_stagger(sz1,m,cyclic)
+%
+% Create a sparse operator for staggering a field in dimension m.
+% The field is a "collapsed" matrix of the size sz1.
+%
+% Input:
+% sz1: size of rhs
+% m: dimension to stagger
+% cyclic: true if domain is cyclic along dimension m. False is the
+% default value
+
+function S = sparse_stagger(sz1,m,cyclic)
+
+
+if nargin == 2
+ cyclic = false;
+end
+
+n1 = prod(sz1);
+
+sz2 = sz1;
+
+if ~cyclic
+ sz2(m) = sz2(m)-1;
+end
+n2 = prod(sz2);
+
+n = length(sz1);
+
+for i=1:n
+ vi{i} = 1:sz2(i);
+end
+
+IJ = cell(n,1);
+
+[IJ{:}] = ndgrid(vi{:});
+
+for i=1:n
+ IJ{i} = IJ{i}(:);
+end
+
+L1 = [1:n2]';
+
+L2 = sub2ind(sz1,IJ{:});
+v = ones(size(L1))/2;
+
+IJ{m} = IJ{m} + 1;
+
+if cyclic
+ IJ{m} = mod(IJ{m}-1,sz1(m))+1;
+end
+
+L2o = sub2ind(sz1,IJ{:});
+
+S = sparse( ...
+ [L1; L1; ], ...
+ [L2; L2o; ], ...
+ [v; v ], n2 , n1 );
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/sparse_trim.m b/inst/sparse_trim.m
new file mode 100644
index 0000000..fff9d9a
--- /dev/null
+++ b/inst/sparse_trim.m
@@ -0,0 +1,61 @@
+% Sparse operator for trimming.
+%
+% T = sparse_trim(sz1,m)
+%
+% Create a sparse operator which trim first and last row (or column) in
+% The field is a "collapsed" matrix of the size sz1.
+%
+% Input:
+% sz1: size of rhs
+% m: dimension to trim
+
+function T = sparse_trim(sz1,m)
+
+n1 = prod(sz1);
+
+sz2 = sz1;
+sz2(m) = sz2(m)-2;
+n2 = prod(sz2);
+
+n = length(sz1);
+
+% L1
+
+for i=1:n
+ vi{i} = 1:sz2(i);
+end
+
+IJ = cell(n,1);
+
+[IJ{:}] = ndgrid(vi{:});
+
+for i=1:n
+ IJ{i} = IJ{i}(:);
+end
+
+L1 = sub2ind(sz2,IJ{:});
+
+IJ{m}=IJ{m}+1;
+
+L2 = sub2ind(sz1,IJ{:});
+
+one = ones(size(L1));
+
+T = sparse(L1, L2, one, n2, n1);
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
diff --git a/inst/statevector_init.m b/inst/statevector_init.m
new file mode 100644
index 0000000..3c0d24e
--- /dev/null
+++ b/inst/statevector_init.m
@@ -0,0 +1,74 @@
+% Initialize structure for packing and unpacking given their mask.
+%
+% s = statevector_init(mask1, mask2, ...)
+%
+% Initialize structure for packing and unpacking
+% multiple variables given their corresponding land-sea mask.
+%
+% Input:
+% mask1, mask2,...: land-sea mask for variable 1,2,... Sea grid points correspond to one and land grid points to zero.
+% Every mask can have a different shape.
+%
+% Output:
+% s: structure to be used with statevector_pack and statevector_unpack.
+%
+% Note:
+% see also statevector_pack, statevector_unpack
+
+% Author: Alexander Barth, 2009 <a.barth at ulg.ac.be>
+% License: GPL 2 or later
+
+function s = statevector_init(varargin)
+
+s.nvar = nargin;
+
+
+for i=1:nargin
+ mask = varargin{i};
+ s.mask{i} = mask;
+ s.numels(i) = sum(mask(:) == 1);
+ s.numels_all(i) = numel(mask);
+ s.size{i} = size(mask);
+end
+
+s.ind = [0 cumsum(s.numels)];
+
+s.n = s.ind(end);
+
+%!test
+%! mask = rand(10,10) > .5;
+%! mask_u = rand(9,10) > .5;
+%! mask_v = rand(10,9) > .5;
+%! sv = statevector_init(mask,mask_u,mask_v);
+%! var = rand(10,10);
+%! var(mask==0) = 0;
+%! var_u = rand(9,10);
+%! var_u(mask_u==0) = 0;
+%! var(mask==0) = 0;
+%! var_v = rand(10,9);
+%! var_v(mask_v==0) = 0;
+%! [E] = statevector_pack(sv,var,var_u,var_v);
+%! [Ezeta2,Eu2,Ev2] = statevector_unpack(sv,E);
+%! assert(Ezeta2,var)
+%! assert(Eu2,var_u)
+%! assert(Ev2,var_v)
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
+
+% LocalWords: statevector init GPL
diff --git a/inst/statevector_pack.m b/inst/statevector_pack.m
new file mode 100644
index 0000000..956d371
--- /dev/null
+++ b/inst/statevector_pack.m
@@ -0,0 +1,75 @@
+% Pack a series of variables into a vector under the control of a mask.
+%
+% x = statevector_pack(s,var1, var2, ...)
+%
+% Pack the different variables var1, var2, ... into the vector x.
+% Only sea grid points are retained.
+%
+% Input:
+% s: structure generated by statevector_init.
+% var1, var2,...: variables to pack (with the same shape as the corresponding masks).
+%
+% Output:
+% x: vector of the packed elements. The size of this vector is the number of elements of all masks equal to 1.
+%
+% Notes:
+% If var1, var2, ... have an additional trailing dimension, then this dimension is assumed
+% to represent the different ensemble members. In this case x is a matrix and its last dimension
+% is the number of ensemble members.
+
+% Author: Alexander Barth, 2009 <a.barth at ulg.ac.be>
+% License: GPL 2 or later
+
+function x = statevector_pack(varargin)
+
+s = varargin{1};
+
+if isvector(varargin{2}) && isvector(s.mask{1})
+ % handle odd case when varargin{1} is a vector of size n x 1 and
+ % s.mask{1} is 1 x n
+ k = 1;
+else
+ k = size(varargin{2},my_ndims(s.mask{1})+1);
+end
+
+x = zeros(s.n,k);
+
+for i=1:s.nvar
+ tmp = reshape(varargin{i+1},s.numels_all(i),k);
+
+ ind = find(s.mask{i}==1);
+
+ x(s.ind(i)+1:s.ind(i+1),:) = tmp(ind,:);
+end
+
+
+function d = my_ndims(v)
+ if isvector(v)
+ d = 1;
+ else
+ d = ndims(v);
+ end
+
+
+
+
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
+
+% LocalWords: statevector init GPL
diff --git a/inst/statevector_unpack.m b/inst/statevector_unpack.m
new file mode 100644
index 0000000..e1d492e
--- /dev/null
+++ b/inst/statevector_unpack.m
@@ -0,0 +1,64 @@
+% Unpack a vector into different variables under the control of a mask.
+%
+% [var1, var2, ...] = statevector_unpack(s,x)
+% [var1, var2, ...] = statevector_unpack(s,x,fillvalue)
+%
+% Unpack the vector x into the different variables var1, var2, ...
+%
+% Input:
+% s: structure generated by statevector_init.
+% x: vector of the packed elements. The size of this vector is the number of elements equal to 1
+% in all masks.
+%
+% Optional input parameter:
+% fillvalue: The value to fill in var1, var2,... where the masks correspond to a land grid point. The default is zero.
+%
+% Output:
+% var1, var2,...: unpacked variables.
+%
+% Notes:
+% If x is a matrix, then the second dimension is assumed
+% to represent the different ensemble members. In this case,
+% var1, var2, ... have also an additional trailing dimension.
+
+% Author: Alexander Barth, 2009 <a.barth at ulg.ac.be>
+% License: GPL 2 or later
+
+function varargout = statevector_unpack(s,x,fillvalue)
+
+if (nargin == 2)
+ fillvalue = 0;
+end
+
+k = size(x,2);
+
+for i=1:s.nvar
+ v = zeros(s.numels_all(i),k);
+ v(:) = fillvalue;
+
+ ind = find(s.mask{i}==1);
+
+ v(ind,:) = x(s.ind(i)+1:s.ind(i+1),:);
+
+ varargout{i} = reshape(v,[s.size{i} k]);
+end
+
+
+
+% Copyright (C) 2009 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
+
+
+% LocalWords: statevector fillvalue init GPL varargout
diff --git a/inst/test_1dvar.m b/inst/test_1dvar.m
new file mode 100644
index 0000000..ef73fab
--- /dev/null
+++ b/inst/test_1dvar.m
@@ -0,0 +1,30 @@
+% Testing divand in 1 dimension.
+
+% grid of background field
+[xi] = linspace(0,1,30)';
+
+x = [.4 .6]';
+f = [.4 .6]';
+
+mask = ones(size(xi));
+mask([1 end]) = 0;
+
+pm = ones(size(xi)) / (xi(2)-xi(1));
+
+[fi,err] = divand(mask,pm,{xi},{x},f,.1,2);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar.m b/inst/test_2dvar.m
new file mode 100644
index 0000000..068fc4f
--- /dev/null
+++ b/inst/test_2dvar.m
@@ -0,0 +1,56 @@
+% Testing divand in 2 dimensions.
+
+%try
+% grid of background field
+[xi,yi] = ndgrid(linspace(0,1,30));
+fi_ref = sin( xi*6 ) .* cos( yi*6);
+
+% grid of observations
+[x,y] = ndgrid(linspace(eps,1-eps,20));
+x = x(:);
+y = y(:);
+
+f = sin( x*6 ) .* cos( y*6);
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+%fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,.1,20);
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,.1,20,'factorize',0);
+rms = divand_rms(fi_ref,fi);
+
+if (rms > 0.005)
+ error('unexpected large difference with reference field');
+end
+
+%fprintf('OK (rms=%g)\n',rms);
+
+
+fi_dual = divand(mask,{pm,pn},{xi,yi},{x,y},f,.1,20,'primal',0);
+rms = divand_rms(fi_dual,fi);
+
+%fprintf(1,'Testing dual 2D-optimal variational inverse: ');
+
+if (rms > 1e-6)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',rms);
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_adv.m b/inst/test_2dvar_adv.m
new file mode 100644
index 0000000..c30642c
--- /dev/null
+++ b/inst/test_2dvar_adv.m
@@ -0,0 +1,42 @@
+% Testing divand in 2 dimensions with advection.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(-1,1,30));
+
+x = .4;
+y = .4;
+f =1;
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+u = ones(size(xi));
+v = u;
+
+a = 5;
+u = a*yi;
+v = -a*xi;
+
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,.2,200,'velocity',{u,v});
+
+if abs( fi(18,24) - 0.89935) > 1e-3
+ error('unexpected large difference');
+end
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_check.m b/inst/test_2dvar_check.m
new file mode 100644
index 0000000..64ee289
--- /dev/null
+++ b/inst/test_2dvar_check.m
@@ -0,0 +1,64 @@
+% Testing divand in 2 dimensions with independent verification.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(0,1,10));
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+% grid of observations
+[x,y] = ndgrid(linspace(eps,1-eps,10));
+x = x(:);
+y = y(:);
+v = sin( x*6 ) .* cos( y*6);
+
+
+lenx = .15;
+leny = .15;
+
+lambda = 20;
+
+[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,{lenx,leny},lambda,'diagnostics',1,'primal',1);
+
+
+
+iR = inv(full(s.R));
+iB = full(s.iB);
+H = full(s.H);
+sv = s.sv;
+
+iP = iB + H'*iR*H;
+
+P = inv(iP);
+
+xa2 = P* (H'*iR*v(:));
+
+[fi2] = statevector_unpack(sv,xa2);
+fi2(~s.mask) = NaN;
+
+rms_diff = [];
+rms_diff(end+1) = divand_rms(va,fi2);
+rms_diff(end+1) = divand_rms(diag(P),err(:));
+
+
+if any(rms_diff > 1e-6)
+ error('unexpected large difference');
+end
+
+fprintf('(max rms=%g) ',max(rms_diff));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_check_correrr.m b/inst/test_2dvar_check_correrr.m
new file mode 100644
index 0000000..30876de
--- /dev/null
+++ b/inst/test_2dvar_check_correrr.m
@@ -0,0 +1,70 @@
+% Testing divand in 2 dimensions with correlated errors.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(0,1,10));
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+% grid of observations
+[x,y] = ndgrid(linspace(eps,1-eps,3));
+x = x(:);
+y = y(:);
+v = sin( x*6 ) .* cos( y*6);
+
+xx = repmat(x,[1 numel(x)]);
+yy = repmat(y,[1 numel(y)]);
+
+len = 0.1;
+variance = 0.1;
+R = variance * exp(- ((xx - xx').^2 + (yy - yy').^2 ) / len^2);
+RC = CovarParam({x,y},variance,len);
+
+len = .1;
+lambda = RC;
+
+[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',0);
+
+iR = inv(full(R));
+iB = full(s.iB);
+H = full(s.H);
+sv = s.sv;
+
+iP = iB + H'*iR*H;
+
+Pa = inv(iP);
+
+xa2 = Pa* (H'*iR*v(:));
+
+[fi2] = statevector_unpack(sv,xa2);
+fi2(~s.mask) = NaN;
+
+rms_diff = divand_rms(va,fi2);
+
+if rms_diff > 1e-6
+ disp('FAILED');
+end
+
+rms_diff = divand_rms(err(:),diag(Pa));
+
+if rms_diff > 1e-6
+ error('unexpected large difference');
+end
+
+fprintf('(max rms=%g) ',max(rms_diff));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_constrain.m b/inst/test_2dvar_constrain.m
new file mode 100644
index 0000000..f3b6552
--- /dev/null
+++ b/inst/test_2dvar_constrain.m
@@ -0,0 +1,84 @@
+% Testing divand in 2 dimensions with a custom constrain.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(0,1,10));
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+% grid of observations
+[x,y] = ndgrid(linspace(eps,1-eps,10));
+x = x(:);
+y = y(:);
+v = sin( x*6 ) .* cos( y*6);
+
+
+len = 3;
+lambda = 20;
+n = sum(mask(:));
+
+% constrain
+fun = @(x) mean(x,1);
+funt = @(x) repmat(x,[n 1])/n;
+
+H = ones(1,n)/n;
+c.H = MatFun([1 n],fun,funt);
+%c.H = H;
+c.R = 0.001;
+c.yo = 1;
+
+% bug
+%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0,'constraint',c,'inversion','pcg');
+
+[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,...
+ 'diagnostics',1,'primal',1,'factorize',0,'inversion','pcg',...
+ 'tol',1e-6,'minit',100,'keepLanczosVectors',1);
+
+
+% ok
+%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0);
+
+
+
+
+iR = inv(full(s.R));
+iB = full(s.iB);
+H = double(s.H);
+yo = s.yo;
+sv = s.sv;
+
+iB = iB + H'*iR*H;
+
+Pa = inv(iB);
+
+xa2 = Pa* (H'*iR*yo);
+
+[fi2] = statevector_unpack(sv,xa2);
+fi2(~s.mask) = NaN;
+
+rms_diff = [];
+rms_diff(end+1) = sqrt(mean((va(:) - fi2(:)).^2));
+rms_diff(end+1) = sqrt(mean((diag(Pa) - err(:)).^2));
+
+
+if any(rms_diff > 1e-4)
+ error('unexpected large difference');
+end
+
+fprintf('(max rms=%g) ',max(rms_diff));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_cyclic.m b/inst/test_2dvar_cyclic.m
new file mode 100644
index 0000000..4d87977
--- /dev/null
+++ b/inst/test_2dvar_cyclic.m
@@ -0,0 +1,97 @@
+% Testing divand in 2 dimensions in a cyclic domain.
+
+% grid of background field
+n = 30;
+
+[xi,yi] = ndgrid(0:n-1);
+
+x = xi(2,15);
+y = yi(2,15);
+f = 1;
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+
+
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0]);
+rms_diff = [];
+
+rms_diff(end+1) = divand_rms(fi(end,:),fi(4,:));
+
+
+x = xi(2,2);
+y = yi(2,2);
+
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n n]);
+
+rms_diff(end+1) = abs(fi(4,4) - fi(end,end));
+
+
+fi2 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.]);
+fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[30.5; 15.]);
+rms_diff(end+1) = divand_rms(fi2,fi3);
+
+
+
+fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.]);
+fi4 = cat(1,fi3(11:end,:),fi3(1:10,:));
+
+rms_diff(end+1) = divand_rms(fi2,fi4);
+
+
+
+
+
+% test with mask
+mask(11,4) = 0;
+mask(3:15,20) = 0;
+
+fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.]);
+fi4 = cat(1,fi3(11:end,:),fi3(1:10,:));
+
+mask = cat(1,mask(11:end,:),mask(1:10,:));
+
+fi2 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.]);
+
+
+rms_diff(end+1) = divand_rms(fi2,fi4);
+rms_diff(end+1) = divand_rms(isnan(fi2),isnan(fi4));
+
+% advection
+
+a = 5;
+u = a*cos(2*pi*yi/n);
+v = a*cos(2*pi*xi/n);
+
+mask = ones(size(xi));
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.],'velocity',{u,v});
+fi2 = cat(1,fi(11:end,:),fi(1:10,:));
+
+ut = cat(1,u(11:end,:),u(1:10,:));
+vt = cat(1,v(11:end,:),v(1:10,:));
+fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.],'velocity',{ut,vt});
+
+rms_diff(end+1) = divand_rms(fi2,fi3);
+
+if any(rms_diff > 1e-10) || any(isnan(rms_diff))
+ error('unexpected large difference');
+end
+
+fprintf('(max rms=%g) ',max(rms_diff));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_eof_check.m b/inst/test_2dvar_eof_check.m
new file mode 100644
index 0000000..499b289
--- /dev/null
+++ b/inst/test_2dvar_eof_check.m
@@ -0,0 +1,97 @@
+% Testing divand in 2 dimensions with EOF constraints.
+
+[xi,yi] = ndgrid(linspace(0,1,30));
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+
+x = .5;
+y = .5;
+v = 1;
+
+len = .1;
+lambda = 10;
+
+clear U2
+U2(:,1) = mask(:);
+U2(:,2) = xi(:);
+U2(:,3) = yi(:);
+
+%U2 = randn(numel(mask),10);
+r = size(U2,2);
+[U3,S3,V3] = svds(U2,r);
+U3 = reshape(U3,[size(mask) r]);
+
+[va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1);
+beta = 10;
+
+[va_eof,s_eof] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'EOF',U3,'EOF_scaling',beta*ones(r,1),'diagnostics',1);
+
+[Jb,Jo,Jeof] = divand_diagnose(s,va_eof,v);
+
+%assert(s_eof.Jb,Jb,1e-6) -> should no longer be the case since s.iB is modified to ensure that B - EE is positive defined
+assert(abs(s_eof.Jo-Jo) < 1e-6)
+
+
+
+s = s_eof;
+
+iR = inv(full(s.R));
+iB = s.iB;
+H = s.H;
+E = s.E;
+sv = s.sv;
+
+iP = iB + H'*iR*H;
+
+Pa = inv(iP - E*E');
+
+xa2 = Pa* (H'*iR*v(:));
+
+[fi2] = statevector_unpack(sv,xa2);
+fi2(~s.mask) = NaN;
+
+rms_diff = divand_rms(va_eof,fi2);
+
+if rms_diff > 1e-6
+ error('unexpected large difference');
+end
+
+Pf = inv(iB - E*E');
+R = inv(iR);
+xa3 = Pf*H' * inv(H*Pf*H' + R) * v;
+
+[fi3] = statevector_unpack(sv,xa3);
+fi3(~s.mask) = NaN;
+
+rms_diff = divand_rms(va_eof,fi3);
+
+if rms_diff > 1e-6
+ error('unexpected large difference');
+end
+
+[U,L] = eig(Pf);
+L = real(diag(L));
+
+if any(L < 0)
+ error('some eingenvalues are negative');
+end
+
+Pfvar = inv(iB);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_lenxy.m b/inst/test_2dvar_lenxy.m
new file mode 100644
index 0000000..7dd24f7
--- /dev/null
+++ b/inst/test_2dvar_lenxy.m
@@ -0,0 +1,54 @@
+% Testing divand in 2 dimensions with lenx /= leny.
+
+n2 = 15;
+% grid of background field
+[xi,yi] = ndgrid(linspace(-1,1,2*n2+1));
+
+x = 0;
+y = 0;
+f = 1;
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+u = ones(size(xi));
+v = u;
+
+a = 5;
+u = a*yi;
+v = -a*xi;
+
+L = 0.2 * ones(size(xi));
+%L = 0.2 * ones(size(xi));
+L = 0.4 * (2 - xi);
+L = ones([size(mask) 2]);
+% must be L(:,:,1) < L(:,:,2)
+
+L(:,:,1) = 0.2;
+L(:,:,2) = 0.5;
+
+Lx = 0.2 * ones(size(mask));
+Ly = 0.5 * ones(size(mask));
+
+
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,{Lx,Ly},2);
+
+if fi(n2,1) < fi(1,n2)
+ error('unexpected values');
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_pcg.m b/inst/test_2dvar_pcg.m
new file mode 100644
index 0000000..577682b
--- /dev/null
+++ b/inst/test_2dvar_pcg.m
@@ -0,0 +1,65 @@
+% Testing divand in 2 dimensions with pcg solver.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(0,1,20));
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+% grid of observations
+[x,y] = ndgrid(linspace(eps,1-eps,10));
+x = x(:);
+y = y(:);
+v = sin( x*6 ) .* cos( y*6);
+
+
+len = 0.1;
+lambda = 1;
+n = sum(mask(:));
+
+% constrain
+fun = @(x) mean(x,1);
+funt = @(x) repmat(x,[n 1])/n;
+
+H = ones(1,n)/n;
+c.H = MatFun([1 n],fun,funt);
+%c.H = H;
+c.R = 0.001;
+c.yo = 1;
+
+% bug
+%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0,'constraint',c,'inversion','pcg');
+
+[va0,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,...
+ 'diagnostics',1,'primal',1);
+
+[va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,...
+ 'diagnostics',1,'primal',1,'inversion','pcg',...
+ 'tol',1e-6,'maxit',10000);
+
+
+rms_diff = [];
+rms_diff(end+1) = sqrt(mean((va(:) - va0(:)).^2));
+
+
+if any(rms_diff > 1e-4)
+ error('unexpected large difference');
+end
+
+fprintf('(max rms=%g) ',max(rms_diff));
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_2dvar_rellen.m b/inst/test_2dvar_rellen.m
new file mode 100644
index 0000000..3f8e70f
--- /dev/null
+++ b/inst/test_2dvar_rellen.m
@@ -0,0 +1,47 @@
+% Testing divand in 2 dimensions with relative correlation length.
+
+% grid of background field
+[xi,yi] = ndgrid(linspace(-1,1,30));
+
+x = 0;
+y = 0;
+f = 1;
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
+pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
+
+
+
+u = ones(size(xi));
+v = u;
+
+a = 5;
+u = a*yi;
+v = -a*xi;
+
+L = 0.2 * ones(size(xi));
+%L = 0.2 * ones(size(xi));
+L = 0.4 * (2 - xi);
+
+fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,L,200);
+
+
+if any(fi(1,:) < fi(end,:))
+ error('unexpected values');
+end
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_3dvar.m b/inst/test_3dvar.m
new file mode 100644
index 0000000..94bd3d3
--- /dev/null
+++ b/inst/test_3dvar.m
@@ -0,0 +1,58 @@
+% Testing divand in 3 dimensions.
+
+% grid of background field
+[xi,yi,zi] = ndgrid(linspace(0,1.,15));
+fi_ref = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+
+% grid of observations
+[x,y,z] = ndgrid(linspace(eps,1-eps,10));
+x = x(:);
+y = y(:);
+z = z(:);
+
+% reference field
+f = sin(6*x) .* cos(6*y) .* sin(6*z);
+
+% all points are valid points
+mask = ones(size(xi));
+
+% this problem has a simple cartesian metric
+% pm is the inverse of the resolution along the 1st dimension
+% pn is the inverse of the resolution along the 2nd dimension
+% po is the inverse of the resolution along the 3rd dimension
+pm = ones(size(xi)) / (xi(2,1,1)-xi(1,1,1));
+pn = ones(size(xi)) / (yi(1,2,1)-yi(1,1,1));
+po = ones(size(xi)) / (zi(1,1,2)-zi(1,1,1));
+
+% correlation length
+len = 0.1;
+
+% signal-to-noise ratio
+lambda = 100;
+
+% fi is the interpolated field
+fi = divand(mask,{pm,pn,po},{xi,yi,zi},{x,y,z},f,len,lambda);
+
+% compute RMS to background field
+rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+if (rms > 0.04)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',rms);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_3dvar_large_stacked.m b/inst/test_3dvar_large_stacked.m
new file mode 100644
index 0000000..62f14b7
--- /dev/null
+++ b/inst/test_3dvar_large_stacked.m
@@ -0,0 +1,56 @@
+% Testing divand in 3 dimensions without correlation in the 3rd dimension
+% (vertically stacked).
+
+kmax = 4;
+
+% grid of background field
+[xi,yi,zi] = ndgrid(linspace(0,1.,30),linspace(0,1.,30),linspace(0,1.,10));
+fi_ref = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+
+% grid of observations
+[x,y,z] = ndgrid(linspace(eps,1-eps,10),linspace(eps,1-eps,10),linspace(eps,1-eps,10));
+x = reshape(x,100,10);
+y = reshape(y,100,10);
+z = reshape(z,100,10);
+
+
+f = sin(6*x) .* cos(6*y) .* sin(6*z);
+
+
+mask = ones(size(xi));
+pm = ones(size(xi)) / (xi(2,1,1)-xi(1,1,1));
+pn = ones(size(xi)) / (yi(1,2,1)-yi(1,1,1));
+po = ones(size(xi)) / (zi(1,1,2)-zi(1,1,1));
+
+len = {0.1*mask,0.1*mask,0.*mask};
+
+fi = divand(mask,{pm,pn,po},{xi,yi,zi},{x,y,z},f,len,20,'primal',1,'alpha',[1 2 1]);
+len = {0.1*mask(:,:,1),0.1*mask(:,:,1)};
+
+fi2 = zeros(size(fi));
+for k = 1:size(mask,3);
+ fi2(:,:,k) = divand(mask(:,:,k),{pm(:,:,k),pn(:,:,k)},{xi(:,:,k),yi(:,:,k)},{x(:,k),y(:,k)},f(:,k),len,20);
+end
+
+rms = divand_rms(fi2,fi);
+
+if (rms > 1e-10)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',rms);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_4dvar.m b/inst/test_4dvar.m
new file mode 100644
index 0000000..a065cf6
--- /dev/null
+++ b/inst/test_4dvar.m
@@ -0,0 +1,53 @@
+% Testing divand in 4 dimensions.
+
+% grid of background field
+[xi,yi,zi,ti] = ndgrid(linspace(0,1,7));
+len = 1;
+k = pi/2;
+k = pi;
+fun = @(x,y,z,t) sin(k*x) .* sin(k*y) .* sin(k*z) .* sin(k*t);
+fi_ref = fun(xi,yi,zi,ti);
+
+% grid of observations
+[x,y,z,t] = ndgrid(linspace(2*eps,1-2*eps,5));
+x = x(:);
+y = y(:);
+z = z(:);
+t = t(:);
+
+pm = ones(size(xi)) / (xi(2,1,1,1)-xi(1));
+pn = ones(size(xi)) / (yi(1,2,1,1)-yi(1));
+po = ones(size(xi)) / (zi(1,1,2,1)-zi(1));
+pp = ones(size(xi)) / (ti(1,1,1,2)-ti(1));
+
+
+f = fun(x,y,z,t);
+mask = ones(size(xi));
+
+%f = x = y = z = t = 0.5
+
+fi = divand(mask,{pm,pn,po,pp},{xi,yi,zi,ti},{x,y,z,t},f,.12,10);
+rms = divand_rms(fi_ref,fi);
+
+% rms should be 0.0111807
+
+if (rms > 0.012)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',rms);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_divand.m b/inst/test_divand.m
new file mode 100644
index 0000000..b502ed1
--- /dev/null
+++ b/inst/test_divand.m
@@ -0,0 +1,55 @@
+% Test if divand is working correctly.
+
+tests = {...
+ 'test_interp_1d',...
+ 'test_interp_2d',...
+ 'test_interp_regular',...
+ 'test_sparse_diff',...
+ 'test_1dvar',...
+ 'test_2dvar',...
+ 'test_2dvar_check',...
+ 'test_2dvar_adv',...
+ 'test_2dvar_rellen',...
+ 'test_2dvar_lenxy',...
+ 'test_2dvar_check_correrr',...
+ 'test_2dvar_pcg',...
+ 'test_2dvar_constrain',...
+ 'test_2dvar_cyclic',...
+ 'test_2dvar_eof_check',...
+ 'test_3dvar',...
+ 'test_3dvar_large_stacked',...
+ 'test_4dvar'};
+
+% disable warning for full matrice product
+state = warning('query','divand:fullmat');
+warning('off','divand:fullmat')
+
+for iindex=1:length(tests);
+
+ fprintf('run %s: ',tests{iindex});
+ try
+ eval(tests{iindex});
+ colordisp(' OK ','green');
+ catch
+ colordisp(' FAIL ','red');
+ disp(lasterr)
+ end
+end
+
+% restore warning state
+warning(state.state,'divand:fullmat')
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_interp_1d.m b/inst/test_interp_1d.m
new file mode 100644
index 0000000..1a8bc6a
--- /dev/null
+++ b/inst/test_interp_1d.m
@@ -0,0 +1,37 @@
+% Testing 1D linear interpolation.
+
+[xi] = linspace(0,1,30);
+fi = sin( xi*6 );
+
+no = 20;
+x = rand(no,1);
+
+mask = ones(size(xi));
+
+[H,out] = interp_regular(mask,{xi},{x});
+
+f1 = interp1(xi,fi,x);
+
+f2 = H * fi(:);
+Difference = max(abs(f1(:) - f2(:)));
+
+if (Difference > 1e-6)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',Difference);
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_interp_2d.m b/inst/test_interp_2d.m
new file mode 100644
index 0000000..1ef1530
--- /dev/null
+++ b/inst/test_interp_2d.m
@@ -0,0 +1,39 @@
+% Testing 2D linear interpolation.
+
+[xi,yi] = ndgrid(linspace(0,1,30));
+fi = sin( xi*6 ) .* cos( yi*6);
+
+no = 20;
+x = rand(no,1);
+y = rand(no,1);
+
+mask = ones(size(xi));
+
+[H,out] = interp_regular(mask,{xi,yi},{x,y});
+
+f1 = interpn(xi,yi,fi,x,y);
+
+f2 = H * fi(:);
+Difference = max(abs(f1(:) - f2(:)));
+
+if (Difference > 1e-6)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',Difference);
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_interp_regular.m b/inst/test_interp_regular.m
new file mode 100644
index 0000000..7b62638
--- /dev/null
+++ b/inst/test_interp_regular.m
@@ -0,0 +1,83 @@
+% Testing linear interpolation on regular grid.
+
+D = [];
+
+[xi,yi,zi] = ndgrid(linspace(0,1,11));
+fi = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+
+no = 20;
+x = rand(no,1);
+y = rand(no,1);
+z = rand(no,1);
+
+[H,out] = interp_regular(ones(size(xi)),{xi,yi,zi},{x,y,z});
+f1 = interpn(xi,yi,zi,fi,x,y,z);
+
+f2 = H * fi(:);
+D(end+1) = max(abs(f1(:) - f2(:)));
+
+
+% 4d
+[xi,yi,zi,ti] = ndgrid(linspace(0,1,5));
+fi = sin(6*xi) .* cos(6*yi) .* sin(6*zi) .* cos(6*ti);
+no = 20;
+x = rand(no,1);
+y = rand(no,1);
+z = rand(no,1);
+t = rand(no,1);
+
+[H,out] = interp_regular(ones(size(xi)),{xi,yi,zi,ti},{x,y,z,t});
+f1 = interpn(xi,yi,zi,ti,fi,x,y,z,t);
+
+f2 = H * fi(:);
+D(end+1) = max(abs(f1(:) - f2(:)));
+
+% 2d with outside points
+
+[xi,yi] = ndgrid(linspace(0,1,30));
+fi = sin( xi*6 ) .* cos( yi*6);
+fi(1:4,:) = NaN;
+
+no = 200;
+x = 2*rand(no,1);
+y = 2*rand(no,1);
+
+
+%x = 0.805224734876956
+%y = 1.00822133283472
+
+%x = 0.134579335939626
+%y = 0.716156681034316
+
+f = sin( x*6 ) .* cos( y*6);
+
+mask = ~isnan(fi);
+
+[H,out] = interp_regular(mask,{xi,yi},{x,y});
+
+f1 = interpn(xi,yi,fi,x,y);
+f2 = H * fi(:);
+f2 = reshape(f2,size(x));
+
+D(end+1) = max(abs(f1(~out) - f2(~out)));
+
+
+if any(D > 1e-6) || any(isnan(f1(:) ~= out(:)))
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',Difference);
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
diff --git a/inst/test_sparse_diff.m b/inst/test_sparse_diff.m
new file mode 100644
index 0000000..73dca23
--- /dev/null
+++ b/inst/test_sparse_diff.m
@@ -0,0 +1,81 @@
+% Testing sparse operators.
+
+clear d
+f = randn(21,30,10);
+
+S = sparse_diff(size(f),1);
+f1 = f(2:end,:,:) - f(1:end-1,:,:);
+f2 = S*f(:);
+d(1) = max(abs(f1(:) - f2(:)));
+
+S = sparse_diff(size(f),2);
+f1 = f(:,2:end,:) - f(:,1:end-1,:);
+f2 = S*f(:);
+d(2) = max(abs(f1(:) - f2(:)));
+
+
+S = sparse_diff(size(f),3);
+f1 = f(:,:,2:end) - f(:,:,1:end-1);
+f2 = S*f(:);
+d(3) = max(abs(f1(:) - f2(:)));
+
+
+% cyclic
+
+% dim 1 cyclic
+fc = cat(1,f,f(1,:,:));
+dfc = fc(2:end,:,:) - fc(1:end-1,:,:);
+
+S = sparse_diff(size(f),1,true);
+f2 = S*f(:);
+d(4) = max(abs(dfc(:) - f2(:)));
+
+
+% dim 2 cyclic
+fc = cat(2,f,f(:,1,:));
+f1 = fc(:,2:end,:) - fc(:,1:end-1,:);
+
+S = sparse_diff(size(f),2,true);
+f2 = S*f(:);
+d(5) = max(abs(f1(:) - f2(:)));
+
+
+% stagger
+
+S = sparse_stagger(size(f),1);
+f1 = (f(2:end,:,:) + f(1:end-1,:,:))/2;
+f2 = S*f(:);
+d(end+1) = max(abs(f1(:) - f2(:)));
+
+
+% dim 1 cyclic
+fc = cat(1,f,f(1,:,:));
+f1 = (fc(2:end,:,:) + fc(1:end-1,:,:))/2;
+
+S = sparse_stagger(size(f),1,true);
+f2 = S*f(:);
+d(end+1) = max(abs(f1(:) - f2(:)));
+
+
+if any(d > 1e-6)
+ error('unexpected large difference with reference field');
+end
+
+fprintf('(max difference=%g) ',max(d));
+
+
+
+% Copyright (C) 2014 Alexander Barth <a.barth at ulg.ac.be>
+%
+% 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, see <http://www.gnu.org/licenses/>.
--
Alioth's /home/groups/pkg-octave/bin/git-commit-notice on /srv/git.debian.org/git/pkg-octave/octave-divand.git
More information about the Pkg-octave-commit
mailing list