--- a/extra/divand/CITATION	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,22 +0,0 @@
-Please cite divand as follows if you use this software in a publication:
-
-  Barth, A., Beckers, J.-M., Troupin, C., Alvera-Azcárate, A., and
-  Vandenbulcke, L.: divand-1.0: n-dimensional variational data analysis
-  for ocean observations, Geosci. Model Dev., 7, 225-241,
-  doi:10.5194/gmd-7-225-2014, 2014.  
-  URL http://www.geosci-model-dev.net/7/225/2014/
-
-A BibTeX entry for LaTeX users is:
-
-  @Article{gmd-7-225-2014,
-    AUTHOR  = {Barth, A. and Beckers, J.-M. and Troupin, C. and Alvera-Azc\'arate, A. and Vandenbulcke, L.},
-    TITLE   = {divand-1.0: n-dimensional variational data analysis for ocean observations},
-    JOURNAL = {Geoscientific Model Development},
-    VOLUME  = {7},
-    YEAR    = {2014},
-    NUMBER  = {1},
-    PAGES   = {225--241},
-    URL     = {http://www.geosci-model-dev.net/7/225/2014/},
-    DOI     = {10.5194/gmd-7-225-2014}
-  }
-

--- a/extra/divand/COPYING	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,337 +0,0 @@
-                    GNU GENERAL PUBLIC LICENSE
-                       Version 2, June 1991
-
- Copyright (C) 1989, 1991 Free Software Foundation, Inc. <http://fsf.org/>
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
-                            Preamble
-
-  The licenses for most software are designed to take away your
-freedom to share and change it.  By contrast, the GNU General Public
-License is intended to guarantee your freedom to share and change free
-software--to make sure the software is free for all its users.  This
-General Public License applies to most of the Free Software
-Foundation's software and to any other program whose authors commit to
-using it.  (Some other Free Software Foundation software is covered by
-the GNU Library General Public License instead.)  You can apply it to
-your programs, too.
-
-  When we speak of free software, we are referring to freedom, not
-price.  Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-this service if you wish), that you receive source code or can get it
-if you want it, that you can change the software or use pieces of it
-in new free programs; and that you know you can do these things.
-
-  To protect your rights, we need to make restrictions that forbid
-anyone to deny you these rights or to ask you to surrender the rights.
-These restrictions translate to certain responsibilities for you if you
-distribute copies of the software, or if you modify it.
-
-  For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must give the recipients all the rights that
-you have.  You must make sure that they, too, receive or can get the
-source code.  And you must show them these terms so they know their
-rights.
-
-  We protect your rights with two steps: (1) copyright the software, and
-(2) offer you this license which gives you legal permission to copy,
-distribute and/or modify the software.
-
-  Also, for each author's protection and ours, we want to make certain
-that everyone understands that there is no warranty for this free
-software.  If the software is modified by someone else and passed on, we
-want its recipients to know that what they have is not the original, so
-that any problems introduced by others will not reflect on the original
-authors' reputations.
-
-  Finally, any free program is threatened constantly by software
-patents.  We wish to avoid the danger that redistributors of a free
-program will individually obtain patent licenses, in effect making the
-program proprietary.  To prevent this, we have made it clear that any
-patent must be licensed for everyone's free use or not licensed at all.
-
-  The precise terms and conditions for copying, distribution and
-modification follow.
-
-                    GNU GENERAL PUBLIC LICENSE
-   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
-
-  0. This License applies to any program or other work which contains
-a notice placed by the copyright holder saying it may be distributed
-under the terms of this General Public License.  The "Program", below,
-refers to any such program or work, and a "work based on the Program"
-means either the Program or any derivative work under copyright law:
-that is to say, a work containing the Program or a portion of it,
-either verbatim or with modifications and/or translated into another
-language.  (Hereinafter, translation is included without limitation in
-the term "modification".)  Each licensee is addressed as "you".
-
-Activities other than copying, distribution and modification are not
-covered by this License; they are outside its scope.  The act of
-running the Program is not restricted, and the output from the Program
-is covered only if its contents constitute a work based on the
-Program (independent of having been made by running the Program).
-Whether that is true depends on what the Program does.
-
-  1. You may copy and distribute verbatim copies of the Program's
-source code as you receive it, in any medium, provided that you
-conspicuously and appropriately publish on each copy an appropriate
-copyright notice and disclaimer of warranty; keep intact all the
-notices that refer to this License and to the absence of any warranty;
-and give any other recipients of the Program a copy of this License
-along with the Program.
-
-You may charge a fee for the physical act of transferring a copy, and
-you may at your option offer warranty protection in exchange for a fee.
-
-  2. You may modify your copy or copies of the Program or any portion
-of it, thus forming a work based on the Program, and copy and
-distribute such modifications or work under the terms of Section 1
-above, provided that you also meet all of these conditions:
-
-    a) You must cause the modified files to carry prominent notices
-    stating that you changed the files and the date of any change.
-
-    b) You must cause any work that you distribute or publish, that in
-    whole or in part contains or is derived from the Program or any
-    part thereof, to be licensed as a whole at no charge to all third
-    parties under the terms of this License.
-
-    c) If the modified program normally reads commands interactively
-    when run, you must cause it, when started running for such
-    interactive use in the most ordinary way, to print or display an
-    announcement including an appropriate copyright notice and a
-    notice that there is no warranty (or else, saying that you provide
-    a warranty) and that users may redistribute the program under
-    these conditions, and telling the user how to view a copy of this
-    License.  (Exception: if the Program itself is interactive but
-    does not normally print such an announcement, your work based on
-    the Program is not required to print an announcement.)
-
-These requirements apply to the modified work as a whole.  If
-identifiable sections of that work are not derived from the Program,
-and can be reasonably considered independent and separate works in
-themselves, then this License, and its terms, do not apply to those
-sections when you distribute them as separate works.  But when you
-distribute the same sections as part of a whole which is a work based
-on the Program, the distribution of the whole must be on the terms of
-this License, whose permissions for other licensees extend to the
-entire whole, and thus to each and every part regardless of who wrote it.
-
-Thus, it is not the intent of this section to claim rights or contest
-your rights to work written entirely by you; rather, the intent is to
-exercise the right to control the distribution of derivative or
-collective works based on the Program.
-
-In addition, mere aggregation of another work not based on the Program
-with the Program (or with a work based on the Program) on a volume of
-a storage or distribution medium does not bring the other work under
-the scope of this License.
-
-  3. You may copy and distribute the Program (or a work based on it,
-under Section 2) in object code or executable form under the terms of
-Sections 1 and 2 above provided that you also do one of the following:
-
-    a) Accompany it with the complete corresponding machine-readable
-    source code, which must be distributed under the terms of Sections
-    1 and 2 above on a medium customarily used for software interchange; or,
-
-    b) Accompany it with a written offer, valid for at least three
-    years, to give any third party, for a charge no more than your
-    cost of physically performing source distribution, a complete
-    machine-readable copy of the corresponding source code, to be
-    distributed under the terms of Sections 1 and 2 above on a medium
-    customarily used for software interchange; or,
-
-    c) Accompany it with the information you received as to the offer
-    to distribute corresponding source code.  (This alternative is
-    allowed only for noncommercial distribution and only if you
-    received the program in object code or executable form with such
-    an offer, in accord with Subsection b above.)
-
-The source code for a work means the preferred form of the work for
-making modifications to it.  For an executable work, complete source
-code means all the source code for all modules it contains, plus any
-associated interface definition files, plus the scripts used to
-control compilation and installation of the executable.  However, as a
-special exception, the source code distributed need not include
-anything that is normally distributed (in either source or binary
-form) with the major components (compiler, kernel, and so on) of the
-operating system on which the executable runs, unless that component
-itself accompanies the executable.
-
-If distribution of executable or object code is made by offering
-access to copy from a designated place, then offering equivalent
-access to copy the source code from the same place counts as
-distribution of the source code, even though third parties are not
-compelled to copy the source along with the object code.
-
-  4. You may not copy, modify, sublicense, or distribute the Program
-except as expressly provided under this License.  Any attempt
-otherwise to copy, modify, sublicense or distribute the Program is
-void, and will automatically terminate your rights under this License.
-However, parties who have received copies, or rights, from you under
-this License will not have their licenses terminated so long as such
-parties remain in full compliance.
-
-  5. You are not required to accept this License, since you have not
-signed it.  However, nothing else grants you permission to modify or
-distribute the Program or its derivative works.  These actions are
-prohibited by law if you do not accept this License.  Therefore, by
-modifying or distributing the Program (or any work based on the
-Program), you indicate your acceptance of this License to do so, and
-all its terms and conditions for copying, distributing or modifying
-the Program or works based on it.
-
-  6. Each time you redistribute the Program (or any work based on the
-Program), the recipient automatically receives a license from the
-original licensor to copy, distribute or modify the Program subject to
-these terms and conditions.  You may not impose any further
-restrictions on the recipients' exercise of the rights granted herein.
-You are not responsible for enforcing compliance by third parties to
-this License.
-
-  7. If, as a consequence of a court judgment or allegation of patent
-infringement or for any other reason (not limited to patent issues),
-conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License.  If you cannot
-distribute so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you
-may not distribute the Program at all.  For example, if a patent
-license would not permit royalty-free redistribution of the Program by
-all those who receive copies directly or indirectly through you, then
-the only way you could satisfy both it and this License would be to
-refrain entirely from distribution of the Program.
-
-If any portion of this section is held invalid or unenforceable under
-any particular circumstance, the balance of the section is intended to
-apply and the section as a whole is intended to apply in other
-circumstances.
-
-It is not the purpose of this section to induce you to infringe any
-patents or other property right claims or to contest validity of any
-such claims; this section has the sole purpose of protecting the
-integrity of the free software distribution system, which is
-implemented by public license practices.  Many people have made
-generous contributions to the wide range of software distributed
-through that system in reliance on consistent application of that
-system; it is up to the author/donor to decide if he or she is willing
-to distribute software through any other system and a licensee cannot
-impose that choice.
-
-This section is intended to make thoroughly clear what is believed to
-be a consequence of the rest of this License.
-
-  8. If the distribution and/or use of the Program is restricted in
-certain countries either by patents or by copyrighted interfaces, the
-original copyright holder who places the Program under this License
-may add an explicit geographical distribution limitation excluding
-those countries, so that distribution is permitted only in or among
-countries not thus excluded.  In such case, this License incorporates
-the limitation as if written in the body of this License.
-
-  9. The Free Software Foundation may publish revised and/or new versions
-of the General Public License from time to time.  Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
-Each version is given a distinguishing version number.  If the Program
-specifies a version number of this License which applies to it and "any
-later version", you have the option of following the terms and conditions
-either of that version or of any later version published by the Free
-Software Foundation.  If the Program does not specify a version number of
-this License, you may choose any version ever published by the Free Software
-Foundation.
-
-  10. If you wish to incorporate parts of the Program into other free
-programs whose distribution conditions are different, write to the author
-to ask for permission.  For software which is copyrighted by the Free
-Software Foundation, write to the Free Software Foundation; we sometimes
-make exceptions for this.  Our decision will be guided by the two goals
-of preserving the free status of all derivatives of our free software and
-of promoting the sharing and reuse of software generally.
-
-                            NO WARRANTY
-
-  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
-FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
-OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
-PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
-OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
-TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
-PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
-REPAIR OR CORRECTION.
-
-  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
-REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
-INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
-OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
-TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
-YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
-PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
-POSSIBILITY OF SUCH DAMAGES.
-
-                     END OF TERMS AND CONDITIONS
-
-            How to Apply These Terms to Your New Programs
-
-  If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
-  To do so, attach the following notices to the program.  It is safest
-to attach them to the start of each source file to most effectively
-convey the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
-    <one line to give the program's name and a brief idea of what it does.>
-    Copyright (C) <year>  <name of author>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with this program; if not, see <http://www.gnu.org/licenses/>.
-
-Also add information on how to contact you by electronic and paper mail.
-
-If the program is interactive, make it output a short notice like this
-when it starts in an interactive mode:
-
-    Gnomovision version 69, Copyright (C) year name of author
-    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
-    This is free software, and you are welcome to redistribute it
-    under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License.  Of course, the commands you use may
-be called something other than `show w' and `show c'; they could even be
-mouse-clicks or menu items--whatever suits your program.
-
-You should also get your employer (if you work as a programmer) or your
-school, if any, to sign a "copyright disclaimer" for the program, if
-necessary.  Here is a sample; alter the names:
-
-  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
-  `Gnomovision' (which makes passes at compilers) written by James Hacker.
-
-  <signature of Ty Coon>, 1 April 1989
-  Ty Coon, President of Vice
-
-This General Public License does not permit incorporating your program into
-proprietary programs.  If your program is a subroutine library, you may
-consider it more useful to permit linking proprietary applications with the
-library.  If this is what you want to do, use the GNU Library General
-Public License instead of this License.

--- a/extra/divand/DESCRIPTION	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-Name: divand
-Version: 1.1.1
-Date: 2014-02-12
-Author: Alexander Barth <barth.alexander@gmail.com>
-Maintainer: Alexander Barth <barth.alexander@gmail.com>
-Title: divand
-Description: divand performs an n-dimensional variational analysis (interpolation) of arbitrarily located observations.
-Depends: octave (>= 3.4.0), mapping
-License: GPLv2+
-Url: http://modb.oce.ulg.ac.be/mediawiki/index.php/divand, http://www.geosci-model-dev.net/7/225/2014/gmd-7-225-2014.html

--- a/extra/divand/INDEX	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,123 +0,0 @@
-divand >> divand
-High-level functions
-  divand
-  divand_metric
-  divand_kernel
-Functions for advanced usage
-  divand_addc
-  divand_background_components
-  divand_background
-  divand_constr_advec
-  divand_diagnose
-  divand_eof_contraint
-  divand_error
-  divand_factorize
-  divand_laplacian
-  divand_obs
-  divand_operators
-  divand_orthogonalize_modes
-  divand_pc_none
-  divand_pc_michol
-  divand_pc_sqrtiB
-  divand_rms
-  divand_solve_eof
-  divand_solve
-Utility functions
-  cat_cell_array
-  colordisp
-  colormsg
-  conjugategradient
-  interp_regular
-  localize_regular_grid
-  localize_separable_grid
-  mtimescorr
-  statevector_init
-  statevector_pack
-  statevector_unpack
-Special matrices
-  @CatBlockMat/CatBlockMat
-  @CatBlockMat/append
-  @CatBlockMat/ctranspose
-  @CatBlockMat/double
-  @CatBlockMat/end
-  @CatBlockMat/full
-  @CatBlockMat/mldivide
-  @CatBlockMat/mtimes
-  @CatBlockMat/size
-  @CatBlockMat/subsref
-  @CovarFun/CovarFun
-  @CovarFun/mldivide
-  @CovarIS/CovarIS
-  @CovarIS/diag
-  @CovarIS/factorize
-  @CovarIS/full
-  @CovarIS/inv
-  @CovarIS/mldivide
-  @CovarIS/mtimes
-  @CovarIS/subsref
-  @CovarLanczos/CovarLanczos
-  @CovarLanczos/diag
-  @CovarLanczos/mtimes
-  @CovarParam/CovarParam
-  @CovarParam/diag
-  @CovarParam/full
-  @CovarParam/isscalar
-  @CovarParam/isvector
-  @CovarParam/mldivide
-  @CovarParam/mtimescorr
-  @CovarParam/mtimes
-  @CovarParam/size
-  @CovarParam/subsref
-  @CovarSMW/CovarSMW
-  @CovarSMW/diag
-  @CovarSMW/full
-  @CovarSMW/mldivide
-  @CovarSMW/mtimes
-  @CovarSMW/size
-  @DiagBlockMat/append
-  @DiagBlockMat/DiagBlockMat
-  @DiagBlockMat/diag
-  @DiagBlockMat/full
-  @DiagBlockMat/inv
-  @DiagBlockMat/mldivide
-  @DiagBlockMat/mtimes
-  @DiagBlockMat/size
-  @MatFun/ctranspose
-  @MatFun/double
-  @MatFun/end
-  @MatFun/full
-  @MatFun/MatFun
-  @MatFun/mtimes
-  @MatFun/size
-Sparse matrix operations
-  sparse_diag
-  sparse_diff
-  sparse_gradient
-  sparse_interp
-  sparse_pack
-  sparse_shift
-  sparse_stagger
-  sparse_trim
-Test script
-  test_1dvar
-  test_2dvar_adv
-  test_2dvar_check_correrr
-  test_2dvar_check
-  test_2dvar_constrain
-  test_2dvar_cyclic
-  test_2dvar_eof_check
-  test_2dvar_lenxy
-  test_2dvar
-  test_2dvar_pcg
-  test_2dvar_rellen
-  test_3dvar_large_stacked
-  test_3dvar
-  test_4dvar
-  test_divand
-  test_interp_1d
-  test_interp_2d
-  test_interp_regular
-  test_sparse_diff
-Example
-  divand_simple_example
-  divand_realistic_example
\ No newline at end of file

--- a/extra/divand/NEWS	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,15 +0,0 @@
-Summary of important user-visible changes for divand 1.1.1:
-------------------------------------------------------------
-
- ** Improved documentation
-
-Summary of important user-visible changes for divand 1.1:
-------------------------------------------------------------
-
- ** new preconditioner: divand_pc_michol and divand_pc_sqrtiB 
-    for iterative method
-
-Summary of important user-visible changes for divand 1.0:
-------------------------------------------------------------
-
- ** Initial release

--- a/extra/divand/doc/README	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1 +0,0 @@
-/home/abarth/matlab/diva-nd/README
\ No newline at end of file

Binary file extra/divand/doc/gmd-7-225-2014.pdf has changed

--- a/extra/divand/inst/@CatBlockMat/CatBlockMat.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,49 +0,0 @@
-% Create a matrix which represents the concatenation of a series of matrices.
-% 
-% C = CatBlockMat(dim,X1,X2,...)
-%
-% Create a matrix C which represents the concatenation of a series of matrices
-% X1,X2,... over dimension dim. The matrices X1,X2,... can be special matrix
-% objects. The resulting matrix C behaves as a regular matrix.
-
-function retval = CatBlockMat(dim,varargin)
-
-self.d = 1:2; % only for matrices now
-self.dim = dim;
-self.B = varargin;
-self.N = length(varargin);
-self.i = zeros(self.N+1,1);
-
-self.sz = size(self.B{1});
-self.sz(dim) = 0;
-
-
-for i=1:self.N
-  sz = size(self.B{i});
-  
-    if any(self.sz(self.d ~= dim) ~= sz(self.d ~= dim))
-        error('arguments dimensions are not consistent.');
-    end
-
-  self.i(i+1) = self.i(i) + sz(dim);
-end
-
-self.sz(dim) = self.i(end);
-
-retval = class(self,'CatBlockMat');
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/append.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,35 +0,0 @@
-% Append a matrix to a CatBlockMat matrix.
-%
-% A = append(B,M)
-%
-% Append the matrix M to a CatBlockMat matrix B.
-
-
-function self = append(self,M)
-
-sz = size(M);
-
-if any(self.sz(self.d ~= self.dim) ~= sz(self.d ~= self.dim))
-    error('arguments dimensions are not consistent.');
-end
-
-self.B{end+1} = M;
-self.N = self.N+1;
-
-
-self.i(end+1) = self.i(end) + sz(self.dim);
-self.sz(self.dim) = self.i(end);
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/ctranspose.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-% Conjugate transpose of a CatBlockMat matrix.
-%
-% T = ctranspose(C)
-%
-% Returns the conjugate transpose of matrix C.
-
-
-function T = ctranspose(self)
-
-B = cell(1,self.N);
-
-for l=1:self.N
-    B{l} = self.B{l}';
-end
-
-d = 3-self.dim;
-T = CatBlockMat(d,B{:});
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/double.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,29 +0,0 @@
-% Convert a CatBlockMat matrix to a regular matrix.
-%
-% F = double(C)
-%
-% Convert the CatBlockMat C matrix to a regular matrix F.
-
-function F = double(self)
-
-b = {};
-for l=1:self.N
-    b{l} = double(self.B{l});
-end
-
-F = cat(self.dim,b{:});
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/end.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Last index of a CatBlockMat matrix.
-%
-% e = end(C,k)
-%
-% Return the last index of the CatBlockMat C along dimension k.
-
-
-function e = end(self,k,n)
-
-e = self.sz(k);
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a CatBlockMat matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the CatBlockMat C matrix to a regular matrix F.
-
-function F = full(self)
-
-F = cat(self.dim,self.B{:});
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of CatBlockMat A and another
-% matrix B.
-
-
-function p = mldivide(self,b)
-
-if self.sz(2) ~= size(b,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-p = full(self) \ b;
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the CatBlockMat A and another matrix B.
-
-function p = mtimes(A,B)
-
-if size(A,2) ~= size(B,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-if isa(A,'CatBlockMat')
-    
-    if A.dim == 1
-        p = zeros(size(A,1),size(B,2));
-        
-        for l=1:A.N
-            i = A.i(l)+1:A.i(l+1);
-            p(i,:) = A.B{l} * B;
-        end
-    else
-        for l=1:A.N
-            i = (A.i(l)+1) :A.i(l+1);
-            
-            % must use subsref explicetly
-            % subsref(B,S) = B(i,:)
-            % http://www.mathworks.com/support/solutions/en/data/1-18J9B/
-            
-            S.type = '()';
-            S.subs = {i, ':'};
-            tmp = A.B{l} * subsref(B,S);
-            
-            if l==1
-                p = tmp;
-            else
-                p = p + tmp;
-            end
-        end        
-    end
-    
-elseif isa(B,'CatBlockMat') && B.dim == 2
-    p = zeros(size(A,1),size(B,2));
-    
-    for l=1:B.N
-        j = B.i(l)+1:B.i(l+1);
-        p(:,j) = A * B.B{l};
-    end    
-else
-    warning('divand:fullmat','use full matrices for product');
-    p = full(A) * full(B);
-end
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/size.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Size of a CatBlockMat matrix.
-%
-% sz = size(C)
-% sz = size(C,dim)
-%
-% Return the size of the CatBlockMat C matrix. If dim is specified, then only 
-% the size along dimension dim is returned.
-
-function sz = size(self,dim);
-
-sz = self.sz;
-
-if nargin == 2
-    sz = sz(dim);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CatBlockMat/subsref.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-% Subreference of a CatBlockMat matrix.
-%
-% B = subsref(C,idx)
-%
-% Perform the subscripted element selection operation according to
-% the subscript specified by idx.
-%
-% See also: subsref.
-
-function B = subsref(self,S)
-
-if strcmp(S.type,'()')           
-    i = S.subs{self.dim};
-    l = find(self.i+1 == i(1));
-    
-    itest = self.i(l)+1:self.i(l+1);
-    
-    if any(i ~= itest) || ~all(strcmp(':',S.subs(self.d ~= self.dim)))
-        error('only whole blocks can be referenced');
-    end
-    
-    B = self.B{l};    
-else
-    error('Attempt to reference field of non-structure array');
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarFun/CovarFun.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,60 +0,0 @@
-% Create a covariance matrix based on a function handler.
-%
-% CF = = CovarFun(n,fun)
-% CF = = CovarFun(n,fun,...)
-%
-% Create the covariance matrix CF of size n-by-n based on the function handler 
-% fun. Optional key word parameters are:
-%  'isvec': 1 if function is vectorized and 0 (default) if not.
-%  'pc': function hander of the preconditioner for pcg
-%  'M1','M2','matit','tol': parameters for pcg.
-%
-% See also: pcg
-
-
-function retval = CovarFun(n,fun,varargin)
-
-isvec = 0;
-self.pc = [];
-self.M1 = [];
-self.M2 = [];
-self.tol = 1e-6;
-self.maxit = 100;
-
-for i=1:2:length(varargin)
-    if strcmp(varargin{i},'isvec')
-        isvec = varargin{i+1};
-    elseif strcmp(varargin{i},'pc')
-        self.pc = varargin{i+1};
-    elseif strcmp(varargin{i},'M1')
-        self.M1 = varargin{i+1};
-    elseif strcmp(varargin{i},'M2')
-        self.M2 = varargin{i+1};
-    elseif strcmp(varargin{i},'maxit')
-        self.maxit = varargin{i+1};
-    elseif strcmp(varargin{i},'tol')
-        self.tol = varargin{i+1};
-    else
-        error('CovarFun:param','unknown param');
-    end
-end
-
-self.fun = fun;
-
-retval = class(self,'CovarFun',MatFun([n,n],fun,fun,isvec));
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarFun/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,48 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of CovarParam A and another
-% matrix B using pcg.
-
-function p = mldivide(self,b)
-
-if size(self,2) ~= size(b,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-if isempty(self.M1)
-    args = {self.pc};
-else
-    args = {self.M1,self.M2};
-end
-
-for i= 1:size(b,2)
-    
-    [p(:,i),flag,relres,iter] = pcg(self.fun, b(:,i), self.tol,...
-        self.maxit,args{:});    
-    
-    %whos b
-    self.pc
-    fprintf('iter %d %g \n',iter,relres);
-    if (flag ~= 0)
-        error('CovarFun:pcg', ['Preconditioned conjugate gradients method'...
-            ' did not converge %d %g %g'],flag,relres,iter);
-    end
-    
-end
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/CovarIS.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Covariance matrix with a sparse inverse matrix.
-% 
-% C = CovarIS(IS)
-%
-% Create the covariance matrix C whose inverse is the sparse matrix IS.
-% To accelerate matrix product C*x, IS can be factorized. 
-
-function retval = CovarIS(IS);
-
-self.IS = IS;
-self.f = 0;
-self.RP = [];
-self.q = [];
-self.PP = [];
-
-retval = class(self,'CovarIS');
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/diag.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,45 +0,0 @@
-% Diagonal elements.
-% 
-% d = diag(C)
-%
-% Return the diagonal elements of the CovarIS matrix C. 
-
-function d = diag(self)
-
-N = size(self.IS,1);
-d = zeros(N,1);
-
-
-if self.f  
-  t0 = cputime;
-
-  for i=1:N
-    % show progress every 5 seconds
-    if (t0 + 5 < cputime)
-      t0 = cputime;
-      fprintf('%d/%d\n',i,N);
-    end
-    
-    d(i) = self.PP(i,:) * (self.RP \ (self.RP' \ self.PP(i,:)'));
-  end
-else
-  disp('factorize CovarIS');
-  self_fact = factorize(self);
-  d = diag(self_fact);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/factorize.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-% Factorize matrix.
-% 
-% self = factorize(self)
-%
-% Factorize matrix self using the Cholesky decomposition.
-
-function self = factorize(self)
-
-
-if issparse(self.IS)
-  [self.RP, self.q, self.PP] = chol (self.IS);
-   
-  if self.q ~= 0
-    error('factoziation failed (matrix is not positive definite)');
-  end
-
-  
-  %length(find(self.PP))/numel(self.PP)
-  %length(find(self.RP))/numel(self.RP)
-
-else
-  self.RP = chol(self.IS);
-  self.PP = speye(size(self.IS));
-end
-
-self.f = 1;
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a CovarIS matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the CovarIS C matrix to a regular matrix F.
-
-function M = full(self)
-
-M = inv(self.IS);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/inv.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-% Inverse of a CovarIS matrix.
-%
-% M = inv(C)
-%
-% Return the inverse of the CovarIS matrix C.
-
-function M = inv(self)
-
-M = self.IS;
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of CovarIS A and another
-% matrix B.
-
-function q = mldivide(self,b)
-
-q = self.IS * b;
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,37 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the CovarIS matrix A and another matrix B.
-
-function p = mtimes(self,b)
-
-%keyboard
-if self.f
-  %'fact'
-  p = self.PP * (self.RP \ (self.RP' \ (self.PP' * b)));
-else  
-%  'direct'
-%  tic
-%  A = (self.IS + self.IS')/2;
-%  p = A \ b;
-%  toc
-%  tic
-  p = self.IS \ b;  
-%  toc
-%  'end'
-end
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarIS/subsref.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,67 +0,0 @@
-% Subreference of a CovarIS matrix.
-%
-% B = subsref(C,idx)
-%
-% Perform the subscripted element selection operation according to
-% the subscript specified by idx.
-%
-% See also: subsref.
-
-function x = subsref(self,idx)
-
-assert(length(idx) == 1)
-assert(strcmp(idx.type,'()'))
-N = size(self.IS,1);
-
-i = idx.subs{1};
-
-if strcmp(i,':')
-    i = 1:N;
-end
-
-j = idx.subs{2};
-
-if strcmp(j,':')
-    j = 1:N;
-end
-
-x = zeros(length(i),length(j));
-
-if self.f  
-  t0 = cputime;
-
-  k = 0;
-  for jj=1:length(j)
-    for ii=1:length(i)
-      
-      % show progress every 5 seconds
-      if (t0 + 5 < cputime)
-        t0 = cputime;
-        fprintf('%d/%d\n',k,numel(x));
-      end
-    
-      x(ii,jj) = self.PP(i(ii),:) * (self.RP \ (self.RP' \ self.PP(j(jj),:)'));
-      k = k+1;
-    end
-  end
-else
-  disp('factorize CovarIS');
-  self_fact = factorize(self);
-  x = subsref(self_fact,idx);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarLanczos/CovarLanczos.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,28 +0,0 @@
-% Covariance matrix reconstructed from Lanczos vectors.
-% 
-% A = CovarLanczos(Q,T)
-%
-% Return a covariance matrix, factorized into A = Q * T * Q'
-% where Q'*Q = I and T is a tridiagonal matrix.
-
-function retval = CovarLanczos(Q,T)
-
-self.Q = Q;
-self.T = T;
-
-retval = class(self,'CovarLanczos');
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarLanczos/diag.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-% Diagonal elements.
-% 
-% d = diag(C)
-%
-% Return the diagonal elements of the CovarLanczos matrix C. 
-
-function d = diag(self)
-
-N = size(self.Q,1);
-d = zeros(N,1);
-
-for i=1:N
-  d(i) = self.Q(i,:) * (self.T \ self.Q(i,:)');
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarLanczos/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the CovarLanczos matrix A and another matrix
-% B.
-
-function p = mtimes(self,b)
-
-p = self.Q * (self.T \ (self.Q' * b));
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/CovarParam.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-% Parametric covariance matrix with a guassian kernel.
-%
-% C = CovarParam({x,y,...},variance,len)
-%
-% Creates a parametric covariance matrix C with a guassian 
-% kernel, a correlation length of len and a variance of 
-% variance. The first parameter is a cell array with the coordinates.
-
-function retval = CovarParam(xi,variance,len,varargin)
-
-xi = squeeze(cat_cell_array(xi));
-m = size(xi,1);
-
-
-% fun = @(xi1,xi2) exp(- sum( (xi1 - xi2).^2 , 2) / len^2);
-% 
-% tic
-% for j=1:m    
-%     for i=1:m    
-%         R2(i,j) = fun(xi(i,:),xi(j,:));
-%     end
-% end
-% toc
-% 
-% maxdiff(R,R2);
-
-self.tol = 1e-6;
-self.maxit = 100;
-self.var_add = 1;
-
-for i=1:2:length(varargin)
-    if strcmp(varargin{i},'tolerance')
-        self.tol = varargin{i+1};
-    elseif strcmp(varargin{i},'maxit')
-        self.maxit = varargin{i+1};
-    elseif strcmp(varargin{i},'fraccorr')
-        self.var_add = varargin{i+1};
-    end
-end
-    
-self.m = m;
-self.xi = xi;
-self.len = len;
-
-if isscalar(variance)
-  variance = ones(self.m,1) * variance;
-end
-  
-self.variance = variance;
-self.std = sqrt(variance);
-self.S = sparse_diag(self.std);
-
-retval = class(self,'CovarParam');
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/diag.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Diagonal elements.
-% 
-% d = diag(C)
-%
-% Return the diagonal elements of the CovarParam matrix C. 
-
-function d = diag(self);
-
-d = self.variance;
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Convert a CovarParam matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the CovarParam C matrix to a regular matrix F.
-
-
-function F = full(self)
-
-F = self * eye(size(self));
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/isscalar.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Return true if argument is a scalar.
-%
-% b = isscalar(X)
-%
-% Return true if X is a scalar.     
-
-function b = isscalar(self)
-
-b = isequal(size(self),[1 1]);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/isvector.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Return true if argument is a vector.
-%
-% b = isscalar(X)
-%
-% Return true if X is a vector. A vector is a N-by-1 or 1-by-N matrix.
-
-function b = isvector(self)
-
-sz = size(self);
-b = length(sz) == 2 && (sz(1) == 1 || sz(2) == 1);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of CovarParam matrix A and 
-% another matrix B.
-
-function q = mldivide(self,b)
-
-q = zeros(size(b));
-b = self.S \ b; 
-
-for i=1:size(b,2)
-  M = [];
-  %M = @(x) ifft(x);
-  
-  if 1
-    [q(:,i),flag,relres,iter] = pcg(@(x) mtimescorr(self,x), b(:,i),self.tol,self.maxit,M);
-  else
-    
-    [q(:,i),flag,relres,iter] = pcg(@(x) fft(mtimescorr(self,ifft(x))), fft(b(:,i)),self.tol,self.maxit,M);
-    q(:,i) = real(ifft(q(:,i)));
-  end
-
-  
-    if flag > 0
-        relres,iter
-        warning('divand:pcgFlag','pcg flag %d',flag);
-    end
-end
-
-q = self.S \ q; 
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,39 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the CovarParam matrix A and another matrix 
-% B.
-
-function p = mtimes(self,b)
-
-m = self.m;
-p = zeros(size(b));
-
-% for j=1:m
-%     for i=1:m
-%         d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 );
-%         p(i) = p(i) + exp(-d2 / self.len^2) * b(j);
-%     end
-% end
-%
-
-b = self.S*b;
-p = mtimescorr(self,b);
-p = self.S*p;
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/mtimescorr.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,40 +0,0 @@
-% Matrix product with correlation matrix.
-%
-% p = mtimescorr(A,B)
-%
-% Return the matrix product between the CovarParam matrix A (with unit variance)
-% and another matrix B.
-
-
-function p = mtimescorr(self,b)
-
-m = self.m;
-p = zeros(size(b));
-
-% for j=1:m
-%     for i=1:m
-%         d2 = sum( (self.xi(i,:) - self.xi(j,:)).^2 );
-%         p(i) = p(i) + exp(-d2 / self.len^2) * b(j);
-%     end
-% end
-%
-
-p = mtimescorr(self.xi,self.len,b);
-%p = mtimescorr_f(self.xi,self.len,b);
-p = self.var_add * p + (1-self.var_add) * b;
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/size.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Size of a CovarParam matrix.
-%
-% sz = size(C)
-% sz = size(C,dim)
-%
-% Return the size of the CovarParam C matrix. If dim is specified, then only 
-% the size along dimension dim is returned.
-
-function sz = size(self,dim);
-
-sz = [self.m, self.m];
-
-if nargin == 2
-    sz = sz(dim);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarParam/subsref.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,35 +0,0 @@
-% Subreference of a CovarParam matrix.
-%
-% B = subsref(C,idx)
-%
-% Perform the subscripted element selection operation according to
-% the subscript specified by idx.
-%
-% See also: subsref.
-
-function C = subsref(self,idx)
-
-assert(strcmp(idx.type,'()'))
-
-i = idx.subs{1};
-j = idx.subs{2};
-
-assert(isequal(i,j));
-
-
-C = CovarParam(self.xi(i,:),self.variance(i),self.len);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/CovarSMW.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,46 +0,0 @@
-% Covariance matrix than can be inverted using the Sherman-Morrison-Woodbury 
-% forumla.
-%
-% CSMW = CovarSMW(C,B)
-%
-% Create the covariance matrix CSMW = C +  B*B' than can be inverted efficiently
-% using the Sherman-Morrison-Woodbury formula, if size(B,2) is much smaller than
-% size(C,1):
-%
-% inv(C +  B*B') = inv(C) -  inv(C)*B*inv(B'*inv(C)*B +  I)*B'*inv(C) 
-%
-% The  symetric matrix C should implement the methods: size, diag, mtimes, 
-% mldivde and the matrix B should implement the methods: size, tranpose, mtimes,
-% sum.
-
-function retval = CovarSMW(C,B)
-
-if size(C,1) ~= size(C,2)
-  error('C should be square matrix');
-end
-
-if size(C,1) ~= size(B,1)
-  error('size of C and B are not compatible');
-end
-
-self.C = C;
-self.B = B;
-self.D = inv(B' * (C \ B) +  eye(size(B,2)));
-
-retval = class(self,'CovarSMW');
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/diag.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Diagonal elements.
-% 
-% d = diag(C)
-%
-% Return the diagonal elements of the CovarSWM matrix C. 
-
-function d = diag(self);
-
-d  = diag(self.C) +  sum(self.B.^2,2);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a CovarSWM matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the CovarSWM C matrix to a regular matrix F.
-
-function F = full(self)
-
-F  = self.C +  self.B * self.B';
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of CovarSWM A and another
-% matrix B.
-
-function Q = mldivide(self,A)
-
-Q = self.C \ A -  self.C \ (self.B * (self.D * (self.B' * (self.C \ A))));
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the CovarSWM matrix A and another matrix B.
-
-function P = mtimes(self,A)
-
-P  = self.C*A +  self.B *  (self.B'*A);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@CovarSMW/size.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Size of a CovarSWM matrix.
-%
-% sz = size(C)
-% sz = size(C,dim)
-%
-% Return the size of the CovarSWM C matrix. If dim is specified, then only 
-% the size along dimension dim is returned.
-
-function sz = size(self,dim);
-
-sz = size(self.C);
-
-if nargin == 2
-    sz = sz(dim);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/DiagBlockMat.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,44 +0,0 @@
-% Block diagonal matrix.
-%
-% DB = DiagBlockMat(A1,A2,...,An)
-%
-% Create the block diagonal matrix DB from the matrices A1, A2,...,An
-%
-%      /  A1   0    0    ... 0 \
-%      |  0    A2   0        : |
-% DB = |  0    0    A3         |
-%      |  :              .   0 |
-%      \  0  ...         0  An /  
-
-function retval = DiagBlockMat(varargin)
-
-self.B = varargin;
-self.N = length(varargin);
-self.i = zeros(self.N+1,1);
-self.j = zeros(self.N+1,1);
-
-for i=1:self.N
-  sz = size(self.B{i});
-  self.i(i+1) = self.i(i) + sz(1);
-  self.j(i+1) = self.j(i) + sz(2);  
-end
-
-self.sz = [self.i(end), self.j(end)];
-
-retval = class(self,'DiagBlockMat');
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/append.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,32 +0,0 @@
-% Append a matrix to a DiagBlockMat matrix.
-%
-% A = append(B,M)
-%
-% Append the matrix M to a DiagBlockMat matrix B along its diagonal.
-
-function self = append(self,M)
-
-self.B{end+1} = M;
-self.N = self.N+1;
-
-sz = size(M);
-
-self.i(end+1) = self.i(end) + sz(1);
-self.j(end+1) = self.j(end) + sz(2);  
-
-self.sz = [self.i(end), self.j(end)];
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/diag.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-% Diagonal elements.
-% 
-% d = diag(C)
-%
-% Return the diagonal elements of the matrix C. 
-
-function d = diag(self)
-
-t = cell(1,self.N);
-
-for l=1:self.N
-    t{l} = diag(self.B{l});
-end
-
-d = cat(1,t{:});
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a DiagBlockMat matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the DiagBlockMat C matrix to a regular matrix F.
-
-function F = full(self)
-
-F = blkdiag(self.B{:});
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/inv.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,28 +0,0 @@
-% Inverse of a DiagBlockMat matrix.
-%
-% M = inv(C)
-%
-% Return the inverse of the DiagBlockMat matrix C.
-
-function iA = inv(A)
-
-iA = A;
-
-for l=1:A.N
-    iA.B{l} = inv(A.B{l});
-end
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/mldivide.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,43 +0,0 @@
-% Matrix left division A \ B.
-%
-% p = mldivide(A,B)
-%
-% Return the matrix product between the inverse of DiagBlockMat A and another
-% matrix B
-
-function p = mldivide(self,b)
-
-if size(self,2) ~= size(b,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-
-if isa(b,'CatBlockMat')
-    B = {};
-    for l=1:self.N
-        B{l} = self.B{l} \ b(self.j(l)+1:self.j(l+1),:);
-    end
-    
-   p = CatBlockMat(1,B{:});
-else
-    
-    p = zeros(size(self,1),size(b,2));
-    
-    for l=1:self.N
-        p(self.i(l)+1:self.i(l+1),:) = self.B{l} \ b(self.j(l)+1:self.j(l+1),:);
-    end
-end
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the DiagBlockMat matrix A and another matrix
-% B.
-
-function p = mtimes(self,b)
-
-if size(self,2) ~= size(b,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-p = zeros(size(self,1),size(b,2));
-
-for l=1:self.N
-  p(self.i(l)+1:self.i(l+1),:) = self.B{l} * b(self.j(l)+1:self.j(l+1),:);
-end
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@DiagBlockMat/size.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Size of a DiagBlockMat matrix.
-%
-% sz = size(C)
-% sz = size(C,dim)
-%
-% Return the size of the DiagBlockMat C matrix. If dim is specified, then only 
-% the size along dimension dim is returned.
-
-function sz = size(self,dim);
-
-sz = self.sz;
-
-if nargin == 2
-    sz = sz(dim);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/MatFun.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,39 +0,0 @@
-% Matrix operator object based on a function handel.
-%
-% MF = MatFun(sz,fun,funt,isvec)
-%
-% Create a operator of size sz based on the linear function fun and its adjoint 
-% funt. For a matrix X of approriate size:
-%
-% MF*X is fun(X) and MF'*X is funt(X)
-% 
-% If isvec is 0 (default), then the function must be applied on the columns of X
-% individually, otherwise fun and funt are assumed to be "vectorized". 
-
-function retval = MatFun(sz,fun,funt,isvec)
-
-if nargin == 3
-  isvec = 0;
-end
-
-self.sz = sz;
-self.fun = fun;
-self.funt = funt;
-self.isvec = isvec;
-
-retval = class(self,'MatFun');
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/ctranspose.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-% Conjugate transpose of a MatFun matrix.
-%
-% T = ctranspose(C)
-%
-% Returns the conjugate transpose of matrix C.
-
-function T = ctranspose(self)
-
-T = MatFun(self.sz([2 1]),self.funt,self.fun);
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/double.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a MatFun matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the MatFun C matrix to a regular matrix F.
-
-function d = double(self)
-
-d = full(self);
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/end.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-% Last index of a MatFun matrix.
-%
-% e = end(C,k)
-%
-% Return the last index of the MatFun matrix C along dimension k.
-
-function e = end(self,k,n)
-
-e = self.sz(k);
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/full.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-% Convert a MatFun matrix to a regular matrix.
-%
-% F = full(C)
-%
-% Convert the MatFun C matrix to a regular matrix F.
-
-function M = full(self);
-
-M = self * eye(self.sz(2));
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/mtimes.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,40 +0,0 @@
-% Matrix product A*B.
-%
-% p = mtimes(A,B)
-%
-% Return the matrix product between the MatFun matrix A and another matrix B.
-
-function p = mtimes(A,B)
-
-if size(A,2) ~= size(B,1)
-    error('Inner matrix dimensions must agree.');
-end
-
-if isa(A,'MatFun')  
-  if A.isvec
-    p = A.fun(B);
-  else    
-    p = zeros(size(A,1),size(B,2));
-  
-    for l=1:size(B,2)
-      p(:,l) = A.fun(B(:,l));
-    end
-  end
-else
-  error('not implemented');
-end
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/@MatFun/size.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Size of a MatFun matrix.
-%
-% sz = size(C)
-% sz = size(C,dim)
-%
-% Return the size of the MatFun matrix C. If dim is specified, then only 
-% the size along dimension dim is returned.
-
-function sz = size(self,dim);
-
-sz = self.sz;
-
-if nargin == 2
-    sz = sz(dim);
-end
-
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/cat_cell_array.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,56 +0,0 @@
-% Concatenate a cell array.
-%
-% x = cat_cell_array(c)
-%
-% Concatenate all arrays of the cell array c along the dimension length(c)+1
-
-
-function x = cat_cell_array(x)
-
-do_grid = 0;
-
-if iscell(x)
-  n = length(x);
-  
-  if do_grid
-    all_vectors = 1;
-  
-    for i=1:n
-      all_vectors = all_vectors & isvector(x{i});
-    end
-
-    if all_vectors
-      [x{:}] = ndgrid(x{:});
-    end
-  end
-
-  sz = size(x{1});
-  for i=2:n
-    if ~isequal(sz,size(x{i}))
-      %sz
-      %size(x{i})
-      error('all arguments in cell array should have the same size');
-    end
-  end
-    
-  x = cat(n+1,x{:});
-  end
-end
-
-  
-
-% Copyright (C) 2009 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation; either version 2 of the License, or
-% (at your option) any later version.
-%
-% This program is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program; If not, see <http://www.gnu.org/licenses/>.
-

--- a/extra/divand/inst/colordisp.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,27 +0,0 @@
-% Display a message in a color (followed by a newline).
-%
-% colordisp(msg,color)
-%
-% The same as disp but in color.
-%
-% See also
-%   colormsg
-
-function colordisp(msg,color)
-colormsg (msg,color)
-fprintf(1,'\n');
-
-% Copyright (C) 2004 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/colormsg.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,71 +0,0 @@
-% Display a message in a color (without a newline).
-%
-% colormsg(msg,color)
-%
-% The message msg in printed on the terminal the specified color. Acceptable 
-% values color color are:
-%   black  
-%   red    
-%   green  
-%   yellow 
-%   blue   
-%   magenta
-%   cyan   
-%   white  
-%
-% In Matlab the message is displayed in black (due to the lack of ANSI escape 
-% code support).
-%
-% See also
-%   colordisp
-
-function colormsg (msg,color)
-
-%getenv('TERM')
-%if strcmp(getenv('TERM'),'xterm') % && exist('puts') > 1
-% only in octave
-if exist('puts') > 1
-  esc = char(27);          
-
-  % ANSI escape codes
-  colors.black   = [esc, '[40m'];
-  colors.red     = [esc, '[41m'];
-  colors.green   = [esc, '[42m'];
-  colors.yellow  = [esc, '[43m'];
-  colors.blue    = [esc, '[44m'];
-  colors.magenta = [esc, '[45m'];
-  colors.cyan    = [esc, '[46m'];
-  colors.white   = [esc, '[47m'];
-
-  reset = [esc, '[0m'];
-  
-  c = getfield(colors,color);
-           
-  %oldpso = page_screen_output (0);
-
-  try
-    fprintf([c, msg, reset]);
-    %puts (reset); % reset terminal
-  catch
-    %page_screen_output (oldpso);
-  end
-else
-  fprintf(msg);
-end
-           
-
-
-% Copyright (C) 2004 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/conjugategradient.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,198 +0,0 @@
-% Solve a linear system with the preconditioned conjugated-gradient method.
-%
-% [x] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
-% [x,Q,T] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
-%
-% solves for x
-% A x = b
-% using the preconditioned conjugated-gradient method.
-% It provides also an approximation of A:
-% A \sim Q*T*Q'
-%
-% J(x) = 1/2 (x' b - x' A x)
-% ∇ J = b - A x
-% A x = b - ∇ J
-% b = ∇ J(0)
-%
-% the columns of Q are the Lanczos vectors
-
-% Alexander Barth 2010,2012-2013
-
-function [x,Q,T,diag] = conjugategradient(fun,b,varargin)
-
-n = length(b);
-
-% default parameters
-maxit = min(n,20);
-minit = 0;
-tol = 1e-6;
-pc = [];
-x0 = [];
-renorm = 0;
-
-prop = varargin;
-
-for i=1:2:length(prop)
-    if strcmp(prop{i},'minit')
-        minit = prop{i+1};
-    elseif strcmp(prop{i},'maxit')
-        maxit = prop{i+1};
-    elseif strcmp(prop{i},'tol')
-        tol = prop{i+1};
-    elseif strcmp(prop{i},'x0')
-        x0 = prop{i+1};
-    elseif strcmp(prop{i},'pc')
-        pc = prop{i+1};
-    elseif strcmp(prop{i},'renorm')
-        renorm = prop{i+1};
-    else
-        
-    end
-end
-
-
-if isempty(x0)
-    % random initial vector
-    x0 = randn(n,1);
-end
-
-
-if isempty(pc)
-    % no pre-conditioning
-    pc = @(x) x;
-elseif isnumeric(pc)
-    invM = pc;
-    pc = @(x)( invM *x);
-end
-
-
-tol2 = tol^2;
-
-delta = [];
-gamma = [];
-q = NaN*zeros(n,1);
-
-%M = inv(invM);
-%E = chol(M);
-
-
-% initial guess
-x = x0;
-
-% gradient at initial guess
-r = b - fun(x);
-
-% quick exit
-if r'*r < tol2
-    return
-end
-
-% apply preconditioner
-z = pc(r);
-
-% first search direction == gradient
-p = z;
-
-% compute: r' * inv(M) * z (we will need this product at several
-% occasions)
-
-zr_old = r'*z;
-
-% r_old: residual at previous iteration
-r_old = r;
-
-for k=1:maxit
-    if k <= n && nargout > 1
-        % keep at most n vectors
-        Q(:,k) = r/sqrt(zr_old);
-    end
-    
-    % compute A*p
-    Ap = fun(p);
-    %maxdiff(A*p,Ap)
-    
-    % how far do we need to go in direction p?
-    % alpha is determined by linesearch
-    
-    % alpha z'*r / (p' * A * p)
-    alpha(k) = zr_old / ( p' * Ap);
-    
-    % get new estimate of x
-    x = x + alpha(k)*p;
-    
-    % recompute gradient at new x. Could be done by
-    % r = b-fun(x);
-    % but this does require an new call to fun
-    r = r - alpha(k)*Ap;
-    
-    if renorm
-        r = r - Q(:,1:k) * Q(:,1:k)' * r ;
-    end
-    %maxdiff(r,b-fun(x))
-    
-    % apply pre-conditionner
-    z = pc(r);
-    
-    
-    zr_new = r'*z;
-        
-    if r'*r < tol2 && k >= minit
-        break
-    end
-    
-    %Fletcher-Reeves
-    beta(k+1) = zr_new / zr_old;
-    %Polak-Ribiere
-    %beta(k+1) = r'*(r-r_old) / zr_old;
-    %Hestenes-Stiefel
-    %beta(k+1) = r'*(r-r_old) / (p'*(r-r_old));
-    %beta(k+1) = r'*(r-r_old) / (r_old'*r_old);
-    
-    
-    % norm(p)
-    p = z + beta(k+1)*p;
-    zr_old = zr_new;
-    r_old = r;
-end
-
-
-%disp('alpha and beta')
-%figure,plot(beta(2:end))
-%rg(alpha)
-%rg(beta(2:end))
-
-if nargout > 1
-    kmax = size(Q,2);
-    
-    delta(1) = 1/alpha(1);
-    
-    %delta(1) - Q(:,1)'*invM*A*invM*Q(:,1)
-    
-    
-    for k=1:kmax-1
-        delta(k+1) = 1/alpha(k+1) + beta(k+1)/alpha(k);
-        gamma(k) = -sqrt(beta(k+1))/alpha(k);
-    end
-    
-    T = sparse([1:kmax   1:kmax-1 2:kmax  ],...
-        [1:kmax   2:kmax   1:kmax-1],...
-        [delta    gamma    gamma]);
-    
-    diag.iter = k;
-    diag.relres = sqrt(r'*r); 
-end
-
-% Copyright (C) 2004 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/divand.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,290 +0,0 @@
-% Compute a variational analysis of arbitrarily located observations.
-%
-% [fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
-% [fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...);
-% [fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
-%
-% Perform an n-dimensional variational analysis of the observations f located at
-% the coordinates x. The array fi represent the interpolated field at the grid
-% defined by the coordinates xi and the scales factors pmn.
-%
-% Input:
-%   mask: binary mask delimiting the domain. 1 is inside and 0 outside.
-%         For oceanographic application, this is the land-sea mask.
-%
-%   pmn: scale factor of the grid. pmn is a cell array with n elements. Every 
-%        element represents the scale factor of the corresponding dimension. Its
-%        inverse is the local resolution of the grid in a particular dimension.
-%
-%   xi: cell array with n elements. Every element represents a coordinate
-%   of the final grid on which the observations are interpolated
-%   x: cell array with n elements. Every element represents a coordinate of
-%   the observations
-%   f: value of the observations *minus* the background estimate (m-by-1 array).
-%     (see note)
-%   len: correlation length
-%   lambda: signal-to-noise ratio of observations (if lambda is a scalar).
-%     The larger this value is, the closer is the field fi to the
-%     observation.
-%     if lambda is a scalar:
-%        R = 1/lambda I, where R is the observation error covariance
-%        matrix),
-%     if lambda is a vector
-%     a vector (R = diag(lambda)) or a matrix or a CovarParam object (R = 
-%     lambda).
-%
-%
-% Optional input arguments specified as pairs of keyword and values:
-%  'velocity', vel: velocity of advection constraint. The default is 
-%        no-advection constraint
-%
-%  'alpha': alpha is vector of coefficients multiplying various terms in the 
-%        cost function. The first element multiplies the norm.
-%        The other i-th element of alpha multiplies the (i+1)-th derivative. 
-%        Per default, the highest derivative is m = ceil(1+n/2) where n is the 
-%        dimension of the problem.
-%
-%        The values of alpha is the (m+1)th row of the Pascal triangle:
-%           m=0         1
-%           m=1       1   1
-%           m=1     1   2   1     (n=1,2)
-%           m=2   1   3   3   1   (n=3,4)
-%           ...
-%
-%  'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or 
-%        off (0, default). If on, they will be returned as the last output 
-%        argument
-%
-%  'EOF', EOF: sub-space constraint. Orthogonal (EOF' WE^2 EOF = I) (units of 
-%        EOF: m^(-n/2))
-%
-%  'EOF_scaling', EOF_scaling: (dimensional)
-%
-%  'constraint': a structure with user specified constrain
-%
-%  'moddim': modulo for cyclic dimension (vector with n elements). 
-%      Zero is used for non-cyclic dimensions. Halo points should 
-%      not be included for cyclic dimensions. For example if the first dimension
-%      is cyclic, then the grid point corresponding to mask(1,j) should be 
-%      between mask(end,1) (left neighbor) and mask(2,j) (right neighbor)
-%
-%  'fracdim': fractional indices (n-by-m array). If this array is specified, 
-%      then x and xi are not used.
-%
-%  'inversion': direct solver ('chol' for Cholesky factorization) or a 
-%      interative solver ('pcg' for preconditioned conjugate gradient) can be 
-%      used.
-%
-%  'compPC': function that returns a preconditioner for the primal formulation 
-%      if inversion is set to 'pcg'. The function has the following arguments:
-%
-%            [M1,M2] = compPC(iB,H,R)
-%
-%     where iB is the inverse background error covariance, H the observation 
-%     operator and R the error covariance of the observation. The used 
-%     preconditioner M will be M = M1 * M2 (or just M = M1 if M2 is empty).
-%     Per default a modified incomplete Cholesky factorization will be used a 
-%     preconditioner.
-%
-%  Note: 'velocity' and 'constraint' may appear multiple times
-%
-% Output:
-%   fi: the analysed field
-%   err: error variance of the analysis field relative to the error variance of 
-%     the background
-%   s: structure
-%     s.iB: adimensional
-%     s.E: scaled EOF (dimensional)
-%
-% Note:
-%   If zero is not a valid first guess for your variable (as it is the case for 
-%   e.g. ocean temperature), you have to subtract the first guess from the 
-%   observations before calling divand and then add the first guess back in.
-%
-% Example:
-%   see divand_simple_example.m
-%
-
-
-function varargout =  divand(mask,pmn,xi,x,f,len,lambda,varargin)
-
-% default values
-
-velocity = {};
-EOF = [];
-diagnostics = 0;
-EOF_lambda = 0;
-primal = 1;
-factorize = 1;
-tol = 1e-6;
-maxit = 100;
-minit = 10;
-constraints = {};
-inversion = 'chol';
-moddim = [];
-fracindex = [];
-alpha = [];
-keepLanczosVectors = 0;
-compPC = @divand_pc_sqrtiB;
-
-prop = varargin;
-for i=1:length(prop)/2
-    if strcmp(prop{2*i-1},'velocity')
-        velocity{end+1}  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'EOF')
-        EOF  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'EOF_scaling')
-        EOF_lambda  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'primal')
-        primal  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'factorize')
-        factorize  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'inversion')
-        inversion  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'tol')
-        tol  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'maxit')
-        maxit  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'minit')
-        minit  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'diagnostics')
-        diagnostics  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'constraint')
-        constraints{end+1}  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'fracindex')
-        fracindex  = prop{2*i};
-    elseif strcmpi(prop{2*i-1},'moddim')
-        moddim  = prop{2*i};        
-    elseif strcmpi(prop{2*i-1},'alpha')
-        alpha  = prop{2*i};        
-    elseif strcmpi(prop{2*i-1},'keepLanczosVectors')
-        keepLanczosVectors  = prop{2*i};        
-    elseif strcmpi(prop{2*i-1},'compPC')
-        compPC  = prop{2*i};        
-    else
-      warning('divand:unknownProperty','unknown property "%s" is irgnored',prop{2*i-1});
-    end
-end
-
-% check inputs
-
-if ~any(mask(:))
-  error('no sea points in mask');
-end
-
-s = divand_background(mask,pmn,len,alpha,moddim);
-s.betap = 0;
-s.EOF_lambda = EOF_lambda;
-s.primal = primal;
-s.factorize = factorize;
-s.tol = tol;
-s.maxit = maxit;
-s.minit = minit;
-s.inversion = inversion;
-s.keepLanczosVectors = keepLanczosVectors;
-s.compPC = compPC;
-
-% remove non-finite elements from observations
-f = f(:);
-valid = isfinite(f);
-x = cat_cell_array(x);
-
-if ~all(valid)  
-  fprintf(1,'remove %d (out of %d) non-finite elements from observation vector\n',sum(~valid),numel(f));
-  x = reshape(x,[length(f) s.n]);
-  f = f(valid);
-  x = reshape(x(repmat(valid,[1 s.n])),[length(f) s.n]);
-  
-  if ~isempty(fracindex)
-    fracindex = fracindex(:,valid);
-  end
-
-  if isscalar(lambda)
-    % do nothing
-  elseif isvector(lambda)
-    lambda = lambda(valid);
-  elseif ismatrix(lambda)
-    lambda = lambda(valid,valid);
-  end
-end
-
-apply_EOF_contraint = ~(isempty(EOF) | all(EOF_lambda == 0));
-
-s.mode = 1;
-
-if ~apply_EOF_contraint
-    s.betap = 0;
-else
-    if s.mode==0
-        s.betap = max(EOF_lambda)/s.coeff;  % units m^(-n)
-    elseif s.mode==1
-        s.betap = max(max(EOF_lambda)-1,0)/s.coeff;
-    end
-end
-
-%assert(s.betap,0,1e-8)
-% increase contraint on total enegery to ensure system is still positive defined
-%s.betap
-s.iB = s.iB + s.betap * s.WE'*s.WE;
-
-% add observation constrain to cost function
-s = divand_addc(s,divand_obs(s,xi,x,f,lambda,fracindex));
-
-% add advection constraint to cost function
-for i=1:length(velocity)
-    %s = divand_advection(s,velocity);
-    s = divand_addc(s,divand_constr_advec(s,velocity{i}));
-end
-
-
-% add all additional constrains
-for i=1:length(constraints)
-    s = divand_addc(s,constraints{i});
-end
-
-
-if apply_EOF_contraint
-    s = divand_eof_contraint(s,EOF_lambda,EOF);
-end
-
-% factorize a posteori error covariance matrix
-% or compute preconditioner
-s = divand_factorize(s);
-
-if ~apply_EOF_contraint
-    [fi,s] = divand_solve(s,f);
-else
-    [fi,s] = divand_solve_eof(s,f);
-end
-
-varargout{1} = fi;
-
-if nargout-diagnostics >= 2
-    err = divand_error(s);
-    varargout{2} = err;
-end
-
-if diagnostics
-    s.B = CovarIS(s.iB);
-    [s.Jb,s.Jo,s.Jeof,s.J] = divand_diagnose(s,fi,f);
-    s.valid = valid;
-    varargout{nargout} = s;
-end
-
-
-% Copyright (C) 2008-2014 Alexander Barth <barth.alexander@gmail.com>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.
-
-% LocalWords:  fi divand pmn len diag CovarParam vel ceil moddim fracdim
\ No newline at end of file

--- a/extra/divand/inst/divand_addc.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-% Add a constraint to the cost function.
-%
-% s = divand_addc(s,constrain)
-%
-% Include in the structure s the specified constrain. 
-%
-% Input:
-%   s: structure created by divand_background
-%   constrain: The parameter constrain has the following fields: R (a covariance
-%     matrix), H (extraction operator) and yo (specified value for the 
-%     constrain).
-%
-% Output:
-%   s: structure to be used by divand_factorize
-
-function s = divand_addc(s,constrain)
-
-if isfield(s,'R')
-  s.R = append(s.R,constrain.R);
-  s.H = append(s.H,constrain.H);
-  s.yo = cat(1,s.yo,constrain.yo);
-else
-   s.H = CatBlockMat(1,constrain.H);
-   s.R = DiagBlockMat(constrain.R);
-   s.yo = constrain.yo; 
-end
-
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.

--- a/extra/divand/inst/divand_background.m	Mon Feb 17 10:02:10 2014 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,265 +0,0 @@
-% Form the inverse of the background error covariance matrix.
-%
-% s = divand_background(mask,pmn,Labs,alpha,moddim)
-%
-% Form the inverse of the background error covariance matrix with
-% finite-difference operators on a curvilinear grid
-%
-% Input:
-%   mask: binary mask delimiting the domain. 1 is inside and 0 outside.
-%         For oceanographic application, this is the land-sea mask.
-%
-%   pmn: scale factor of the grid.
-%
-%   Labs: correlation length
-%
-%   alpha: a dimensional coefficients for norm, gradient, laplacian,...
-%      alpha is usually [1 2 1]. 
-%       
-%
-% Output:
-%   s: stucture containing
-%   s.iB: inverse of the background error covariance
-%   s.L: spatial average correlation length
-%   s.n: number of dimenions
-%   s.coeff: scaling coefficient such that the background variance
-%     diag(inv(iB)) is one far away from the boundary.
-
-function s = divand_background(mask,pmn,Labs,alpha,moddim,mapindex)
-
-sz = size(mask);
-% number of dimensions
-pmn = cat_cell_array(pmn);
-n = size(pmn,ndims(pmn));
-
-if isempty(moddim)
-  moddim = zeros(1,n);
-end
-
-if nargin == 5 
-    mapindex = [];
-end
-
-iscyclic = moddim > 0;
-
-
-if isscalar(Labs)
-    Labs = Labs * ones([size(mask) n]);
-elseif iscell(Labs)
-    for i=1:n
-        if isscalar(Labs{i})
-            Labs{i} = Labs{i} * ones(size(mask));
-        else
-            if ~isequal(size(mask),size(Labs{i}))
-                error('mask (%s) and correlation length (%s) have incompatible size',...
-                      formatsize(size(mask)),formatsize(size(Labs{i})));
-            end            
-        end
-    end
-    
-    Labs = cat_cell_array(Labs);
-else
-    if ~isequal(size(mask),size(Labs))
-        error('mask (%s) and correlation length (%s) have incompatible size',...
-              formatsize(size(mask)),formatsize(size(Labs)));
-    end
-    
-    Labs = repmat(Labs,[ones(1,ndims(mask)) n]);
-end
-
-if isempty(alpha)
-  % kernel should has be continuous derivative
-  
-  % highest derivative in cost function
-  m = ceil(1+n/2);
-  
-  % alpha is the (m+1)th row of the Pascal triangle:
-  % m=0         1
-  % m=1       1   1
-  % m=1     1   2   1
-  % m=2   1   3   3   1
-  % ...
-  
-  alpha = arrayfun (@(k) nchoosek(m,k), 0:m);
-end
-
-%if ~isequal([size(mask) n],size(pmn))
-%  error('mask (%s) and metric (%s) have incompatible size',formatsize(size(mask)),formatsize(size(pmn)));
-%end
-
-s = divand_operators(mask,pmn,Labs.^2,iscyclic,mapindex);
-D = s.D; % laplacian (a dimensional, since nu = Labs.^2)
-sv = s.sv;
-n = s.n;
-
-
-% Labsp: 1st index represents the dimensions
-Labsp = permute(Labs,[n+1 1:n]);
-pmnp = permute(pmn,[n+1 1:n]);
-
-% must handle the case when Labs is zero in some dimension
-% thus reducing the effective dimension
-
-% mean correlation length in every dimension
-Ld = mean(reshape(Labsp,[n sv.numels_all]),2);
-neff = sum(Ld > 0);
-%L = mean(Ld(Ld > 0));
-
-% geometric mean
-L = geomean(Ld(Ld > 0));
-
-
-% norm taking only dimension into account with non-zero correlation
-% WE: units length^(neff/2)
-d = prod(pmnp(find(Ld > 0),:),1)';
-WE = sparse_diag(statevector_pack(sv,1./sqrt(d)));
-
-
-% scale iB such that the diagonal of inv(iB) is 1 far from
-% the boundary
-% we use the effective dimension neff to take into account that the
-% correlation length-scale might be zero in some directions
-
-Ln = prod(Ld(Ld > 0));
-
-if any(Ld <= 0)   
-   pmnd = mean(reshape(pmnp,[n sv.numels_all]),2);
-   %Ln = Ln * prod(pmnd(Ld <= 0)); 
-end
-
-
-coeff = divand_kernel(neff,alpha) * Ln; % units length^n
-
-pmnv = reshape(pmn,numel(mask),n);
-pmnv(:,Ld == 0) = 1;
-
-% staggered version of norm
-
-for i=1:n
-  S = sparse_stagger(sz,i,iscyclic(i));
-  ma = (S * mask(:) == 1);
-  d = sparse_pack(ma) * prod(S * pmnv,2);
-  d = 1./d;
-  s.WEs{i} = sparse_diag(sqrt(d));
-end
-
-% staggered version of norm scaled by length-scale
-
-s.WEss = cell(n,1);
-s.Dxs = cell(n,1);
-for i=1:n
-  Li2 = Labsp(i,:).^2;
-  
-  S = sparse_stagger(sz,i,iscyclic(i));
-  % mask for staggered variable
-  m = (S * mask(:) == 1);
-  
-  tmp = sparse_pack(m) * sqrt(S*Li2(:));
-  
-  s.WEss{i} = sparse_diag(tmp) * s.WEs{i};
-%  s.Dxs{i} = sparse_diag(sqrt(tmp)) * s.Dx{i};
-end
-
-% adjust weight of halo points
-if ~isempty(mapindex)
-  % ignore halo points at the center of the cell
-
-  WE = sparse_diag(s.isinterior) * WE;
-
-  % divide weight be two at the edged of halo-interior cell
-    % weight of the grid points between halo and interior points 
-    % are 1/2 (as there are two) and interior points are 1    
-  for i=1:n
-    s.WEs{i} = sparse_diag(sqrt(s.isinterior_stag{i})) * s.WEs{i};  
-  end
-end
-
-s.WE = WE;
-s.coeff = coeff;
-% number of dimensions
-s.n = n;
-
-[iB_,iB] = divand_background_components(s,alpha);
-
-if 0
-iB_ = cell(length(alpha),1);
-
-% constrain of total norm
-
-iB_{1} =  (1/coeff) * (WE'*WE);
-
-
-% loop over all derivatives
-
-for j=2:length(alpha)    
-    % exponent of laplacian
-    k = floor((j-2)/2);
-        
-    if mod(j,2) == 0
-        % constrain of derivative with uneven order (j-1)
-        % (gradient, gradient*laplacian,...)
-        % normalized by surface
-                
-        iB_{j} = sparse(size(D,1),size(D,1));
-        
-        for i=1:n
-            Dx = s.WEss{i} * s.Dx{i} * D^k;
-            
-            iB_{j} = iB_{j} + Dx'*Dx;
-        end                
-    else
-        % constrain of derivative with even order (j-1)
-        % (laplacian, biharmonic,...)
-        
-        % normalize by surface of each cell
-        % such that inner produces (i.e. WE'*WE)
-        % become integrals
-        % WD: units length^(n/2)
-        
-        WD = WE * D^(k+1);
-        iB_{j} = WD'*WD;
-    end
-    
-    iB_{j} = iB_{j}/coeff;
-end
-
-
-% sum all terms of iB
-% iB is adimentional
-iB = alpha(1) * iB_{1};
-for j=2:length(alpha)    
-  iB = iB + alpha(j) * iB_{j};
-end
-
-s.coeff = coeff;
-end
-
-% inverse of background covariance matrix
-s.iB = iB;
-
-% individual terms of background covariance matrix (corresponding alpha)
-s.iB_ = iB_;
-
-s.L = L;
-s.Ln = Ln;
-
-s.moddim = moddim;
-s.iscyclic = iscyclic;
-
-s.alpha = alpha;
-s.neff = neff;
-s.WE = WE; % units length^(n/2)
-% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program is distributed in the hope that it will be useful, but WITHOUT
-% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-% details.
-%
-% You should have received a copy of the GNU General Public License along with
-% this program; if not, see <http://www.gnu.org/licenses/>.