Mercurial > forge
changeset 12377:a56cd865541e octave-forge
maint: moved divand package to separate individual hg repository
author | carandraug |
---|---|
date | Mon, 17 Feb 2014 16:38:14 +0000 |
parents | c0015cd07975 |
children | 355d7a8e6382 |
files | extra/divand/CITATION extra/divand/COPYING extra/divand/DESCRIPTION extra/divand/INDEX extra/divand/NEWS extra/divand/doc/README extra/divand/doc/gmd-7-225-2014.pdf extra/divand/inst/@CatBlockMat/CatBlockMat.m extra/divand/inst/@CatBlockMat/append.m extra/divand/inst/@CatBlockMat/ctranspose.m extra/divand/inst/@CatBlockMat/double.m extra/divand/inst/@CatBlockMat/end.m extra/divand/inst/@CatBlockMat/full.m extra/divand/inst/@CatBlockMat/mldivide.m extra/divand/inst/@CatBlockMat/mtimes.m extra/divand/inst/@CatBlockMat/size.m extra/divand/inst/@CatBlockMat/subsref.m extra/divand/inst/@CovarFun/CovarFun.m extra/divand/inst/@CovarFun/mldivide.m extra/divand/inst/@CovarIS/CovarIS.m extra/divand/inst/@CovarIS/diag.m extra/divand/inst/@CovarIS/factorize.m extra/divand/inst/@CovarIS/full.m extra/divand/inst/@CovarIS/inv.m extra/divand/inst/@CovarIS/mldivide.m extra/divand/inst/@CovarIS/mtimes.m extra/divand/inst/@CovarIS/subsref.m extra/divand/inst/@CovarLanczos/CovarLanczos.m extra/divand/inst/@CovarLanczos/diag.m extra/divand/inst/@CovarLanczos/mtimes.m extra/divand/inst/@CovarParam/CovarParam.m extra/divand/inst/@CovarParam/diag.m extra/divand/inst/@CovarParam/full.m extra/divand/inst/@CovarParam/isscalar.m extra/divand/inst/@CovarParam/isvector.m extra/divand/inst/@CovarParam/mldivide.m extra/divand/inst/@CovarParam/mtimes.m extra/divand/inst/@CovarParam/mtimescorr.m extra/divand/inst/@CovarParam/size.m extra/divand/inst/@CovarParam/subsref.m extra/divand/inst/@CovarSMW/CovarSMW.m extra/divand/inst/@CovarSMW/diag.m extra/divand/inst/@CovarSMW/full.m extra/divand/inst/@CovarSMW/mldivide.m extra/divand/inst/@CovarSMW/mtimes.m extra/divand/inst/@CovarSMW/size.m extra/divand/inst/@DiagBlockMat/DiagBlockMat.m extra/divand/inst/@DiagBlockMat/append.m extra/divand/inst/@DiagBlockMat/diag.m extra/divand/inst/@DiagBlockMat/full.m extra/divand/inst/@DiagBlockMat/inv.m extra/divand/inst/@DiagBlockMat/mldivide.m extra/divand/inst/@DiagBlockMat/mtimes.m extra/divand/inst/@DiagBlockMat/size.m extra/divand/inst/@MatFun/MatFun.m extra/divand/inst/@MatFun/ctranspose.m extra/divand/inst/@MatFun/double.m extra/divand/inst/@MatFun/end.m extra/divand/inst/@MatFun/full.m extra/divand/inst/@MatFun/mtimes.m extra/divand/inst/@MatFun/size.m extra/divand/inst/cat_cell_array.m extra/divand/inst/colordisp.m extra/divand/inst/colormsg.m extra/divand/inst/conjugategradient.m extra/divand/inst/divand.m extra/divand/inst/divand_addc.m extra/divand/inst/divand_background.m extra/divand/inst/divand_background_components.m extra/divand/inst/divand_constr_advec.m extra/divand/inst/divand_diagnose.m extra/divand/inst/divand_eof_contraint.m extra/divand/inst/divand_error.m extra/divand/inst/divand_factorize.m extra/divand/inst/divand_kernel.m extra/divand/inst/divand_laplacian.m extra/divand/inst/divand_metric.m extra/divand/inst/divand_obs.m extra/divand/inst/divand_operators.m extra/divand/inst/divand_orthogonalize_modes.m extra/divand/inst/divand_pc_michol.m extra/divand/inst/divand_pc_none.m extra/divand/inst/divand_pc_sqrtiB.m extra/divand/inst/divand_realistic_example.m extra/divand/inst/divand_rms.m extra/divand/inst/divand_simple_example.m extra/divand/inst/divand_solve.m extra/divand/inst/divand_solve_eof.m extra/divand/inst/interp_regular.m extra/divand/inst/localize_regular_grid.m extra/divand/inst/localize_separable_grid.m extra/divand/inst/mtimescorr.m extra/divand/inst/private/geomean.m extra/divand/inst/sparse_diag.m extra/divand/inst/sparse_diff.m extra/divand/inst/sparse_gradient.m extra/divand/inst/sparse_interp.m extra/divand/inst/sparse_pack.m extra/divand/inst/sparse_shift.m extra/divand/inst/sparse_stagger.m extra/divand/inst/sparse_trim.m extra/divand/inst/statevector_init.m extra/divand/inst/statevector_pack.m extra/divand/inst/statevector_unpack.m extra/divand/inst/test_1dvar.m extra/divand/inst/test_2dvar.m extra/divand/inst/test_2dvar_adv.m extra/divand/inst/test_2dvar_check.m extra/divand/inst/test_2dvar_check_correrr.m extra/divand/inst/test_2dvar_constrain.m extra/divand/inst/test_2dvar_cyclic.m extra/divand/inst/test_2dvar_eof_check.m extra/divand/inst/test_2dvar_lenxy.m extra/divand/inst/test_2dvar_pcg.m extra/divand/inst/test_2dvar_rellen.m extra/divand/inst/test_3dvar.m extra/divand/inst/test_3dvar_large_stacked.m extra/divand/inst/test_4dvar.m extra/divand/inst/test_divand.m extra/divand/inst/test_interp_1d.m extra/divand/inst/test_interp_2d.m extra/divand/inst/test_interp_regular.m extra/divand/inst/test_sparse_diff.m |
diffstat | 123 files changed, 0 insertions(+), 6911 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/divand/CITATION Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -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} - } -
--- a/extra/divand/COPYING Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,337 +0,0 @@ - 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.
--- a/extra/divand/DESCRIPTION Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -Name: divand -Version: 1.1.1 -Date: 2014-02-12 -Author: Alexander Barth <barth.alexander@gmail.com> -Maintainer: Alexander Barth <barth.alexander@gmail.com> -Title: divand -Description: divand performs an n-dimensional variational analysis (interpolation) of arbitrarily located observations. -Depends: octave (>= 3.4.0), mapping -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
--- a/extra/divand/INDEX Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,123 +0,0 @@ -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 - @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
--- a/extra/divand/NEWS Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -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
--- a/extra/divand/doc/README Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -/home/abarth/matlab/diva-nd/README \ No newline at end of file
--- a/extra/divand/inst/@CatBlockMat/CatBlockMat.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,49 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/append.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,35 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/ctranspose.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/double.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,29 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/end.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/size.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CatBlockMat/subsref.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarFun/CovarFun.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarFun/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,48 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/CovarIS.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,45 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/factorize.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/inv.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarIS/subsref.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarLanczos/CovarLanczos.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,28 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarLanczos/diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarLanczos/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/CovarParam.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/isscalar.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/isvector.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,39 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/mtimescorr.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/size.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarParam/subsref.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,35 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/CovarSMW.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,46 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@CovarSMW/size.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/DiagBlockMat.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,44 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/append.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,32 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/inv.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,28 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/mldivide.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@DiagBlockMat/size.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/MatFun.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,39 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/ctranspose.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/double.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/end.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/full.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/mtimes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/@MatFun/size.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/cat_cell_array.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,56 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/colordisp.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/colormsg.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,71 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/conjugategradient.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,198 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,290 +0,0 @@ -% 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@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
--- a/extra/divand/inst/divand_addc.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_background.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,265 +0,0 @@ -% 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 - - -coeff = divand_kernel(neff,alpha) * Ln; % units length^n - -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@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/>.
--- a/extra/divand/inst/divand_background_components.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,80 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_constr_advec.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_diagnose.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_eof_contraint.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,68 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_error.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,88 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_factorize.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,82 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_kernel.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,83 +0,0 @@ -% 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('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@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/>.
--- a/extra/divand/inst/divand_laplacian.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,119 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_metric.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,56 +0,0 @@ -% 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,:)); -dx = cat(1,dx(1,:),dx,dx(end,:)); - -dy = distance(lat(:,j-1),lon(:,j-1),lat(:,j+1),lon(:,j+1)); -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@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/>.
--- a/extra/divand/inst/divand_obs.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,83 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_operators.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,112 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_orthogonalize_modes.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,48 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_pc_michol.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,50 +0,0 @@ -% Compute a preconditioner using a modified incomplete Cholesky decomposition. -% -% [M1,M2] = diavnd_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] = diavnd_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@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/>.
--- a/extra/divand/inst/divand_pc_none.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% No preconditioner is used. -% -% [M1,M2] = diavnd_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] = diavnd_pc_none(iB,H,R) - -M1 = []; -M2 = []; - - - -% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- a/extra/divand/inst/divand_pc_sqrtiB.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% Compute a preconditioner using the Cholesky decomposition. -% -% [M1,M2] = diavnd_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, diavnd_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@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/>.
--- a/extra/divand/inst/divand_realistic_example.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,99 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_rms.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_simple_example.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,68 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_solve.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,87 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/divand_solve_eof.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/interp_regular.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,39 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/localize_regular_grid.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/localize_separable_grid.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,90 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/mtimescorr.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,46 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/private/geomean.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -% for compatibility with octave - -function m = geomean(x,varargin) - -m = prod(x,varargin{:})^(1/length(x)); -% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- a/extra/divand/inst/sparse_diag.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/sparse_diff.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -% 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@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
--- a/extra/divand/inst/sparse_gradient.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/sparse_interp.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,127 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/sparse_pack.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,35 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/sparse_shift.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/sparse_stagger.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/sparse_trim.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -% 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@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/>. -
--- a/extra/divand/inst/statevector_init.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -% 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@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@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
--- a/extra/divand/inst/statevector_pack.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -% 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@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@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
--- a/extra/divand/inst/statevector_unpack.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,64 +0,0 @@ -% 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@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@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
--- a/extra/divand/inst/test_1dvar.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -% 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(:); - -on = numel(x); -var = 0.01 * ones(on,1); -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@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/>.
--- a/extra/divand/inst/test_2dvar_adv.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,42 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_check.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,64 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_check_correrr.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_constrain.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,84 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_cyclic.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_eof_check.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_lenxy.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,54 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_pcg.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,65 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_2dvar_rellen.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_3dvar.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,54 +0,0 @@ -% 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(:); - -on = numel(x); -var = 0.01 * ones(on,1); -f = sin(6*x) .* cos(6*y) .* sin(6*z); - -if 0 - x = .5; - y = .5; - z = .5; - f = 1; - var = 0.01; -end - -m = 20; - -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)); - -fi = divand(mask,{pm,pn,po},{xi,yi,zi},{x,y,z},f,.1,100); -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@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/>.
--- a/extra/divand/inst/test_3dvar_large_stacked.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,71 +0,0 @@ -% 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); - - -on = numel(x); -var = 0.01 * ones(on,1); -f = sin(6*x) .* cos(6*y) .* sin(6*z); - -if 0 - x = .5; - y = .5; - z = .5; - f = 1; - var = 0.01; -end - -m = 20; - -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 = [.1 .1 .1]; -len = {.1,.1,.1}; - -len = zeros(size(mask)); -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@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/>.
--- a/extra/divand/inst/test_4dvar.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_divand.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_interp_1d.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_interp_2d.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,39 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_interp_regular.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,83 +0,0 @@ -% 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@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/>.
--- a/extra/divand/inst/test_sparse_diff.m Mon Feb 17 10:02:10 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -% 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@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/>.