Mercurial > forge
changeset 12365:626f1fdc3153 octave-forge
initial import to sourceforge
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/CITATION Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,22 @@ +Please cite divand as follows if you use this software in a publication: + + Barth, A., Beckers, J.-M., Troupin, C., Alvera-Azcárate, A., and + Vandenbulcke, L.: divand-1.0: n-dimensional variational data analysis + for ocean observations, Geosci. Model Dev., 7, 225-241, + doi:10.5194/gmd-7-225-2014, 2014. + URL http://www.geosci-model-dev.net/7/225/2014/ + +A BibTeX entry for LaTeX users is: + + @Article{gmd-7-225-2014, + AUTHOR = {Barth, A. and Beckers, J.-M. and Troupin, C. and Alvera-Azc\'arate, A. and Vandenbulcke, L.}, + TITLE = {divand-1.0: n-dimensional variational data analysis for ocean observations}, + JOURNAL = {Geoscientific Model Development}, + VOLUME = {7}, + YEAR = {2014}, + NUMBER = {1}, + PAGES = {225--241}, + URL = {http://www.geosci-model-dev.net/7/225/2014/}, + DOI = {10.5194/gmd-7-225-2014} + } +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/COPYING Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,337 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. <http://fsf.org/> + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + <one line to give the program's name and a brief idea of what it does.> + Copyright (C) <year> <name of author> + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + <signature of Ty Coon>, 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/DESCRIPTION Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,10 @@ +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 of arbitrarily located observations. +Depends: octave (>= 3.4.0) +License: GPLv2+ +Url: http://modb.oce.ulg.ac.be/mediawiki/index.php/divand, http://www.geosci-model-dev.net/7/225/2014/gmd-7-225-2014.html
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/INDEX Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,123 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/NEWS Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,15 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/doc/README Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,1 @@ +/home/abarth/matlab/diva-nd/README \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/CatBlockMat.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,49 @@ +% Create a matrix which represents the concatenation of a series of matrices. +% +% C = CatBlockMat(dim,X1,X2,...) +% +% Create a matrix C which represents the concatenation of a series of matrices +% X1,X2,... over dimension dim. The matrices X1,X2,... can be special matrix +% objects. The resulting matrix C behaves as a regular matrix. + +function retval = CatBlockMat(dim,varargin) + +self.d = 1:2; % only for matrices now +self.dim = dim; +self.B = varargin; +self.N = length(varargin); +self.i = zeros(self.N+1,1); + +self.sz = size(self.B{1}); +self.sz(dim) = 0; + + +for i=1:self.N + sz = size(self.B{i}); + + if any(self.sz(self.d ~= dim) ~= sz(self.d ~= dim)) + error('arguments dimensions are not consistent.'); + end + + self.i(i+1) = self.i(i) + sz(dim); +end + +self.sz(dim) = self.i(end); + +retval = class(self,'CatBlockMat'); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/append.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,35 @@ +% Append a matrix to a CatBlockMat matrix. +% +% A = append(B,M) +% +% Append the matrix M to a CatBlockMat matrix B. + + +function self = append(self,M) + +sz = size(M); + +if any(self.sz(self.d ~= self.dim) ~= sz(self.d ~= self.dim)) + error('arguments dimensions are not consistent.'); +end + +self.B{end+1} = M; +self.N = self.N+1; + + +self.i(end+1) = self.i(end) + sz(self.dim); +self.sz(self.dim) = self.i(end); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/ctranspose.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,33 @@ +% Conjugate transpose of a CatBlockMat matrix. +% +% T = ctranspose(C) +% +% Returns the conjugate transpose of matrix C. + + +function T = ctranspose(self) + +B = cell(1,self.N); + +for l=1:self.N + B{l} = self.B{l}'; +end + +d = 3-self.dim; +T = CatBlockMat(d,B{:}); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/double.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,29 @@ +% Convert a CatBlockMat matrix to a regular matrix. +% +% F = double(C) +% +% Convert the CatBlockMat C matrix to a regular matrix F. + +function F = double(self) + +b = {}; +for l=1:self.N + b{l} = double(self.B{l}); +end + +F = cat(self.dim,b{:}); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/end.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Last index of a CatBlockMat matrix. +% +% e = end(C,k) +% +% Return the last index of the CatBlockMat C along dimension k. + + +function e = end(self,k,n) + +e = self.sz(k); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a CatBlockMat matrix to a regular matrix. +% +% F = full(C) +% +% Convert the CatBlockMat C matrix to a regular matrix F. + +function F = full(self) + +F = cat(self.dim,self.B{:}); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,30 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of CatBlockMat A and another +% matrix B. + + +function p = mldivide(self,b) + +if self.sz(2) ~= size(b,1) + error('Inner matrix dimensions must agree.'); +end + +p = full(self) \ b; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,66 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the CatBlockMat A and another matrix B. + +function p = mtimes(A,B) + +if size(A,2) ~= size(B,1) + error('Inner matrix dimensions must agree.'); +end + +if isa(A,'CatBlockMat') + + if A.dim == 1 + p = zeros(size(A,1),size(B,2)); + + for l=1:A.N + i = A.i(l)+1:A.i(l+1); + p(i,:) = A.B{l} * B; + end + else + for l=1:A.N + i = (A.i(l)+1) :A.i(l+1); + + % must use subsref explicetly + % subsref(B,S) = B(i,:) + % http://www.mathworks.com/support/solutions/en/data/1-18J9B/ + + S.type = '()'; + S.subs = {i, ':'}; + tmp = A.B{l} * subsref(B,S); + + if l==1 + p = tmp; + else + p = p + tmp; + end + end + end + +elseif isa(B,'CatBlockMat') && B.dim == 2 + p = zeros(size(A,1),size(B,2)); + + for l=1:B.N + j = B.i(l)+1:B.i(l+1); + p(:,j) = A * B.B{l}; + end +else + warning('divand:fullmat','use full matrices for product'); + p = full(A) * full(B); +end +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/size.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Size of a CatBlockMat matrix. +% +% sz = size(C) +% sz = size(C,dim) +% +% Return the size of the CatBlockMat C matrix. If dim is specified, then only +% the size along dimension dim is returned. + +function sz = size(self,dim); + +sz = self.sz; + +if nargin == 2 + sz = sz(dim); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CatBlockMat/subsref.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,41 @@ +% Subreference of a CatBlockMat matrix. +% +% B = subsref(C,idx) +% +% Perform the subscripted element selection operation according to +% the subscript specified by idx. +% +% See also: subsref. + +function B = subsref(self,S) + +if strcmp(S.type,'()') + i = S.subs{self.dim}; + l = find(self.i+1 == i(1)); + + itest = self.i(l)+1:self.i(l+1); + + if any(i ~= itest) || ~all(strcmp(':',S.subs(self.d ~= self.dim))) + error('only whole blocks can be referenced'); + end + + B = self.B{l}; +else + error('Attempt to reference field of non-structure array'); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarFun/CovarFun.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,60 @@ +% Create a covariance matrix based on a function handler. +% +% CF = = CovarFun(n,fun) +% CF = = CovarFun(n,fun,...) +% +% Create the covariance matrix CF of size n-by-n based on the function handler +% fun. Optional key word parameters are: +% 'isvec': 1 if function is vectorized and 0 (default) if not. +% 'pc': function hander of the preconditioner for pcg +% 'M1','M2','matit','tol': parameters for pcg. +% +% See also: pcg + + +function retval = CovarFun(n,fun,varargin) + +isvec = 0; +self.pc = []; +self.M1 = []; +self.M2 = []; +self.tol = 1e-6; +self.maxit = 100; + +for i=1:2:length(varargin) + if strcmp(varargin{i},'isvec') + isvec = varargin{i+1}; + elseif strcmp(varargin{i},'pc') + self.pc = varargin{i+1}; + elseif strcmp(varargin{i},'M1') + self.M1 = varargin{i+1}; + elseif strcmp(varargin{i},'M2') + self.M2 = varargin{i+1}; + elseif strcmp(varargin{i},'maxit') + self.maxit = varargin{i+1}; + elseif strcmp(varargin{i},'tol') + self.tol = varargin{i+1}; + else + error('CovarFun:param','unknown param'); + end +end + +self.fun = fun; + +retval = class(self,'CovarFun',MatFun([n,n],fun,fun,isvec)); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarFun/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,48 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of CovarParam A and another +% matrix B using pcg. + +function p = mldivide(self,b) + +if size(self,2) ~= size(b,1) + error('Inner matrix dimensions must agree.'); +end + +if isempty(self.M1) + args = {self.pc}; +else + args = {self.M1,self.M2}; +end + +for i= 1:size(b,2) + + [p(:,i),flag,relres,iter] = pcg(self.fun, b(:,i), self.tol,... + self.maxit,args{:}); + + %whos b + self.pc + fprintf('iter %d %g \n',iter,relres); + if (flag ~= 0) + error('CovarFun:pcg', ['Preconditioned conjugate gradients method'... + ' did not converge %d %g %g'],flag,relres,iter); + end + +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/CovarIS.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Covariance matrix with a sparse inverse matrix. +% +% C = CovarIS(IS) +% +% Create the covariance matrix C whose inverse is the sparse matrix IS. +% To accelerate matrix product C*x, IS can be factorized. + +function retval = CovarIS(IS); + +self.IS = IS; +self.f = 0; +self.RP = []; +self.q = []; +self.PP = []; + +retval = class(self,'CovarIS'); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,45 @@ +% Diagonal elements. +% +% d = diag(C) +% +% Return the diagonal elements of the CovarIS matrix C. + +function d = diag(self) + +N = size(self.IS,1); +d = zeros(N,1); + + +if self.f + t0 = cputime; + + for i=1:N + % show progress every 5 seconds + if (t0 + 5 < cputime) + t0 = cputime; + fprintf('%d/%d\n',i,N); + end + + d(i) = self.PP(i,:) * (self.RP \ (self.RP' \ self.PP(i,:)')); + end +else + disp('factorize CovarIS'); + self_fact = factorize(self); + d = diag(self_fact); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/factorize.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,41 @@ +% Factorize matrix. +% +% self = factorize(self) +% +% Factorize matrix self using the Cholesky decomposition. + +function self = factorize(self) + + +if issparse(self.IS) + [self.RP, self.q, self.PP] = chol (self.IS); + + if self.q ~= 0 + error('factoziation failed (matrix is not positive definite)'); + end + + + %length(find(self.PP))/numel(self.PP) + %length(find(self.RP))/numel(self.RP) + +else + self.RP = chol(self.IS); + self.PP = speye(size(self.IS)); +end + +self.f = 1; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a CovarIS matrix to a regular matrix. +% +% F = full(C) +% +% Convert the CovarIS C matrix to a regular matrix F. + +function M = full(self) + +M = inv(self.IS); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/inv.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,23 @@ +% Inverse of a CovarIS matrix. +% +% M = inv(C) +% +% Return the inverse of the CovarIS matrix C. + +function M = inv(self) + +M = self.IS; +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of CovarIS A and another +% matrix B. + +function q = mldivide(self,b) + +q = self.IS * b; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,37 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the CovarIS matrix A and another matrix B. + +function p = mtimes(self,b) + +%keyboard +if self.f + %'fact' + p = self.PP * (self.RP \ (self.RP' \ (self.PP' * b))); +else +% 'direct' +% tic +% A = (self.IS + self.IS')/2; +% p = A \ b; +% toc +% tic + p = self.IS \ b; +% toc +% 'end' +end +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarIS/subsref.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,67 @@ +% Subreference of a CovarIS matrix. +% +% B = subsref(C,idx) +% +% Perform the subscripted element selection operation according to +% the subscript specified by idx. +% +% See also: subsref. + +function x = subsref(self,idx) + +assert(length(idx) == 1) +assert(strcmp(idx.type,'()')) +N = size(self.IS,1); + +i = idx.subs{1}; + +if strcmp(i,':') + i = 1:N; +end + +j = idx.subs{2}; + +if strcmp(j,':') + j = 1:N; +end + +x = zeros(length(i),length(j)); + +if self.f + t0 = cputime; + + k = 0; + for jj=1:length(j) + for ii=1:length(i) + + % show progress every 5 seconds + if (t0 + 5 < cputime) + t0 = cputime; + fprintf('%d/%d\n',k,numel(x)); + end + + x(ii,jj) = self.PP(i(ii),:) * (self.RP \ (self.RP' \ self.PP(j(jj),:)')); + k = k+1; + end + end +else + disp('factorize CovarIS'); + self_fact = factorize(self); + x = subsref(self_fact,idx); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarLanczos/CovarLanczos.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,28 @@ +% Covariance matrix reconstructed from Lanczos vectors. +% +% A = CovarLanczos(Q,T) +% +% Return a covariance matrix, factorized into A = Q * T * Q' +% where Q'*Q = I and T is a tridiagonal matrix. + +function retval = CovarLanczos(Q,T) + +self.Q = Q; +self.T = T; + +retval = class(self,'CovarLanczos'); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarLanczos/diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,30 @@ +% Diagonal elements. +% +% d = diag(C) +% +% Return the diagonal elements of the CovarLanczos matrix C. + +function d = diag(self) + +N = size(self.Q,1); +d = zeros(N,1); + +for i=1:N + d(i) = self.Q(i,:) * (self.T \ self.Q(i,:)'); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarLanczos/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the CovarLanczos matrix A and another matrix +% B. + +function p = mtimes(self,b) + +p = self.Q * (self.T \ (self.Q' * b)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/CovarParam.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,69 @@ +% Parametric covariance matrix with a guassian kernel. +% +% C = CovarParam({x,y,...},variance,len) +% +% Creates a parametric covariance matrix C with a guassian +% kernel, a correlation length of len and a variance of +% variance. The first parameter is a cell array with the coordinates. + +function retval = CovarParam(xi,variance,len,varargin) + +xi = squeeze(cat_cell_array(xi)); +m = size(xi,1); + + +% fun = @(xi1,xi2) exp(- sum( (xi1 - xi2).^2 , 2) / len^2); +% +% tic +% for j=1:m +% for i=1:m +% R2(i,j) = fun(xi(i,:),xi(j,:)); +% end +% end +% toc +% +% maxdiff(R,R2); + +self.tol = 1e-6; +self.maxit = 100; +self.var_add = 1; + +for i=1:2:length(varargin) + if strcmp(varargin{i},'tolerance') + self.tol = varargin{i+1}; + elseif strcmp(varargin{i},'maxit') + self.maxit = varargin{i+1}; + elseif strcmp(varargin{i},'fraccorr') + self.var_add = varargin{i+1}; + end +end + +self.m = m; +self.xi = xi; +self.len = len; + +if isscalar(variance) + variance = ones(self.m,1) * variance; +end + +self.variance = variance; +self.std = sqrt(variance); +self.S = sparse_diag(self.std); + +retval = class(self,'CovarParam'); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Diagonal elements. +% +% d = diag(C) +% +% Return the diagonal elements of the CovarParam matrix C. + +function d = diag(self); + +d = self.variance; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Convert a CovarParam matrix to a regular matrix. +% +% F = full(C) +% +% Convert the CovarParam C matrix to a regular matrix F. + + +function F = full(self) + +F = self * eye(size(self)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/isscalar.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Return true if argument is a scalar. +% +% b = isscalar(X) +% +% Return true if X is a scalar. + +function b = isscalar(self) + +b = isequal(size(self),[1 1]); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/isvector.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Return true if argument is a vector. +% +% b = isscalar(X) +% +% Return true if X is a vector. A vector is a N-by-1 or 1-by-N matrix. + +function b = isvector(self) + +sz = size(self); +b = length(sz) == 2 && (sz(1) == 1 || sz(2) == 1); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,47 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of CovarParam matrix A and +% another matrix B. + +function q = mldivide(self,b) + +q = zeros(size(b)); +b = self.S \ b; + +for i=1:size(b,2) + M = []; + %M = @(x) ifft(x); + + if 1 + [q(:,i),flag,relres,iter] = pcg(@(x) mtimescorr(self,x), b(:,i),self.tol,self.maxit,M); + else + + [q(:,i),flag,relres,iter] = pcg(@(x) fft(mtimescorr(self,ifft(x))), fft(b(:,i)),self.tol,self.maxit,M); + q(:,i) = real(ifft(q(:,i))); + end + + + if flag > 0 + relres,iter + warning('divand:pcgFlag','pcg flag %d',flag); + end +end + +q = self.S \ q; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,39 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the CovarParam matrix A and another matrix +% B. + +function p = mtimes(self,b) + +m = self.m; +p = zeros(size(b)); + +% for j=1:m +% for i=1:m +% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 ); +% p(i) = p(i) + exp(-d2 / self.len^2) * b(j); +% end +% end +% + +b = self.S*b; +p = mtimescorr(self,b); +p = self.S*p; + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/mtimescorr.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,40 @@ +% Matrix product with correlation matrix. +% +% p = mtimescorr(A,B) +% +% Return the matrix product between the CovarParam matrix A (with unit variance) +% and another matrix B. + + +function p = mtimescorr(self,b) + +m = self.m; +p = zeros(size(b)); + +% for j=1:m +% for i=1:m +% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 ); +% p(i) = p(i) + exp(-d2 / self.len^2) * b(j); +% end +% end +% + +p = mtimescorr(self.xi,self.len,b); +%p = mtimescorr_f(self.xi,self.len,b); +p = self.var_add * p + (1-self.var_add) * b; + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/size.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Size of a CovarParam matrix. +% +% sz = size(C) +% sz = size(C,dim) +% +% Return the size of the CovarParam C matrix. If dim is specified, then only +% the size along dimension dim is returned. + +function sz = size(self,dim); + +sz = [self.m, self.m]; + +if nargin == 2 + sz = sz(dim); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarParam/subsref.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,35 @@ +% Subreference of a CovarParam matrix. +% +% B = subsref(C,idx) +% +% Perform the subscripted element selection operation according to +% the subscript specified by idx. +% +% See also: subsref. + +function C = subsref(self,idx) + +assert(strcmp(idx.type,'()')) + +i = idx.subs{1}; +j = idx.subs{2}; + +assert(isequal(i,j)); + + +C = CovarParam(self.xi(i,:),self.variance(i),self.len); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/CovarSMW.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,46 @@ +% Covariance matrix than can be inverted using the Sherman-Morrison-Woodbury +% forumla. +% +% CSMW = CovarSMW(C,B) +% +% Create the covariance matrix CSMW = C + B*B' than can be inverted efficiently +% using the Sherman-Morrison-Woodbury formula, if size(B,2) is much smaller than +% size(C,1): +% +% inv(C + B*B') = inv(C) - inv(C)*B*inv(B'*inv(C)*B + I)*B'*inv(C) +% +% The symetric matrix C should implement the methods: size, diag, mtimes, +% mldivde and the matrix B should implement the methods: size, tranpose, mtimes, +% sum. + +function retval = CovarSMW(C,B) + +if size(C,1) ~= size(C,2) + error('C should be square matrix'); +end + +if size(C,1) ~= size(B,1) + error('size of C and B are not compatible'); +end + +self.C = C; +self.B = B; +self.D = inv(B' * (C \ B) + eye(size(B,2))); + +retval = class(self,'CovarSMW'); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Diagonal elements. +% +% d = diag(C) +% +% Return the diagonal elements of the CovarSWM matrix C. + +function d = diag(self); + +d = diag(self.C) + sum(self.B.^2,2); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a CovarSWM matrix to a regular matrix. +% +% F = full(C) +% +% Convert the CovarSWM C matrix to a regular matrix F. + +function F = full(self) + +F = self.C + self.B * self.B'; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of CovarSWM A and another +% matrix B. + +function Q = mldivide(self,A) + +Q = self.C \ A - self.C \ (self.B * (self.D * (self.B' * (self.C \ A)))); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the CovarSWM matrix A and another matrix B. + +function P = mtimes(self,A) + +P = self.C*A + self.B * (self.B'*A); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@CovarSMW/size.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Size of a CovarSWM matrix. +% +% sz = size(C) +% sz = size(C,dim) +% +% Return the size of the CovarSWM C matrix. If dim is specified, then only +% the size along dimension dim is returned. + +function sz = size(self,dim); + +sz = size(self.C); + +if nargin == 2 + sz = sz(dim); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/DiagBlockMat.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,44 @@ +% Block diagonal matrix. +% +% DB = DiagBlockMat(A1,A2,...,An) +% +% Create the block diagonal matrix DB from the matrices A1, A2,...,An +% +% / A1 0 0 ... 0 \ +% | 0 A2 0 : | +% DB = | 0 0 A3 | +% | : . 0 | +% \ 0 ... 0 An / + +function retval = DiagBlockMat(varargin) + +self.B = varargin; +self.N = length(varargin); +self.i = zeros(self.N+1,1); +self.j = zeros(self.N+1,1); + +for i=1:self.N + sz = size(self.B{i}); + self.i(i+1) = self.i(i) + sz(1); + self.j(i+1) = self.j(i) + sz(2); +end + +self.sz = [self.i(end), self.j(end)]; + +retval = class(self,'DiagBlockMat'); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/append.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,32 @@ +% Append a matrix to a DiagBlockMat matrix. +% +% A = append(B,M) +% +% Append the matrix M to a DiagBlockMat matrix B along its diagonal. + +function self = append(self,M) + +self.B{end+1} = M; +self.N = self.N+1; + +sz = size(M); + +self.i(end+1) = self.i(end) + sz(1); +self.j(end+1) = self.j(end) + sz(2); + +self.sz = [self.i(end), self.j(end)]; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,30 @@ +% Diagonal elements. +% +% d = diag(C) +% +% Return the diagonal elements of the matrix C. + +function d = diag(self) + +t = cell(1,self.N); + +for l=1:self.N + t{l} = diag(self.B{l}); +end + +d = cat(1,t{:}); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a DiagBlockMat matrix to a regular matrix. +% +% F = full(C) +% +% Convert the DiagBlockMat C matrix to a regular matrix F. + +function F = full(self) + +F = blkdiag(self.B{:}); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/inv.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,28 @@ +% Inverse of a DiagBlockMat matrix. +% +% M = inv(C) +% +% Return the inverse of the DiagBlockMat matrix C. + +function iA = inv(A) + +iA = A; + +for l=1:A.N + iA.B{l} = inv(A.B{l}); +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/mldivide.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,43 @@ +% Matrix left division A \ B. +% +% p = mldivide(A,B) +% +% Return the matrix product between the inverse of DiagBlockMat A and another +% matrix B + +function p = mldivide(self,b) + +if size(self,2) ~= size(b,1) + error('Inner matrix dimensions must agree.'); +end + + +if isa(b,'CatBlockMat') + B = {}; + for l=1:self.N + B{l} = self.B{l} \ b(self.j(l)+1:self.j(l+1),:); + end + + p = CatBlockMat(1,B{:}); +else + + p = zeros(size(self,1),size(b,2)); + + for l=1:self.N + p(self.i(l)+1:self.i(l+1),:) = self.B{l} \ b(self.j(l)+1:self.j(l+1),:); + end +end +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,33 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the DiagBlockMat matrix A and another matrix +% B. + +function p = mtimes(self,b) + +if size(self,2) ~= size(b,1) + error('Inner matrix dimensions must agree.'); +end + +p = zeros(size(self,1),size(b,2)); + +for l=1:self.N + p(self.i(l)+1:self.i(l+1),:) = self.B{l} * b(self.j(l)+1:self.j(l+1),:); +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@DiagBlockMat/size.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Size of a DiagBlockMat matrix. +% +% sz = size(C) +% sz = size(C,dim) +% +% Return the size of the DiagBlockMat C matrix. If dim is specified, then only +% the size along dimension dim is returned. + +function sz = size(self,dim); + +sz = self.sz; + +if nargin == 2 + sz = sz(dim); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/MatFun.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,39 @@ +% Matrix operator object based on a function handel. +% +% MF = MatFun(sz,fun,funt,isvec) +% +% Create a operator of size sz based on the linear function fun and its adjoint +% funt. For a matrix X of approriate size: +% +% MF*X is fun(X) and MF'*X is funt(X) +% +% If isvec is 0 (default), then the function must be applied on the columns of X +% individually, otherwise fun and funt are assumed to be "vectorized". + +function retval = MatFun(sz,fun,funt,isvec) + +if nargin == 3 + isvec = 0; +end + +self.sz = sz; +self.fun = fun; +self.funt = funt; +self.isvec = isvec; + +retval = class(self,'MatFun'); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/ctranspose.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,25 @@ +% Conjugate transpose of a MatFun matrix. +% +% T = ctranspose(C) +% +% Returns the conjugate transpose of matrix C. + +function T = ctranspose(self) + +T = MatFun(self.sz([2 1]),self.funt,self.fun); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/double.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a MatFun matrix to a regular matrix. +% +% F = full(C) +% +% Convert the MatFun C matrix to a regular matrix F. + +function d = double(self) + +d = full(self); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/end.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,23 @@ +% Last index of a MatFun matrix. +% +% e = end(C,k) +% +% Return the last index of the MatFun matrix C along dimension k. + +function e = end(self,k,n) + +e = self.sz(k); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/full.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,24 @@ +% Convert a MatFun matrix to a regular matrix. +% +% F = full(C) +% +% Convert the MatFun C matrix to a regular matrix F. + +function M = full(self); + +M = self * eye(self.sz(2)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/mtimes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,40 @@ +% Matrix product A*B. +% +% p = mtimes(A,B) +% +% Return the matrix product between the MatFun matrix A and another matrix B. + +function p = mtimes(A,B) + +if size(A,2) ~= size(B,1) + error('Inner matrix dimensions must agree.'); +end + +if isa(A,'MatFun') + if A.isvec + p = A.fun(B); + else + p = zeros(size(A,1),size(B,2)); + + for l=1:size(B,2) + p(:,l) = A.fun(B(:,l)); + end + end +else + error('not implemented'); +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/@MatFun/size.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% Size of a MatFun matrix. +% +% sz = size(C) +% sz = size(C,dim) +% +% Return the size of the MatFun matrix C. If dim is specified, then only +% the size along dimension dim is returned. + +function sz = size(self,dim); + +sz = self.sz; + +if nargin == 2 + sz = sz(dim); +end + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/cat_cell_array.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,56 @@ +% Concatenate a cell array. +% +% x = cat_cell_array(c) +% +% Concatenate all arrays of the cell array c along the dimension length(c)+1 + + +function x = cat_cell_array(x) + +do_grid = 0; + +if iscell(x) + n = length(x); + + if do_grid + all_vectors = 1; + + for i=1:n + all_vectors = all_vectors & isvector(x{i}); + end + + if all_vectors + [x{:}] = ndgrid(x{:}); + end + end + + sz = size(x{1}); + for i=2:n + if ~isequal(sz,size(x{i})) + %sz + %size(x{i}) + error('all arguments in cell array should have the same size'); + end + end + + x = cat(n+1,x{:}); + end +end + + + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/colordisp.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,27 @@ +% Display a message in a color (followed by a newline). +% +% colordisp(msg,color) +% +% The same as disp but in color. +% +% See also +% colormsg + +function colordisp(msg,color) +colormsg (msg,color) +fprintf(1,'\n'); + +% Copyright (C) 2004 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/colormsg.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,71 @@ +% Display a message in a color (without a newline). +% +% colormsg(msg,color) +% +% The message msg in printed on the terminal the specified color. Acceptable +% values color color are: +% black +% red +% green +% yellow +% blue +% magenta +% cyan +% white +% +% In Matlab the message is displayed in black (due to the lack of ANSI escape +% code support). +% +% See also +% colordisp + +function colormsg (msg,color) + +%getenv('TERM') +%if strcmp(getenv('TERM'),'xterm') % && exist('puts') > 1 +% only in octave +if exist('puts') > 1 + esc = char(27); + + % ANSI escape codes + colors.black = [esc, '[40m']; + colors.red = [esc, '[41m']; + colors.green = [esc, '[42m']; + colors.yellow = [esc, '[43m']; + colors.blue = [esc, '[44m']; + colors.magenta = [esc, '[45m']; + colors.cyan = [esc, '[46m']; + colors.white = [esc, '[47m']; + + reset = [esc, '[0m']; + + c = getfield(colors,color); + + %oldpso = page_screen_output (0); + + try + fprintf([c, msg, reset]); + %puts (reset); % reset terminal + catch + %page_screen_output (oldpso); + end +else + fprintf(msg); +end + + + +% Copyright (C) 2004 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/conjugategradient.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,198 @@ +% Solve a linear system with the preconditioned conjugated-gradient method. +% +% [x] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0); +% [x,Q,T] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0); +% +% solves for x +% A x = b +% using the preconditioned conjugated-gradient method. +% It provides also an approximation of A: +% A \sim Q*T*Q' +% +% J(x) = 1/2 (x' b - x' A x) +% ∇ J = b - A x +% A x = b - ∇ J +% b = ∇ J(0) +% +% the columns of Q are the Lanczos vectors + +% Alexander Barth 2010,2012-2013 + +function [x,Q,T,diag] = conjugategradient(fun,b,varargin) + +n = length(b); + +% default parameters +maxit = min(n,20); +minit = 0; +tol = 1e-6; +pc = []; +x0 = []; +renorm = 0; + +prop = varargin; + +for i=1:2:length(prop) + if strcmp(prop{i},'minit') + minit = prop{i+1}; + elseif strcmp(prop{i},'maxit') + maxit = prop{i+1}; + elseif strcmp(prop{i},'tol') + tol = prop{i+1}; + elseif strcmp(prop{i},'x0') + x0 = prop{i+1}; + elseif strcmp(prop{i},'pc') + pc = prop{i+1}; + elseif strcmp(prop{i},'renorm') + renorm = prop{i+1}; + else + + end +end + + +if isempty(x0) + % random initial vector + x0 = randn(n,1); +end + + +if isempty(pc) + % no pre-conditioning + pc = @(x) x; +elseif isnumeric(pc) + invM = pc; + pc = @(x)( invM *x); +end + + +tol2 = tol^2; + +delta = []; +gamma = []; +q = NaN*zeros(n,1); + +%M = inv(invM); +%E = chol(M); + + +% initial guess +x = x0; + +% gradient at initial guess +r = b - fun(x); + +% quick exit +if r'*r < tol2 + return +end + +% apply preconditioner +z = pc(r); + +% first search direction == gradient +p = z; + +% compute: r' * inv(M) * z (we will need this product at several +% occasions) + +zr_old = r'*z; + +% r_old: residual at previous iteration +r_old = r; + +for k=1:maxit + if k <= n && nargout > 1 + % keep at most n vectors + Q(:,k) = r/sqrt(zr_old); + end + + % compute A*p + Ap = fun(p); + %maxdiff(A*p,Ap) + + % how far do we need to go in direction p? + % alpha is determined by linesearch + + % alpha z'*r / (p' * A * p) + alpha(k) = zr_old / ( p' * Ap); + + % get new estimate of x + x = x + alpha(k)*p; + + % recompute gradient at new x. Could be done by + % r = b-fun(x); + % but this does require an new call to fun + r = r - alpha(k)*Ap; + + if renorm + r = r - Q(:,1:k) * Q(:,1:k)' * r ; + end + %maxdiff(r,b-fun(x)) + + % apply pre-conditionner + z = pc(r); + + + zr_new = r'*z; + + if r'*r < tol2 && k >= minit + break + end + + %Fletcher-Reeves + beta(k+1) = zr_new / zr_old; + %Polak-Ribiere + %beta(k+1) = r'*(r-r_old) / zr_old; + %Hestenes-Stiefel + %beta(k+1) = r'*(r-r_old) / (p'*(r-r_old)); + %beta(k+1) = r'*(r-r_old) / (r_old'*r_old); + + + % norm(p) + p = z + beta(k+1)*p; + zr_old = zr_new; + r_old = r; +end + + +%disp('alpha and beta') +%figure,plot(beta(2:end)) +%rg(alpha) +%rg(beta(2:end)) + +if nargout > 1 + kmax = size(Q,2); + + delta(1) = 1/alpha(1); + + %delta(1) - Q(:,1)'*invM*A*invM*Q(:,1) + + + for k=1:kmax-1 + delta(k+1) = 1/alpha(k+1) + beta(k+1)/alpha(k); + gamma(k) = -sqrt(beta(k+1))/alpha(k); + end + + T = sparse([1:kmax 1:kmax-1 2:kmax ],... + [1:kmax 2:kmax 1:kmax-1],... + [delta gamma gamma]); + + diag.iter = k; + diag.relres = sqrt(r'*r); +end + +% Copyright (C) 2004 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,290 @@ +% Compute a variational analysis of arbitrarily located observations. +% +% [fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...); +% [fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...); +% [fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...); +% +% Perform an n-dimensional variational analysis of the observations f located at +% the coordinates x. The array fi represent the interpolated field at the grid +% defined by the coordinates xi and the scales factors pmn. +% +% Input: +% mask: binary mask delimiting the domain. 1 is inside and 0 outside. +% For oceanographic application, this is the land-sea mask. +% +% pmn: scale factor of the grid. pmn is a cell array with n elements. Every +% element represents the scale factor of the corresponding dimension. Its +% inverse is the local resolution of the grid in a particular dimension. +% +% xi: cell array with n elements. Every element represents a coordinate +% of the final grid on which the observations are interpolated +% x: cell array with n elements. Every element represents a coordinate of +% the observations +% f: value of the observations *minus* the background estimate (m-by-1 array). +% (see note) +% len: correlation length +% lambda: signal-to-noise ratio of observations (if lambda is a scalar). +% The larger this value is, the closer is the field fi to the +% observation. +% if lambda is a scalar: +% R = 1/lambda I, where R is the observation error covariance +% matrix), +% if lambda is a vector +% a vector (R = diag(lambda)) or a matrix or a CovarParam object (R = +% lambda). +% +% +% Optional input arguments specified as pairs of keyword and values: +% 'velocity', vel: velocity of advection constraint. The default is +% no-advection constraint +% +% 'alpha': alpha is vector of coefficients multiplying various terms in the +% cost function. The first element multiplies the norm. +% The other i-th element of alpha multiplies the (i+1)-th derivative. +% Per default, the highest derivative is m = ceil(1+n/2) where n is the +% dimension of the problem. +% +% The values of alpha is the (m+1)th row of the Pascal triangle: +% m=0 1 +% m=1 1 1 +% m=1 1 2 1 (n=1,2) +% m=2 1 3 3 1 (n=3,4) +% ... +% +% 'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or +% off (0, default). If on, they will be returned as the last output +% argument +% +% 'EOF', EOF: sub-space constraint. Orthogonal (EOF' WE^2 EOF = I) (units of +% EOF: m^(-n/2)) +% +% 'EOF_scaling', EOF_scaling: (dimensional) +% +% 'constraint': a structure with user specified constrain +% +% 'moddim': modulo for cyclic dimension (vector with n elements). +% Zero is used for non-cyclic dimensions. Halo points should +% not be included for cyclic dimensions. For example if the first dimension +% is cyclic, then the grid point corresponding to mask(1,j) should be +% between mask(end,1) (left neighbor) and mask(2,j) (right neighbor) +% +% 'fracdim': fractional indices (n-by-m array). If this array is specified, +% then x and xi are not used. +% +% 'inversion': direct solver ('chol' for Cholesky factorization) or a +% interative solver ('pcg' for preconditioned conjugate gradient) can be +% used. +% +% 'compPC': function that returns a preconditioner for the primal formulation +% if inversion is set to 'pcg'. The function has the following arguments: +% +% [M1,M2] = compPC(iB,H,R) +% +% where iB is the inverse background error covariance, H the observation +% operator and R the error covariance of the observation. The used +% preconditioner M will be M = M1 * M2 (or just M = M1 if M2 is empty). +% Per default a modified incomplete Cholesky factorization will be used a +% preconditioner. +% +% Note: 'velocity' and 'constraint' may appear multiple times +% +% Output: +% fi: the analysed field +% err: error variance of the analysis field relative to the error variance of +% the background +% s: structure +% s.iB: adimensional +% s.E: scaled EOF (dimensional) +% +% Note: +% If zero is not a valid first guess for your variable (as it is the case for +% e.g. ocean temperature), you have to subtract the first guess from the +% observations before calling divand and then add the first guess back in. +% +% Example: +% see divand_simple_example.m +% + + +function varargout = divand(mask,pmn,xi,x,f,len,lambda,varargin) + +% default values + +velocity = {}; +EOF = []; +diagnostics = 0; +EOF_lambda = 0; +primal = 1; +factorize = 1; +tol = 1e-6; +maxit = 100; +minit = 10; +constraints = {}; +inversion = 'chol'; +moddim = []; +fracindex = []; +alpha = []; +keepLanczosVectors = 0; +compPC = @divand_pc_sqrtiB; + +prop = varargin; +for i=1:length(prop)/2 + if strcmp(prop{2*i-1},'velocity') + velocity{end+1} = prop{2*i}; + elseif strcmpi(prop{2*i-1},'EOF') + EOF = prop{2*i}; + elseif strcmpi(prop{2*i-1},'EOF_scaling') + EOF_lambda = prop{2*i}; + elseif strcmpi(prop{2*i-1},'primal') + primal = prop{2*i}; + elseif strcmpi(prop{2*i-1},'factorize') + factorize = prop{2*i}; + elseif strcmpi(prop{2*i-1},'inversion') + inversion = prop{2*i}; + elseif strcmpi(prop{2*i-1},'tol') + tol = prop{2*i}; + elseif strcmpi(prop{2*i-1},'maxit') + maxit = prop{2*i}; + elseif strcmpi(prop{2*i-1},'minit') + minit = prop{2*i}; + elseif strcmpi(prop{2*i-1},'diagnostics') + diagnostics = prop{2*i}; + elseif strcmpi(prop{2*i-1},'constraint') + constraints{end+1} = prop{2*i}; + elseif strcmpi(prop{2*i-1},'fracindex') + fracindex = prop{2*i}; + elseif strcmpi(prop{2*i-1},'moddim') + moddim = prop{2*i}; + elseif strcmpi(prop{2*i-1},'alpha') + alpha = prop{2*i}; + elseif strcmpi(prop{2*i-1},'keepLanczosVectors') + keepLanczosVectors = prop{2*i}; + elseif strcmpi(prop{2*i-1},'compPC') + compPC = prop{2*i}; + else + warning('divand:unknownProperty','unknown property "%s" is irgnored',prop{2*i-1}); + end +end + +% check inputs + +if ~any(mask(:)) + error('no sea points in mask'); +end + +s = divand_background(mask,pmn,len,alpha,moddim); +s.betap = 0; +s.EOF_lambda = EOF_lambda; +s.primal = primal; +s.factorize = factorize; +s.tol = tol; +s.maxit = maxit; +s.minit = minit; +s.inversion = inversion; +s.keepLanczosVectors = keepLanczosVectors; +s.compPC = compPC; + +% remove non-finite elements from observations +f = f(:); +valid = isfinite(f); +x = cat_cell_array(x); + +if ~all(valid) + fprintf(1,'remove %d (out of %d) non-finite elements from observation vector\n',sum(~valid),numel(f)); + x = reshape(x,[length(f) s.n]); + f = f(valid); + x = reshape(x(repmat(valid,[1 s.n])),[length(f) s.n]); + + if ~isempty(fracindex) + fracindex = fracindex(:,valid); + end + + if isscalar(lambda) + % do nothing + elseif isvector(lambda) + lambda = lambda(valid); + elseif ismatrix(lambda) + lambda = lambda(valid,valid); + end +end + +apply_EOF_contraint = ~(isempty(EOF) | all(EOF_lambda == 0)); + +s.mode = 1; + +if ~apply_EOF_contraint + s.betap = 0; +else + if s.mode==0 + s.betap = max(EOF_lambda)/s.coeff; % units m^(-n) + elseif s.mode==1 + s.betap = max(max(EOF_lambda)-1,0)/s.coeff; + end +end + +%assert(s.betap,0,1e-8) +% increase contraint on total enegery to ensure system is still positive defined +%s.betap +s.iB = s.iB + s.betap * s.WE'*s.WE; + +% add observation constrain to cost function +s = divand_addc(s,divand_obs(s,xi,x,f,lambda,fracindex)); + +% add advection constraint to cost function +for i=1:length(velocity) + %s = divand_advection(s,velocity); + s = divand_addc(s,divand_constr_advec(s,velocity{i})); +end + + +% add all additional constrains +for i=1:length(constraints) + s = divand_addc(s,constraints{i}); +end + + +if apply_EOF_contraint + s = divand_eof_contraint(s,EOF_lambda,EOF); +end + +% factorize a posteori error covariance matrix +% or compute preconditioner +s = divand_factorize(s); + +if ~apply_EOF_contraint + [fi,s] = divand_solve(s,f); +else + [fi,s] = divand_solve_eof(s,f); +end + +varargout{1} = fi; + +if nargout-diagnostics >= 2 + err = divand_error(s); + varargout{2} = err; +end + +if diagnostics + s.B = CovarIS(s.iB); + [s.Jb,s.Jo,s.Jeof,s.J] = divand_diagnose(s,fi,f); + s.valid = valid; + varargout{nargout} = s; +end + + +% Copyright (C) 2008-2014 Alexander Barth <barth.alexander@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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_addc.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,41 @@ +% Add a constraint to the cost function. +% +% s = divand_addc(s,constrain) +% +% Include in the structure s the specified constrain. +% +% Input: +% s: structure created by divand_background +% constrain: The parameter constrain has the following fields: R (a covariance +% matrix), H (extraction operator) and yo (specified value for the +% constrain). +% +% Output: +% s: structure to be used by divand_factorize + +function s = divand_addc(s,constrain) + +if isfield(s,'R') + s.R = append(s.R,constrain.R); + s.H = append(s.H,constrain.H); + s.yo = cat(1,s.yo,constrain.yo); +else + s.H = CatBlockMat(1,constrain.H); + s.R = DiagBlockMat(constrain.R); + s.yo = constrain.yo; +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_background.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,265 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_background_components.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,80 @@ +% Form the different components of the background error covariance matrix. +% +% [iB_,iB] = divand_background_components(s,alpha) +% +% Compute the components of the background error covariance matrix iB_ and +% their sum based on alpha (the a-dimensional coefficients for norm, gradient, +% laplacian,...). + +function [iB_,iB] = divand_background_components(s,alpha) + +WE = s.WE; +D = s.D; +coeff = s.coeff; +n = s.n; + +iB_ = cell(length(alpha),1); + +% constrain of total norm + +iB_{1} = (1/coeff) * (WE'*WE); + + +% loop over all derivatives + +for j=2:length(alpha) + % exponent of laplacian + k = floor((j-2)/2); + + if mod(j,2) == 0 + % constrain of derivative with uneven order (j-1) + % (gradient, gradient*laplacian,...) + % normalized by surface + + iB_{j} = sparse(size(D,1),size(D,1)); + + for i=1:n + Dx = s.WEss{i} * s.Dx{i} * D^k; + %Dx = s.WEs{i} * s.Dxs{i} * D^k; + iB_{j} = iB_{j} + Dx'*Dx; + end + else + % constrain of derivative with even order (j-1) + % (laplacian, biharmonic,...) + + % normalize by surface of each cell + % such that inner produces (i.e. WE'*WE) + % become integrals + % WD: units length^(n/2) + + WD = WE * D^(k+1); + iB_{j} = WD'*WD; + end + + iB_{j} = iB_{j}/coeff; +end + + +% sum all terms of iB +% iB is adimentional +iB = alpha(1) * iB_{1}; +for j=2:length(alpha) + iB = iB + alpha(j) * iB_{j}; +end + +% LocalWords: iB divand + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_constr_advec.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,67 @@ +% Create the advection constrain. +% +% c = divand_constr_advec(s,velocity) +% +% Create the advection constrain using the specified velocity. +% +% Input: +% s: structure created by divand_background +% velocity: cell array of velocity vectors +% +% Output: +% c: structure to be used by divand_addc with the following fields: R (a +% covariance matrix), H (extraction operator) and yo (specified value for +% the constrain). + +function c = divand_constr_advec(s,velocity) + +velocity = cat_cell_array(velocity); + +% check for NaNs +nancount = sum(isnan(velocity(:))); +if nancount > 0 + error('%d velocity values are equal to NaN',nancount); +end + +mask = s.mask; +n = s.n; +iscyclic = s.iscyclic; + +sz = size(mask); + +A = sparse(0); + +vel = reshape(velocity,numel(mask),n); + + + +for i=1:n + S = sparse_stagger(sz,i,iscyclic(i)); + m = (S * mask(:) == 1); + + d = vel(:,i); + + A = A + sparse_diag(d(mask==1)) * sparse_pack(mask) * S' * sparse_pack(m)' * s.Dx{i}; +end + +l = size(A,1); + +c.H = A; +c.yo = sparse(l,1); +c.R = speye(l); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_diagnose.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,47 @@ +% Computes diagnostics for the analysis. +% +% [Jb,Jo,Jeof,J] = divand_diagnose(s,fi,d) +% +% Computes the value of different terms of the cost-function +% +% Input: +% s: structure created by divand_solve +% fi: analysis +% d: observations +% +% Output: +% Jb: background constrain +% Jeof: observation constrain +% Jeof: EOF constrain + +function [Jb,Jo,Jeof,J] = divand_diagnose(s,fi,d) + +d = s.yo; +xa = statevector_pack(s.sv,fi); + +Jb = xa' * s.iB * xa - s.betap * (s.WE'*xa)' * (s.WE'*xa); +Jo = (s.H*xa - d)' * (s.R \ (s.H*xa - d)); + +if all(s.EOF_lambda == 0) + Jeof = 0; +else + Jeof = s.betap * ((s.WE'*xa)' * (s.WE'*xa)) - (s.E'*xa)' * (s.E'*xa); +end + +Jo = full(Jo); + +J = Jb + Jo + Jeof; +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_eof_contraint.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,68 @@ +% Include constraint from EOFs. +% +% s = divand_eof_contraint(s,EOF_lambda,EOF) +% +% Include the constraint from the EOF and their eigenvalues (EOF_lambda) +% in the cost function described by the structure s. +% + +function s = divand_eof_contraint(s,EOF_lambda,EOF) + +iB = s.iB; +sv = s.sv; +WE = s.WE; + +coeff = s.coeff; + +% remove land points +E = statevector_pack(sv,EOF); +% units m^(-n/2) + +% normalize by surface +% for coeff see divand_background.m +%E = WE * E; +%EOF_lambda + +E = WE^2 * E * diag(sqrt(EOF_lambda / coeff)); +%% units: m^(n) m^(-n/2) sqrt(m^(-n)) = 1 + +if 0 + +BE = iB \ E; + + + +A = BE * inv(E'*BE - eye(size(E,2))); +Beof_var = - sum(A .* BE,2); + +% scaling factor to enshure that: +% scaling * inv(iB - E E') has a variance of 1 far from boundaries + +scaling = 1/(1 + mean(Beof_var)) +scaling = 1 + +%disp('apply scaling') +E = 1/sqrt(scaling) * E; +iB = iB/scaling; + +s.scaling = scaling; + +end + +s.E = E; +s.iB = iB; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_error.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,88 @@ +% Compute the expected error variance of the analysis. +% +% err = divand_error(s) +% +% Compute the expected error variance of the analysis based on the +% background error covariance, observation error covariance and data +% distribution. +% +% Input: +% s: structure created by divand_factorize +% +% Output: +% err: expected error variance of the analysis + + +function err = divand_error(s) + +H = s.H; +sv = s.sv; +N = sum(s.mask(:)==1); +n = s.n; + +errp = zeros(N,1); + +if s.primal + errp = diag(s.P); + + %global P_CL + %errp = diag(P_CL); +else + %dual + B = s.B; + R = s.R; + + % fun(x) computes (H B H' + R)*x + fun = @(x) H * (B * (H'*x)) + R*x; + t0 = cputime; + + for i=1:N + % show progress every 5 seconds + if (t0 + 5 < cputime) + t0 = cputime; + fprintf('%d/%d\n',i,N); + end + + e = zeros(N,1); + e(i) = 1; + % analyzed covariance + % P = B - B * H' * inv(H * B * H' + R) * H * B + % analyzed variance + % P = e' B e - e' * B * H' * inv(H * B * H' + R) * H * B * e + + Be = B * e; + errp(i) = e' * Be; + HBe = H*Be; + %keyboard + + % tmp is inv(H * B * H' + R) * H * B * e + %tic + [tmp,flag,relres,iter] = pcg(fun, HBe,s.tol,s.maxit,s.funPC); + %toc + + if (flag ~= 0) + error('divand:pcg', ['Preconditioned conjugate gradients method'... + ' did not converge %d %g %g'],flag,relres,iter); + end + + errp(i) = errp(i) - HBe' * tmp; + end +end + +err = statevector_unpack(sv,errp); +err(~s.mask) = NaN; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_factorize.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,82 @@ +% Factorize some matrices to increase performance. +% +% s = divand_factorize(s) +% +% Create the inverse of the posteriori covariance matrix +% and factorize +% +% Input: +% s: structure created by divand_background +% +% Output: +% s.P: factorized P + + +function s = divand_factorize(s) + +R = s.R; +iB = s.iB; +H = s.H; + +if s.primal + if strcmp(s.inversion,'chol') + iP = iB + H'*(R\H); + P = CovarIS(iP); + + % Cholesky factor of the inverse of a posteriori + % error covariance iP + if s.factorize + P = factorize(P); + end + + s.P = P; + else + tic + [s.M1,s.M2] = s.compPC(iB,H,R); + s.pc_time = toc(); + end +else % dual + %C = H * (iB \ H') + R; + s.B = CovarIS(iB); + % pre-conditioning function for conjugate gradient + s.funPC = []; + + if s.factorize + s.B = factorize(s.B); + + % iM = inv(H * B * H' + sparse_diag(diag(R))) + % iM * (H * B * H' + R) is identify matrix if R is diagonal + + M = H * (s.B * H') + sparse_diag(diag(R)); + %M = H * (s.B * H') + sparse_diag(sum(R,1)); + iM = CovarIS(M); + iM = factorize(iM); + + % pre-conditioning function for conjugate gradient + s.funPC = @(x) iM*x; + + if 0 + C = H * (s.B * H') + sparse_diag(diag(R)); + [RC, q, QC] = chol (C); + s.funPC = @(x) RC' \ (QC*x); + + end + end +end + +% LocalWords: divand + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_kernel.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,83 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_laplacian.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,119 @@ +% Create the laplacian operator. +% +% Lap = divand_laplacian(mask,pmn,nu,iscyclic) +% +% Form a Laplacian using finite differences +% assumes that gradient is zero at "coastline" +% +% Input: +% mask: binary mask delimiting the domain. 1 is inside and 0 outside. +% For oceanographic application, this is the land-sea mask. +% pmn: scale factor of the grid. +% nu: diffusion coefficient of the Laplacian +% field of the size mask or cell arrays of fields +% +% Output: +% Lap: sparce matrix represeting a Laplaciant +% +% +function Lap = divand_laplacian(mask,pmn,nu,iscyclic) + +pmn = cat_cell_array(pmn); + +% number of dimensions +n = size(pmn,ndims(pmn)); + +sz = size(mask); + +if nargin == 2 + nu = ones(sz); +end + +if isequal(sz,size(nu)) + % assume diffusion is the same along all dimensions + nu = repmat(nu,[ones(1,n) n]); +end + +% nu: 1st index represents the dimensions +nu = permute(nu,[n+1 1:n]); + +% extraction operator of sea points +H = sparse_pack(mask==1); + + +sz = size(mask); + +%DD = sparse(0); +DD = sparse(prod(sz),prod(sz)); + +Pmn = reshape(pmn,[prod(sz) n]); +for i=1:n + % operator for staggering in dimension i + S = sparse_stagger(sz,i,iscyclic(i)); + + % d = 1 for interfaces surounded by water, zero otherwise + d = S * mask(:) == 1; + + % metric + for j = 1:n + tmp = S * Pmn(:,j); + + if j==i + d = d .* tmp; + else + d = d ./ tmp; + end + end + + % tmp: "diffusion coefficient" along dimension i + tmp = nu(i,:); + + %whos tmp d S nu + %keyboard + d = d .* (S * tmp(:)); + + % Flux operators D + % zero at "coastline" + + D = sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)); + szt = sz; + + if ~iscyclic(i) + % extx: extension operator + szt(i) = szt(i)+1; + extx = sparse_trim(szt,i)'; + + D = extx * D; + else + % shift back one grid cell along dimension i + D = sparse_shift(sz,i,iscyclic(i))' * D; + end + + % add laplacian along dimension i + DD = DD + sparse_diff(szt,i,iscyclic(i)) * D; +end + +ivol = prod(pmn,n+1); + +% Laplacian on regualar grid DD +DD = sparse_diag(ivol) * DD; + +% Laplacian on grid with on sea points Lap +Lap = H * DD * H'; + + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_metric.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,56 @@ +% Compute metric scale factors. +% +% [pm,pn] = divand_metric(lon,lat) +% +% Compute metric scale factors pm and pn based on +% longitude lon and latitude lat. The variables pm and pn +% represent the inverse of the local resolution in meters using +% the mean Earth radius. + + + +function [pm,pn] = divand_metric(lon,lat) + +sz = size(lon); +i = 2:sz(1)-1; +j = 2:sz(2)-1; + + +dx = distance(lat(i-1,:),lon(i-1,:),lat(i+1,:),lon(i+1,:)); +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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_obs.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,83 @@ +% Include the constrain from the observations. +% +% s = divand_obs(s,xi,x,lambda,I) +% +% Set observations of variational problem. +% It is assumed that the each coordinate depends only on one +% index. If this is not the case, then matrix I must be provided. +% +% Input: +% s: structure created by divand_background +% xi: coordinates of observations* +% x: coordinates of grid* +% lambda: signal-to-noise ratio of observations +% I (optional): fractional indexes of location of observation +% within the grid +% +% Output: +% s: structure to be used by divand_factorize +% +% Note: +% *these parameters can either be specified as a cell +% array of all dimenions: +% xi = {Xi,Yi,Zi} +% or as n+1 dimensional array + +function c = divand_obs(s,xi,x,yo,lambda,I) + +if nargin == 5 + I = []; +end + +xi = cat_cell_array(xi); +x = cat_cell_array(x); + +mask = s.mask; +iscyclic = s.iscyclic; +moddim = s.moddim; + + +if isempty(I) + I = localize_separable_grid(x,mask,xi); +end + +[H,out] = sparse_interp(mask,I,iscyclic); + +nout = sum(out); +if nout ~= 0 + fprintf(1,'Observations out of domain: %d\n',nout); +end + + +H = H * sparse_pack(mask)'; + +% iB is scaled such that diag(inv(iB)) is 1 far from the +% boundary + +if isscalar(lambda) + R = 1/lambda * speye(size(H,1)); +elseif isvector(lambda) + R = sparse_diag(lambda); +else + R = lambda; +end + +c.H = H; +s.out = out; +c.R = R; +c.yo = yo; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_operators.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,112 @@ +% Generate the gradient and Laplacian operators. +% +% s = divand_operators(mask,pmn,nu,iscyclic) +% +% Form sparse matrixes representing the gradient and Laplacian using +% finite differences +% +% Input: +% mask: binary mask delimiting the domain. 1 is inside and 0 outside. +% For oceanographic application, this is the land-sea mask. +% pmn: scale factor of the grid. +% nu: diffusion coefficient of the Laplacian +% +% Output: +% s: stucture containing +% s.Dx: cell array of the gradient +% s.D: Laplaciant +% s.sv: structure describing the state vector +% s.mask: land-sea mask +% s.WE: diagonal matrix where each element is the surface +% of a grid cell + + +function s = divand_operators(mask,pmn,nu,iscyclic,mapindex) + +pmn = cat_cell_array(pmn); + +% number of dimensions +n = size(pmn,ndims(pmn)); + +sz = size(mask); + +sv = statevector_init(mask); + +if ~isempty(mapindex) + % mapindex is unpacked and referers to unpacked indices + + % range of packed indices + % land point map to 1, but those points are remove by statevector_pack + i2 = statevector_unpack(sv,(1:sv.n)',1); + mapindex_packed = statevector_pack(sv,i2(mapindex)); + + % applybc*x applies the boundary conditions to x + i = [1:sv.n]'; + applybc = sparse(i,mapindex_packed(i),ones(1,sv.n),sv.n,sv.n); + + % a halo point is points which maps to a (different) interior point + % a interior point maps to itself + s.isinterior = (i == mapindex_packed(i)); + s.isinterior_unpacked = statevector_unpack(sv,s.isinterior); + + s.mapindex_packed = mapindex_packed; +end + +D = divand_laplacian(mask,pmn,nu,iscyclic); + +% XXX remove this WE +d = statevector_pack(sv,1./(prod(pmn,n+1))); +WE = sparse_diag(sqrt(d)); + +for i=1:n + S = sparse_stagger(sz,i,iscyclic(i)); + + % mask on staggered grid + ma = (S * mask(:) == 1); + s.mask_stag{i} = ma; + + d = sparse_pack(ma) * prod(S * reshape(pmn,numel(mask),n),2); + d = 1./d; + s.WEs{i} = sparse_diag(sqrt(d)); +end + +s.Dx = cell(n,1); +[s.Dx{:}] = sparse_gradient(mask,pmn,iscyclic); + + +if ~isempty(mapindex) + D = applybc * D * applybc; + WE = sparse_diag(s.isinterior) * WE; + + for i=1:n + S = sparse_stagger(sz,i,iscyclic(i)); + s.isinterior_stag{i} = sparse_pack(s.mask_stag{i}) * S * s.isinterior_unpacked(:); + + % the results of 's.Dx{i} * field' satisfies the bc is field does + % there is need to reapply the bc on the result + s.Dx{i} = s.Dx{i} * applybc; + end + + s.applybc = applybc; +end + +s.D = D; +s.sv = sv; +s.mask = mask; +s.WE = WE; +s.n = n; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_orthogonalize_modes.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,48 @@ +% Orthogonalize EOF modes. +% +% [U3,WE] = divand_orthogonalize_mode(mask,pmn,U) +% +% The modes U are orthogonalized using the mask and metric (pmn). WE is square +% root of norm. Diagonal of WE.^2 is the surface of corresponding cell (units: +% m). The output U3 represent the normalized EOFs (units: m-1) + +function [U3,WE] = divand_orthogonalize_modes(mask,pmn,U); + +pmn = cat_cell_array(pmn); +% number of dimensions +n = size(pmn,ndims(pmn)); + +sv = statevector_init(mask); +d = statevector_pack(sv,1./(prod(pmn,n+1))); + +% spare root of the norm +WE = sparse_diag(sqrt(d)); + +U2 = statevector_pack(sv,U); +U2 = WE * U2; + +r = size(U2,2); +[U3,S3,V3] = svds(U2,r); +%U3 = reshape(U3,[size(mask) r]); + +U3 = inv(WE) * U3; + +%assert(U3' * WE^2 * U3, eye(r),1e-8) + +U3 = statevector_unpack(sv,U3); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_pc_michol.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,50 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_pc_none.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,30 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_pc_sqrtiB.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,31 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_realistic_example.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,99 @@ +% A realistic example of divand in 2 dimensions +% with salinity observations in the Mediterranean Sea at 30m. +% +% Data is available at +% http://modb.oce.ulg.ac.be/mediawiki/index.php/Divand + +% load observations + +load data.dat + +% extract columns of data +lonobs = data(:,1); +latobs = data(:,2); +S = data(:,3); + +% load bathymetry +% bath is negative in water and positive in air +bat = ncread('diva_bath.nc','bat'); +lon = ncread('diva_bath.nc','lon'); +lat = ncread('diva_bath.nc','lat'); + +% compute grid metric +% pm and pn are in meters^-1 +[lon,lat] = ndgrid(lon,lat); +[pm,pn] = divand_metric(lon,lat); + +% compute mean and anomalies +Smean = mean(S); +Sanom = S - mean(S); + +% correlation length (in meters) +len = 200e3; +len = 100e3; +% signal-to-noise ratio +lambda = 50.5; + +% land-sea mask +% mask everything below 30 m +mask = bat < -30; + +% mask for the analysis +maska = mask; +% remove Atlantic +maska(1:60,130:end) = 0; +% remove Black Sea +maska(323:end,100:end) = 0; + +% make analysis +Sa = divand(maska,{pm,pn},{lon,lat},{lonobs,latobs},Sanom,len,lambda); + +% add mean back +Sa2 = Sa + Smean; + +% create plots + +% color axis +ca = [36.2520 39.4480]; + +% aspect ratio +ar = [1 cos(mean(lat(:)) * pi/180) 1]; + + +subplot(2,1,1); +scatter(lonobs,latobs,10,S) +caxis(ca); +colorbar +hold on, contour(lon,lat,mask,0.5,'k') +xlim(lon([1 end])) +ylim(lat([1 end])) +title('Observations'); +set(gca,'DataAspectRatio',ar,'Layer','top') + +subplot(2,1,2); +pcolor(lon,lat,Sa2), shading flat,colorbar +hold on +plot(lonobs,latobs,'k.','MarkerSize',1) +caxis(ca) +contour(lon,lat,mask,0.5,'k') +xlim(lon([1 end])) +ylim(lat([1 end])) +title('Analysis'); +set(gca,'DataAspectRatio',ar,'Layer','top') + +set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 15 15 ]); +print -dpng divand_realistic_example.png +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_rms.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,37 @@ +% The RMS error between two variables. +% +% r = divand_rms(x,y,norm) +% +% Returns rms between x and y (taking norm into account if present) + +% Alexander Barth +function r = divand_rms(x,y,norm) + +d = x-y; + +if nargin == 3 + d = d .* sqrt(norm); +end + +m = ~isnan(d); +r = mean(d(m).^2); + +if nargin == 3 + r = r/mean(norm(m)); +end + +r = sqrt(r); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_simple_example.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,68 @@ +% A simple example of divand in 2 dimensions +% with observations from an analytical function. + + +% observations +x = rand(75,1); +y = rand(75,1); +f = sin(x*6) .* cos(y*6); + +% final grid +[xi,yi] = ndgrid(linspace(0,1,30)); + +% reference field +fref = sin(xi*6) .* cos(yi*6); + +% all points are valid points +mask = ones(size(xi)); + +% this problem has a simple cartesian metric +% pm is the inverse of the resolution along the 1st dimension +% pn is the inverse of the resolution along the 2nd dimension + +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +% correlation length +len = 0.1; + +% signal-to-noise ratio +lambda = 20; + +% fi is the interpolated field +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,len,lambda); + +% plotting of results +subplot(1,2,1); +pcolor(xi,yi,fref); +shading flat,colorbar +ca = caxis; +set(findobj(gcf,'FontSize',10),'FontSize',16) +title('Reference field and observation loc.'); +hold on +plot(x,y,'k.'); + +subplot(1,2,2); +hold on +pcolor(xi,yi,fi); +shading flat,colorbar +caxis(ca); +set(findobj(gcf,'FontSize',10),'FontSize',16) +title('Interpolated field'); + +set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 25 12 ]); +print -dpng divand_simple_example.png +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_solve.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,87 @@ +% Solve the variational problem. +% +% fi = divand_solve(s,yo) +% +% Derive the analysis based on all contraints included in s and using the +% observations yo +% +% Input: +% s: structure created by divand_factorize +% yo: value of the observations +% +% Output: +% fi: analyzed field + +function [fi,s] = divand_solve(s,yo) + +H = s.H; +sv = s.sv; +R = s.R; +% fix me: ignore argument +yo = s.yo; + + +if s.primal + if strcmp(s.inversion,'chol') + P = s.P; + fpi = P * (H'* (R \ yo(:))); + else + %H = double(H); + HiRyo = H'* (R \ yo(:)); + + fun = @(x) s.iB*x + H'*(R\ (H * x)); + + if ~s.keepLanczosVectors + [fpi,s.flag,s.relres,s.iter] = pcg(fun,HiRyo,s.tol,s.maxit,s.M1,s.M2); + + if s.flag ~= 0 + warning('divand-noconvergence','pcg did not converge'); + end + + else + %pc = @(x) s.M2 \ (s.M1 \ x); + pc = []; + x0 = zeros(size(s.iB,1),1); + [fpi,Q,T,diag] = conjugategradient(fun,HiRyo,'tol',s.tol,... + 'maxit',s.maxit,... + 'minit',s.minit,... + 'x0',x0,'renorm',1,'pc',pc); + s.iter = diag.iter; + s.relres = diag.relres; + s.P = CovarLanczos(Q,T); + end + end +else % dual + B = s.B; + % fun(x) computes (H B H' + R)*x + fun = @(x) H * (B * (H'*x)) + R*x; + + [tmp,flag,relres,iter] = pcg(fun, yo,s.tol,s.maxit,s.funPC); + + if (flag ~= 0) + error('divand:pcg', ['Preconditioned conjugate gradients method'... + ' did not converge %d %g %g'],flag,relres,iter); + end + + fpi = B * (H'*tmp); +end + + +[fi] = statevector_unpack(sv,fpi); + +fi(~s.mask) = NaN; + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/divand_solve_eof.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,73 @@ +% Solve the variational problem with the contraints from EOFs. +% +% [fi,s] = divand_solve_eof(s,yo,EOF_lambda,EOF) +% +% Derive the analysis based on all contraints included in s and using the +% observations yo and EOFs +% +% Input: +% s: structure created by divand_factorize +% yo: value of the observations +% EOF_lambda: weight of each EOF (adimentional) +% EOF: matrix containing the EOFs (units m^(-n/2)) +% +% Output: +% fi: analyzed field +% s.E: scaled EOFs + + +function [fi,s] = divand_solve_eof(s,yo,EOF_lambda,EOF) + +H = s.H; +sv = s.sv; +R = s.R; +P = s.P; +E = s.E; + + +r = size(E,2); + +%E = zeros(s.sv.n,r); +%for i=1:r +% E(:,i) = statevector_pack(sv,EOF(:,i)); +%end + +yo2 = (H'* (R \ yo(:))); + +% "classical analysis" +xa = P * yo2; + +% apply Pa to all EOFs +PaE = P * E; + +%keyboard + +M = E'*PaE - eye(r); + +% analysis with EOFs +xa = xa - PaE * pinv(M) * (E'*xa); + +[fi] = statevector_unpack(sv,xa); + +fi(~s.mask) = NaN; + + +% debugging information +s.cond = cond(M); +%lam = eig(M); +%s.lambda_min = min(lam); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/interp_regular.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,39 @@ +% Interpolation matrix for a n-dimensional interpolation. +% +% [H,out] = interp_regular(mask,x,xi) +% +% Creates sparse interpolation matrix H for n-dimensional interpolation problem. +% +% Input: +% mask: binary mask delimiting the domain. 1 is inside and 0 outside. +% For oceanographic application, this is the land-sea mask. +% x: cell array with n elements. Every element represents a coordinate of the +% observations +% xi: cell array with n elements. Every element represents a coordinate of the +% final grid on which the observations are interpolated. +% +% Output: +% H: sparse matrix representing the linear interpolation +% out: 1 for each observation out of the grid, 0 otherwise + +function [H,out] = interp_regular(mask,x,xi) + +%I = localize_regular_grid(xi,mask,x); +I = localize_separable_grid(xi,mask,x); +[H,out] = sparse_interp(mask,I); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/localize_regular_grid.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,67 @@ +% Derive fractional indices on a regular grid. +% +% I = localize_regular_grid(xi,mask,x) +% +% Derive fractional indices where xi are the points to localize in the +% regular grid x (constant increment in all dimension). +% The output I is an n-by-m array where n number of dimensions and m number of +% observations + +function I = localize_regular_grid(xi,mask,x) + +x = cat_cell_array(x); +xi = cat_cell_array(xi); + +% m is the number of arbitrarily distributed observations +tmp = size(xi); +mi = prod(tmp(1:end-1)); + +% sz is the size of the grid +tmp = size(x); +sz = tmp(1:end-1); +m = prod(sz); + +% n dimension of the problem +n = tmp(end); + +xi = reshape(xi,mi,n); + +scale = [1 cumprod(sz(1:end-1))]; + +A = zeros(n,n); + +x0 = x([0:n-1]' * m +1)'; + +E = eye(n); + +I = zeros(n,size(xi,1)); + +for i=1:n + I(i,:) = (xi(:,i) - x0(i)); + + % linear index of element (2,1,1,...,:), + % (1,2,1,...,:), (1,1,2,...,:) ... + + l = scale*E(:,i) + [0:n-1]' * m +1; + A(i,:) = x(l)- x0(:); +end + + +I = A \ I + 1; + +% LocalWords: indices sz + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/localize_separable_grid.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,90 @@ +% Derive fractional indices on a separable grid. +% +% I = localize_separable_grid(xi,mask,x) +% +% Derive fractional indices where xi are the points to localize in the +% separable grid x (every dimension in independent on other dimension). +% The output I is an n-by-m array where n number of dimensions and m number of +% observations + +function I = localize_separable_grid(xi,mask,x) + +x = cat_cell_array(x); +xi = cat_cell_array(xi); + +% m is the number of arbitrarily distributed observations +tmp = size(xi); +mi = prod(tmp(1:end-1)); + +% sz is the size of the grid +tmp = size(x); +sz = tmp(1:end-1); + +if isvector(x) + n = 1; + mi = length(xi); + I = zeros(1,mi); + I(1,:) = interp1(x,1:length(x),xi); +else + % n dimension of the problem + n = tmp(end); + m = prod(sz); + + xi = reshape(xi,mi,n); + x = reshape(x,[m n]); + + IJ = cell(n,1); + vi = cell(n,1); + X = cell(n,1); + XI = cell(n,1); + I = zeros(n,mi); + + for i=1:n + vi{i} = 1:sz(i); + X{i} = reshape(x(:,i),sz); + XI{i} = xi(:,i); + end + + [IJ{:}] = ndgrid(vi{:}); + + + for i=1:n + I(i,:) = interpn(X{:},IJ{i},XI{:}); + end + +end + + +% handle rounding errors +% snap to domain bounding box if difference does not exceeds tol +tol = 50*eps; + +for i=1:n + % upper bound + ind = sz(i) < I(i,:) & I(i,:) <= sz(i) + tol; + I(i,ind) = sz(i); + + % lower bound + ind = 1 < I(i,:) & I(i,:) <= 1 + tol; + I(i,ind) = 1; +end + + + + +% LocalWords: indices sz tol + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/mtimescorr.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,46 @@ +% Product between a Gaussian covariance matrix and a vector. +% +% p = mtimescorr(xi,len,b) +% +% Compute the product between a Gaussian covariance with a length scale len and +% grid points located in xi and the vector b. The covariance matrix has a unit +% variance. + +function p = mtimescorr(xi,len,b) + +%'mat' +m = size(xi,1); +p = zeros(size(b)); + +% for j=1:m +% for i=1:m +% d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 ); +% p(i) = p(i) + exp(-d2 / self.len^2) * b(j); +% end +% end +% + + +for i=1:m + X = repmat(xi(i,:),[m 1]); + d2 = sum( (X - xi).^2, 2) ; + ed = exp(-d2 / len^2); + p(i,:) = ed' * b; +end + +% LocalWords: mtimescorr len + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/private/geomean.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,19 @@ +% for compatibility with octave + +function m = geomean(x,varargin) + +m = prod(x,varargin{:})^(1/length(x)); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_diag.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,26 @@ +% Create diagonal sparse matrix. +% +% s = sparse_diag(d) +% +% Create diagonal sparse matrix s based on the diagonal elements d. The +% parameter d is transformed into a vector. + +function s = sparse_diag(d) + +s = spdiags(d(:),0,numel(d),numel(d)); + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_diff.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,81 @@ +% Sparse operator for differentiation. +% +% diffx = sparse_diff(sz1,m,cyclic) +% +% Sparse operator for differentiation along dimension m for "collapsed" matrix +% of the size sz1. +% +% Input: +% sz1: size of rhs +% m: dimension to differentiate +% cyclic: true if domain is cyclic along dimension m. False is the +% default value + +function diffx = sparse_diff(sz1,m,cyclic) + +if nargin == 2 + cyclic = false; +end + +n1 = prod(sz1); +sz2 = sz1; + +if ~cyclic + sz2(m) = sz2(m)-1; +end + +n2 = prod(sz2); + +n = length(sz1); + +vi = cell(1,n); +for i=1:n + vi{i} = 1:sz2(i); +end + +IJ = cell(n,1); + +[IJ{:}] = ndgrid(vi{:}); + +for i=1:n + IJ{i} = IJ{i}(:); +end + +L1 = [1:n2]'; + +L2 = sub2ind(sz1,IJ{:}); +one = ones(size(L1)); + +IJ{m} = IJ{m} + 1; + +if cyclic + IJ{m} = mod(IJ{m}-1,sz1(m))+1; +end + + +L2o = sub2ind(sz1,IJ{:}); + +diffx = sparse( ... + [L1; L1; ], ... + [L2; L2o; ], ... + [-one; one ], n2 , n1 ); + + + +% Copyright (C) 2009,2012 Alexander Barth <a.barth@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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_gradient.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,59 @@ +% Sparse operator for a gradient. +% +% [Dx1,Dx2,...,Dxn] = sparse_gradient(mask,pmn) +% +% Form the gradient using finite differences in all n-dimensions +% +% Input: +% mask: binary mask delimiting the domain. 1 is inside and 0 outside. +% For oceanographic application, this is the land-sea mask. +% +% pmn: scale factor of the grid. +% +% Output: +% Dx1,Dx2,...,Dxn: sparce matrix represeting a gradient along +% different dimensions + +function varargout = sparse_gradient(mask,pmn,iscyclic) + +H = sparse_pack(mask); + +sz = size(mask); +n = size(pmn,ndims(pmn)); +Pmn = reshape(pmn,[prod(sz) n]); + +if nargin == 2 + iscyclic = zeros(n,1); +end + +for i=1:n + % staggering operator + S = sparse_stagger(sz,i,iscyclic(i)); + + % mask for staggered variable + m = (S * mask(:) == 1); + + d = m .* (S * Pmn(:,i)); + + varargout{i} = sparse_pack(m) * sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)) * H'; +% varargout{i} = sparse_diag(d) * sparse_diff(sz,i) * H'; +% size(varargout{i}) +end + + + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_interp.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,127 @@ +% Create a sparse interpolation matrix. +% +% [H,out] = sparse_interp(mask,I) +% +% Create interpolation matrix from mask and fractional indexes I +% +% Input: +% mask: 0 invalid and 1 valid points (n-dimensional array) +% I: fractional indexes (2-dim array n by mi, where mi is the number of points to interpolate) +% Ouput: +% H: sparse matrix with interpolation coefficients +% out: true if value outside of grid + +function [H,out] = sparse_interp(mask,I,iscyclic) + +if ndims(mask) ~= size(I,1) && ~(isvector(mask) && size(I,1) == 1) + error('sparse_interp: inconsistent arguments') +end + +if isvector(mask) + sz = numel(mask); + mask = reshape(mask,[sz 1]); +else + sz = size(mask); +end + +% n dimension of the problem +n = size(I,1); +m = prod(sz); + +if nargin == 2 + iscyclic = zeros(1,n); +end + + +% mi is the number of arbitrarly distributed observations +mi = size(I,2); + +% handle cyclic dimensions +for i = 1:n + if iscyclic(i) + % bring I(i,:) inside the interval [1 sz(i)+1[ + % since the i-th dimension is cyclic + + I(i,:) = mod(I(i,:)-1,sz(i))+1; + end +end + + +scale = [1 cumprod(sz(1:end-1))]; + +% integer index +ind = floor(I); +inside = true(1,mi); + +for i = 1:n + if ~iscyclic(i) + % make a range check only for non-cyclic dimension + + % handle border cases + p = find(I(i,:) == sz(i)); + ind(i,p) = sz(i)-1; + + inside = inside & 1 <= ind(i,:) & ind(i,:) < sz(i); + end +end + +%whos ind + +% interpolation coefficient +coeff(:,:,2) = I - ind; +coeff(:,:,1) = 1 - coeff(:,:,2); + +% consider now only points inside +iind = find(inside); +coeff2 = coeff(:,iind,:); +ind2 = ind(:,iind); +mip = length(iind); % number of obs. inside + +si = repmat(1:mip,[2^n 1]); +sj = ones(2^n, mip); +ss = ones(2^n, mip); + +% loop over all corner of hypercube +for i=1:2^n + + % loop over all dimensions + for j=1:n + bit = bitget(i,j); + + % the j-index of the i-th corner has the index ip + % this index ip is zero-based + ip = (ind2(j,:) + bit - 1); + + % ip must be [0 and sz(j)-1] (zero-based) + % we know aleady that the point is inside the domain + % so, if it is outside this range then it is because of periodicity + ip = mod(ip,sz(j)); + + sj(i,:) = sj(i,:) + scale(j) * ip; + ss(i,:) = ss(i,:) .* coeff2(j,:,bit + 1); + end +end + +% sj must refer to a valid point or its interpolation coefficient +% must be zero + +inside(iind) = all(mask(sj) | ss == 0,1); +H = sparse(iind(si(:)),sj(:),ss(:),mi,m); + +out = ~inside; +% LocalWords: interp + +% Copyright (C) 2004 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_pack.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,35 @@ +% Create a sparse matrix which packs an array under the control of a mask. +% +% H = sparse_pack(mask) +% +% Create a sparse matrix H which packs an array. It selects the elements where +% mask is true. + +function H = sparse_pack(mask) + + +j = find(mask)'; +m = length(j); +i = 1:m; +s = ones(1,m); + +n = numel(mask); + +%whos i j s +H = sparse(i,j,s,m,n); + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_shift.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,75 @@ +% Sparse operator shifting a field in a given dimension. +% +% function S = sparse_shift(sz1,m,cyclic) +% +% Sparse operator shifting a field in the dimension m. The field is a +% "collapsed" matrix of the size sz1. +% +% Input: +% sz1: size of rhs +% m: dimension to shift +% cyclic: true if domain is cyclic along dimension m. False is the +% default value + +function S = sparse_shift(sz1,m,cyclic) + +if nargin == 2 + cyclic = false; +end + +n1 = prod(sz1); + +sz2 = sz1; + +if ~cyclic + sz2(m) = sz2(m)-1; +end + +n2 = prod(sz2); + +n = length(sz1); + +vi = cell(1,n); +for i=1:n + vi{i} = 1:sz2(i); +end + +IJ = cell(n,1); + +[IJ{:}] = ndgrid(vi{:}); + +for i=1:n + IJ{i} = IJ{i}(:); +end + +L1 = [1:n2]'; +one = ones(size(L1)); + +IJ{m} = IJ{m} + 1; + +if cyclic + IJ{m} = mod(IJ{m}-1,sz1(m))+1; +end + + +L2 = sub2ind(sz1,IJ{:}); + +S = sparse(L1,L2,one,n2,n1); + + + +% Copyright (C) 2012 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_stagger.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,78 @@ +% Sparse operator for staggering. +% +% S = sparse_stagger(sz1,m,cyclic) +% +% Create a sparse operator for staggering a field in dimension m. +% The field is a "collapsed" matrix of the size sz1. +% +% Input: +% sz1: size of rhs +% m: dimension to stagger +% cyclic: true if domain is cyclic along dimension m. False is the +% default value + +function S = sparse_stagger(sz1,m,cyclic) + + +if nargin == 2 + cyclic = false; +end + +n1 = prod(sz1); + +sz2 = sz1; + +if ~cyclic + sz2(m) = sz2(m)-1; +end +n2 = prod(sz2); + +n = length(sz1); + +for i=1:n + vi{i} = 1:sz2(i); +end + +IJ = cell(n,1); + +[IJ{:}] = ndgrid(vi{:}); + +for i=1:n + IJ{i} = IJ{i}(:); +end + +L1 = [1:n2]'; + +L2 = sub2ind(sz1,IJ{:}); +v = ones(size(L1))/2; + +IJ{m} = IJ{m} + 1; + +if cyclic + IJ{m} = mod(IJ{m}-1,sz1(m))+1; +end + +L2o = sub2ind(sz1,IJ{:}); + +S = sparse( ... + [L1; L1; ], ... + [L2; L2o; ], ... + [v; v ], n2 , n1 ); + + + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/sparse_trim.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,61 @@ +% Sparse operator for trimming. +% +% T = sparse_trim(sz1,m) +% +% Create a sparse operator which trim first and last row (or column) in +% The field is a "collapsed" matrix of the size sz1. +% +% Input: +% sz1: size of rhs +% m: dimension to trim + +function T = sparse_trim(sz1,m) + +n1 = prod(sz1); + +sz2 = sz1; +sz2(m) = sz2(m)-2; +n2 = prod(sz2); + +n = length(sz1); + +% L1 + +for i=1:n + vi{i} = 1:sz2(i); +end + +IJ = cell(n,1); + +[IJ{:}] = ndgrid(vi{:}); + +for i=1:n + IJ{i} = IJ{i}(:); +end + +L1 = sub2ind(sz2,IJ{:}); + +IJ{m}=IJ{m}+1; + +L2 = sub2ind(sz1,IJ{:}); + +one = ones(size(L1)); + +T = sparse(L1, L2, one, n2, n1); + + +% Copyright (C) 2009 Alexander Barth <a.barth@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/>. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/statevector_init.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,74 @@ +% Initialize structure for packing and unpacking given their mask. +% +% s = statevector_init(mask1, mask2, ...) +% +% Initialize structure for packing and unpacking +% multiple variables given their corresponding land-sea mask. +% +% Input: +% mask1, mask2,...: land-sea mask for variable 1,2,... Sea grid points correspond to one and land grid points to zero. +% Every mask can have a different shape. +% +% Output: +% s: structure to be used with statevector_pack and statevector_unpack. +% +% Note: +% see also statevector_pack, statevector_unpack + +% Author: Alexander Barth, 2009 <a.barth@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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/statevector_pack.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,75 @@ +% Pack a series of variables into a vector under the control of a mask. +% +% x = statevector_pack(s,var1, var2, ...) +% +% Pack the different variables var1, var2, ... into the vector x. +% Only sea grid points are retained. +% +% Input: +% s: structure generated by statevector_init. +% var1, var2,...: variables to pack (with the same shape as the corresponding masks). +% +% Output: +% x: vector of the packed elements. The size of this vector is the number of elements of all masks equal to 1. +% +% Notes: +% If var1, var2, ... have an additional trailing dimension, then this dimension is assumed +% to represent the different ensemble members. In this case x is a matrix and its last dimension +% is the number of ensemble members. + +% Author: Alexander Barth, 2009 <a.barth@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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/statevector_unpack.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,64 @@ +% Unpack a vector into different variables under the control of a mask. +% +% [var1, var2, ...] = statevector_unpack(s,x) +% [var1, var2, ...] = statevector_unpack(s,x,fillvalue) +% +% Unpack the vector x into the different variables var1, var2, ... +% +% Input: +% s: structure generated by statevector_init. +% x: vector of the packed elements. The size of this vector is the number of elements equal to 1 +% in all masks. +% +% Optional input parameter: +% fillvalue: The value to fill in var1, var2,... where the masks correspond to a land grid point. The default is zero. +% +% Output: +% var1, var2,...: unpacked variables. +% +% Notes: +% If x is a matrix, then the second dimension is assumed +% to represent the different ensemble members. In this case, +% var1, var2, ... have also an additional trailing dimension. + +% Author: Alexander Barth, 2009 <a.barth@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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_1dvar.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,30 @@ +% Testing divand in 1 dimension. + +% grid of background field +[xi] = linspace(0,1,30)'; + +x = [.4 .6]'; +f = [.4 .6]'; + +mask = ones(size(xi)); +mask([1 end]) = 0; + +pm = ones(size(xi)) / (xi(2)-xi(1)); + +[fi,err] = divand(mask,pm,{xi},{x},f,.1,2); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,58 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_adv.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,42 @@ +% Testing divand in 2 dimensions with advection. + +% grid of background field +[xi,yi] = ndgrid(linspace(-1,1,30)); + +x = .4; +y = .4; +f =1; + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +u = ones(size(xi)); +v = u; + +a = 5; +u = a*yi; +v = -a*xi; + +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,.2,200,'velocity',{u,v}); + +if abs( fi(18,24) - 0.89935) > 1e-3 + error('unexpected large difference'); +end + + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_check.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,64 @@ +% Testing divand in 2 dimensions with independent verification. + +% grid of background field +[xi,yi] = ndgrid(linspace(0,1,10)); + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +% grid of observations +[x,y] = ndgrid(linspace(eps,1-eps,10)); +x = x(:); +y = y(:); +v = sin( x*6 ) .* cos( y*6); + + +lenx = .15; +leny = .15; + +lambda = 20; + +[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,{lenx,leny},lambda,'diagnostics',1,'primal',1); + + + +iR = inv(full(s.R)); +iB = full(s.iB); +H = full(s.H); +sv = s.sv; + +iP = iB + H'*iR*H; + +P = inv(iP); + +xa2 = P* (H'*iR*v(:)); + +[fi2] = statevector_unpack(sv,xa2); +fi2(~s.mask) = NaN; + +rms_diff = []; +rms_diff(end+1) = divand_rms(va,fi2); +rms_diff(end+1) = divand_rms(diag(P),err(:)); + + +if any(rms_diff > 1e-6) + error('unexpected large difference'); +end + +fprintf('(max rms=%g) ',max(rms_diff)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_check_correrr.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,70 @@ +% Testing divand in 2 dimensions with correlated errors. + +% grid of background field +[xi,yi] = ndgrid(linspace(0,1,10)); + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +% grid of observations +[x,y] = ndgrid(linspace(eps,1-eps,3)); +x = x(:); +y = y(:); +v = sin( x*6 ) .* cos( y*6); + +xx = repmat(x,[1 numel(x)]); +yy = repmat(y,[1 numel(y)]); + +len = 0.1; +variance = 0.1; +R = variance * exp(- ((xx - xx').^2 + (yy - yy').^2 ) / len^2); +RC = CovarParam({x,y},variance,len); + +len = .1; +lambda = RC; + +[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',0); + +iR = inv(full(R)); +iB = full(s.iB); +H = full(s.H); +sv = s.sv; + +iP = iB + H'*iR*H; + +Pa = inv(iP); + +xa2 = Pa* (H'*iR*v(:)); + +[fi2] = statevector_unpack(sv,xa2); +fi2(~s.mask) = NaN; + +rms_diff = divand_rms(va,fi2); + +if rms_diff > 1e-6 + disp('FAILED'); +end + +rms_diff = divand_rms(err(:),diag(Pa)); + +if rms_diff > 1e-6 + error('unexpected large difference'); +end + +fprintf('(max rms=%g) ',max(rms_diff)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_constrain.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,84 @@ +% Testing divand in 2 dimensions with a custom constrain. + +% grid of background field +[xi,yi] = ndgrid(linspace(0,1,10)); + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +% grid of observations +[x,y] = ndgrid(linspace(eps,1-eps,10)); +x = x(:); +y = y(:); +v = sin( x*6 ) .* cos( y*6); + + +len = 3; +lambda = 20; +n = sum(mask(:)); + +% constrain +fun = @(x) mean(x,1); +funt = @(x) repmat(x,[n 1])/n; + +H = ones(1,n)/n; +c.H = MatFun([1 n],fun,funt); +%c.H = H; +c.R = 0.001; +c.yo = 1; + +% bug +%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0,'constraint',c,'inversion','pcg'); + +[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,... + 'diagnostics',1,'primal',1,'factorize',0,'inversion','pcg',... + 'tol',1e-6,'minit',100,'keepLanczosVectors',1); + + +% ok +%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0); + + + + +iR = inv(full(s.R)); +iB = full(s.iB); +H = double(s.H); +yo = s.yo; +sv = s.sv; + +iB = iB + H'*iR*H; + +Pa = inv(iB); + +xa2 = Pa* (H'*iR*yo); + +[fi2] = statevector_unpack(sv,xa2); +fi2(~s.mask) = NaN; + +rms_diff = []; +rms_diff(end+1) = sqrt(mean((va(:) - fi2(:)).^2)); +rms_diff(end+1) = sqrt(mean((diag(Pa) - err(:)).^2)); + + +if any(rms_diff > 1e-4) + error('unexpected large difference'); +end + +fprintf('(max rms=%g) ',max(rms_diff)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_cyclic.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,97 @@ +% Testing divand in 2 dimensions in a cyclic domain. + +% grid of background field +n = 30; + +[xi,yi] = ndgrid(0:n-1); + +x = xi(2,15); +y = yi(2,15); +f = 1; + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + + + +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0]); +rms_diff = []; + +rms_diff(end+1) = divand_rms(fi(end,:),fi(4,:)); + + +x = xi(2,2); +y = yi(2,2); + +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n n]); + +rms_diff(end+1) = abs(fi(4,4) - fi(end,end)); + + +fi2 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.]); +fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[30.5; 15.]); +rms_diff(end+1) = divand_rms(fi2,fi3); + + + +fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.]); +fi4 = cat(1,fi3(11:end,:),fi3(1:10,:)); + +rms_diff(end+1) = divand_rms(fi2,fi4); + + + + + +% test with mask +mask(11,4) = 0; +mask(3:15,20) = 0; + +fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.]); +fi4 = cat(1,fi3(11:end,:),fi3(1:10,:)); + +mask = cat(1,mask(11:end,:),mask(1:10,:)); + +fi2 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.]); + + +rms_diff(end+1) = divand_rms(fi2,fi4); +rms_diff(end+1) = divand_rms(isnan(fi2),isnan(fi4)); + +% advection + +a = 5; +u = a*cos(2*pi*yi/n); +v = a*cos(2*pi*xi/n); + +mask = ones(size(xi)); +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[10.5; 15.],'velocity',{u,v}); +fi2 = cat(1,fi(11:end,:),fi(1:10,:)); + +ut = cat(1,u(11:end,:),u(1:10,:)); +vt = cat(1,v(11:end,:),v(1:10,:)); +fi3 = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,200,'moddim',[n 0],'fracindex',[0.5; 15.],'velocity',{ut,vt}); + +rms_diff(end+1) = divand_rms(fi2,fi3); + +if any(rms_diff > 1e-10) || any(isnan(rms_diff)) + error('unexpected large difference'); +end + +fprintf('(max rms=%g) ',max(rms_diff)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_eof_check.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,97 @@ +% Testing divand in 2 dimensions with EOF constraints. + +[xi,yi] = ndgrid(linspace(0,1,30)); + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + + +x = .5; +y = .5; +v = 1; + +len = .1; +lambda = 10; + +clear U2 +U2(:,1) = mask(:); +U2(:,2) = xi(:); +U2(:,3) = yi(:); + +%U2 = randn(numel(mask),10); +r = size(U2,2); +[U3,S3,V3] = svds(U2,r); +U3 = reshape(U3,[size(mask) r]); + +[va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1); +beta = 10; + +[va_eof,s_eof] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'EOF',U3,'EOF_scaling',beta*ones(r,1),'diagnostics',1); + +[Jb,Jo,Jeof] = divand_diagnose(s,va_eof,v); + +%assert(s_eof.Jb,Jb,1e-6) -> should no longer be the case since s.iB is modified to ensure that B - EE is positive defined +assert(abs(s_eof.Jo-Jo) < 1e-6) + + + +s = s_eof; + +iR = inv(full(s.R)); +iB = s.iB; +H = s.H; +E = s.E; +sv = s.sv; + +iP = iB + H'*iR*H; + +Pa = inv(iP - E*E'); + +xa2 = Pa* (H'*iR*v(:)); + +[fi2] = statevector_unpack(sv,xa2); +fi2(~s.mask) = NaN; + +rms_diff = divand_rms(va_eof,fi2); + +if rms_diff > 1e-6 + error('unexpected large difference'); +end + +Pf = inv(iB - E*E'); +R = inv(iR); +xa3 = Pf*H' * inv(H*Pf*H' + R) * v; + +[fi3] = statevector_unpack(sv,xa3); +fi3(~s.mask) = NaN; + +rms_diff = divand_rms(va_eof,fi3); + +if rms_diff > 1e-6 + error('unexpected large difference'); +end + +[U,L] = eig(Pf); +L = real(diag(L)); + +if any(L < 0) + error('some eingenvalues are negative'); +end + +Pfvar = inv(iB); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_lenxy.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,54 @@ +% Testing divand in 2 dimensions with lenx /= leny. + +n2 = 15; +% grid of background field +[xi,yi] = ndgrid(linspace(-1,1,2*n2+1)); + +x = 0; +y = 0; +f = 1; + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +u = ones(size(xi)); +v = u; + +a = 5; +u = a*yi; +v = -a*xi; + +L = 0.2 * ones(size(xi)); +%L = 0.2 * ones(size(xi)); +L = 0.4 * (2 - xi); +L = ones([size(mask) 2]); +% must be L(:,:,1) < L(:,:,2) + +L(:,:,1) = 0.2; +L(:,:,2) = 0.5; + +Lx = 0.2 * ones(size(mask)); +Ly = 0.5 * ones(size(mask)); + + +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,{Lx,Ly},2); + +if fi(n2,1) < fi(1,n2) + error('unexpected values'); +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_pcg.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,65 @@ +% Testing divand in 2 dimensions with pcg solver. + +% grid of background field +[xi,yi] = ndgrid(linspace(0,1,20)); + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + +% grid of observations +[x,y] = ndgrid(linspace(eps,1-eps,10)); +x = x(:); +y = y(:); +v = sin( x*6 ) .* cos( y*6); + + +len = 0.1; +lambda = 1; +n = sum(mask(:)); + +% constrain +fun = @(x) mean(x,1); +funt = @(x) repmat(x,[n 1])/n; + +H = ones(1,n)/n; +c.H = MatFun([1 n],fun,funt); +%c.H = H; +c.R = 0.001; +c.yo = 1; + +% bug +%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0,'constraint',c,'inversion','pcg'); + +[va0,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,... + 'diagnostics',1,'primal',1); + +[va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,... + 'diagnostics',1,'primal',1,'inversion','pcg',... + 'tol',1e-6,'maxit',10000); + + +rms_diff = []; +rms_diff(end+1) = sqrt(mean((va(:) - va0(:)).^2)); + + +if any(rms_diff > 1e-4) + error('unexpected large difference'); +end + +fprintf('(max rms=%g) ',max(rms_diff)); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_2dvar_rellen.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,47 @@ +% Testing divand in 2 dimensions with relative correlation length. + +% grid of background field +[xi,yi] = ndgrid(linspace(-1,1,30)); + +x = 0; +y = 0; +f = 1; + +mask = ones(size(xi)); +pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); +pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); + + + +u = ones(size(xi)); +v = u; + +a = 5; +u = a*yi; +v = -a*xi; + +L = 0.2 * ones(size(xi)); +%L = 0.2 * ones(size(xi)); +L = 0.4 * (2 - xi); + +fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,L,200); + + +if any(fi(1,:) < fi(end,:)) + error('unexpected values'); +end + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_3dvar.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,54 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_3dvar_large_stacked.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,71 @@ +% 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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_4dvar.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,53 @@ +% Testing divand in 4 dimensions. + +% grid of background field +[xi,yi,zi,ti] = ndgrid(linspace(0,1,7)); +len = 1; +k = pi/2; +k = pi; +fun = @(x,y,z,t) sin(k*x) .* sin(k*y) .* sin(k*z) .* sin(k*t); +fi_ref = fun(xi,yi,zi,ti); + +% grid of observations +[x,y,z,t] = ndgrid(linspace(2*eps,1-2*eps,5)); +x = x(:); +y = y(:); +z = z(:); +t = t(:); + +pm = ones(size(xi)) / (xi(2,1,1,1)-xi(1)); +pn = ones(size(xi)) / (yi(1,2,1,1)-yi(1)); +po = ones(size(xi)) / (zi(1,1,2,1)-zi(1)); +pp = ones(size(xi)) / (ti(1,1,1,2)-ti(1)); + + +f = fun(x,y,z,t); +mask = ones(size(xi)); + +%f = x = y = z = t = 0.5 + +fi = divand(mask,{pm,pn,po,pp},{xi,yi,zi,ti},{x,y,z,t},f,.12,10); +rms = divand_rms(fi_ref,fi); + +% rms should be 0.0111807 + +if (rms > 0.012) + error('unexpected large difference with reference field'); +end + +fprintf('(max difference=%g) ',rms); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_divand.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,55 @@ +% Test if divand is working correctly. + +tests = {... + 'test_interp_1d',... + 'test_interp_2d',... + 'test_interp_regular',... + 'test_sparse_diff',... + 'test_1dvar',... + 'test_2dvar',... + 'test_2dvar_check',... + 'test_2dvar_adv',... + 'test_2dvar_rellen',... + 'test_2dvar_lenxy',... + 'test_2dvar_check_correrr',... + 'test_2dvar_pcg',... + 'test_2dvar_constrain',... + 'test_2dvar_cyclic',... + 'test_2dvar_eof_check',... + 'test_3dvar',... + 'test_3dvar_large_stacked',... + 'test_4dvar'}; + +% disable warning for full matrice product +state = warning('query','divand:fullmat'); +warning('off','divand:fullmat') + +for iindex=1:length(tests); + + fprintf('run %s: ',tests{iindex}); + try + eval(tests{iindex}); + colordisp(' OK ','green'); + catch + colordisp(' FAIL ','red'); + disp(lasterr) + end +end + +% restore warning state +warning(state.state,'divand:fullmat') + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_interp_1d.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,37 @@ +% Testing 1D linear interpolation. + +[xi] = linspace(0,1,30); +fi = sin( xi*6 ); + +no = 20; +x = rand(no,1); + +mask = ones(size(xi)); + +[H,out] = interp_regular(mask,{xi},{x}); + +f1 = interp1(xi,fi,x); + +f2 = H * fi(:); +Difference = max(abs(f1(:) - f2(:))); + +if (Difference > 1e-6) + error('unexpected large difference with reference field'); +end + +fprintf('(max difference=%g) ',Difference); + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_interp_2d.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,39 @@ +% Testing 2D linear interpolation. + +[xi,yi] = ndgrid(linspace(0,1,30)); +fi = sin( xi*6 ) .* cos( yi*6); + +no = 20; +x = rand(no,1); +y = rand(no,1); + +mask = ones(size(xi)); + +[H,out] = interp_regular(mask,{xi,yi},{x,y}); + +f1 = interpn(xi,yi,fi,x,y); + +f2 = H * fi(:); +Difference = max(abs(f1(:) - f2(:))); + +if (Difference > 1e-6) + error('unexpected large difference with reference field'); +end + +fprintf('(max difference=%g) ',Difference); + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_interp_regular.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,83 @@ +% Testing linear interpolation on regular grid. + +D = []; + +[xi,yi,zi] = ndgrid(linspace(0,1,11)); +fi = sin(6*xi) .* cos(6*yi) .* sin(6*zi); + +no = 20; +x = rand(no,1); +y = rand(no,1); +z = rand(no,1); + +[H,out] = interp_regular(ones(size(xi)),{xi,yi,zi},{x,y,z}); +f1 = interpn(xi,yi,zi,fi,x,y,z); + +f2 = H * fi(:); +D(end+1) = max(abs(f1(:) - f2(:))); + + +% 4d +[xi,yi,zi,ti] = ndgrid(linspace(0,1,5)); +fi = sin(6*xi) .* cos(6*yi) .* sin(6*zi) .* cos(6*ti); +no = 20; +x = rand(no,1); +y = rand(no,1); +z = rand(no,1); +t = rand(no,1); + +[H,out] = interp_regular(ones(size(xi)),{xi,yi,zi,ti},{x,y,z,t}); +f1 = interpn(xi,yi,zi,ti,fi,x,y,z,t); + +f2 = H * fi(:); +D(end+1) = max(abs(f1(:) - f2(:))); + +% 2d with outside points + +[xi,yi] = ndgrid(linspace(0,1,30)); +fi = sin( xi*6 ) .* cos( yi*6); +fi(1:4,:) = NaN; + +no = 200; +x = 2*rand(no,1); +y = 2*rand(no,1); + + +%x = 0.805224734876956 +%y = 1.00822133283472 + +%x = 0.134579335939626 +%y = 0.716156681034316 + +f = sin( x*6 ) .* cos( y*6); + +mask = ~isnan(fi); + +[H,out] = interp_regular(mask,{xi,yi},{x,y}); + +f1 = interpn(xi,yi,fi,x,y); +f2 = H * fi(:); +f2 = reshape(f2,size(x)); + +D(end+1) = max(abs(f1(~out) - f2(~out))); + + +if any(D > 1e-6) || any(isnan(f1(:) ~= out(:))) + error('unexpected large difference with reference field'); +end + +fprintf('(max difference=%g) ',Difference); +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/divand/inst/test_sparse_diff.m Wed Feb 12 12:30:06 2014 +0000 @@ -0,0 +1,81 @@ +% Testing sparse operators. + +clear d +f = randn(21,30,10); + +S = sparse_diff(size(f),1); +f1 = f(2:end,:,:) - f(1:end-1,:,:); +f2 = S*f(:); +d(1) = max(abs(f1(:) - f2(:))); + +S = sparse_diff(size(f),2); +f1 = f(:,2:end,:) - f(:,1:end-1,:); +f2 = S*f(:); +d(2) = max(abs(f1(:) - f2(:))); + + +S = sparse_diff(size(f),3); +f1 = f(:,:,2:end) - f(:,:,1:end-1); +f2 = S*f(:); +d(3) = max(abs(f1(:) - f2(:))); + + +% cyclic + +% dim 1 cyclic +fc = cat(1,f,f(1,:,:)); +dfc = fc(2:end,:,:) - fc(1:end-1,:,:); + +S = sparse_diff(size(f),1,true); +f2 = S*f(:); +d(4) = max(abs(dfc(:) - f2(:))); + + +% dim 2 cyclic +fc = cat(2,f,f(:,1,:)); +f1 = fc(:,2:end,:) - fc(:,1:end-1,:); + +S = sparse_diff(size(f),2,true); +f2 = S*f(:); +d(5) = max(abs(f1(:) - f2(:))); + + +% stagger + +S = sparse_stagger(size(f),1); +f1 = (f(2:end,:,:) + f(1:end-1,:,:))/2; +f2 = S*f(:); +d(end+1) = max(abs(f1(:) - f2(:))); + + +% dim 1 cyclic +fc = cat(1,f,f(1,:,:)); +f1 = (fc(2:end,:,:) + fc(1:end-1,:,:))/2; + +S = sparse_stagger(size(f),1,true); +f2 = S*f(:); +d(end+1) = max(abs(f1(:) - f2(:))); + + +if any(d > 1e-6) + error('unexpected large difference with reference field'); +end + +fprintf('(max difference=%g) ',max(d)); + + + +% Copyright (C) 2014 Alexander Barth <a.barth@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/>.