* Bug #12130 fixed - Elementary_functions: block flipping for flipdim. 92/13692/4
Samuel Gougeon [Mon, 10 Feb 2014 07:59:49 +0000 (08:59 +0100)]
flipdim can now flip the input by blocks. jmat is then obsolete.

Change-Id: Id0ebb3adf968b501a291930d4bd5d2fca3d3283d

13 files changed:
SEP/INDEX
SEP/SEP_127_jmat_flipdim.odt [new file with mode: 0644]
scilab/CHANGES_5.5.X
scilab/modules/elementary_functions/help/en_US/matrixmanipulation/flipdim.xml
scilab/modules/elementary_functions/help/fr_FR/matrixmanipulation/flipdim.xml
scilab/modules/elementary_functions/macros/flipdim.sci
scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.dia.ref [new file with mode: 0644]
scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.tst [new file with mode: 0644]
scilab/modules/signal_processing/help/en_US/miscellaneous/jmat.xml
scilab/modules/signal_processing/macros/jmat.sci
scilab/modules/signal_processing/macros/lattn.sci
scilab/modules/signal_processing/macros/lattp.sci
scilab/modules/signal_processing/macros/levin.sci

index 6ea00fa..21b6ed5 100644 (file)
--- a/SEP/INDEX
+++ b/SEP/INDEX
@@ -122,3 +122,4 @@ SEP #123: Grid improvements
 SEP #124: New input argument for csvRead(), to ignore header comments.
 SEP #125: New input argument to matfile_open to open Matlab 7.3 files.
 SEP #126: New function unwrap
+SEP #127: jmat tagged as obsolete. New input argument to flipdim to replace jmat.
diff --git a/SEP/SEP_127_jmat_flipdim.odt b/SEP/SEP_127_jmat_flipdim.odt
new file mode 100644 (file)
index 0000000..a449c19
Binary files /dev/null and b/SEP/SEP_127_jmat_flipdim.odt differ
index 42e2a8f..46a7c3a 100644 (file)
@@ -52,6 +52,9 @@ Obsolete & Removed Functions
 * %asn tagged as obsolete. Will be removed in Scilab 5.5.1.
   Please use delip instead.
 
+* jmat tagged as obsolete. Will be removed in Scilab 5.5.1.
+  Please use flipdim instead
+
 
 Compilation
 ============
@@ -273,6 +276,8 @@ Scilab Bug Fixes
 
 * Bug #12121 fixed - inv function did not work for complex sparse matrices.
 
+* Bug #12130 fixed - flipdim can now flip blocks, thus making jmat obsolete.
+
 * Bug #12145 fixed - Internal function demo_mdialog removed.
 
 * Bug #12156 fixed - Closing a Scilab session in Javasci could lead to a HDF5 error message.
index 4af9927..013ebc9 100644 (file)
@@ -2,7 +2,7 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- * Copyright (C) 2013 - Samuel GOUGEON : restriction to decimal numbers removed. Examples added for booleans, integer-encoded numbers, text, polynomials, rationals
+ * Copyright (C) 2013 - Samuel GOUGEON: restriction to decimal numbers removed. Examples added for booleans, integer-encoded numbers, text, polynomials, rationals
  *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="flipdim" xml:lang="en">
     <refnamediv>
         <refname>flipdim</refname>
-        <refpurpose>flip x components along a given dimension</refpurpose>
+        <refpurpose>flip x block components along a given dimension</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>y = flipdim(x, dim)</synopsis>
+        <synopsis>y = flipdim(x, dim [,sb])</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
         <variablelist>
             <varlistentry>
                 <term>x, y</term>
-                
                 <listitem>
                     <para>scalars, vectors, matrices, or hypermatrices of any type, of same sizes</para>
-                    
                 </listitem>
             </varlistentry>
             <varlistentry>
                     <para>a positive integer</para>
                 </listitem>
             </varlistentry>
+            <varlistentry>
+                <term>sb</term>
+                <listitem>
+                    <para>a positive integer, the size of the blocks to permute</para>
+                </listitem>
+            </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
         <title>Description</title>
         <para>
             Given <literal>x</literal>, a scalar/vector/matrix/hypermatrix of any type and
-            <literal>dim</literal> a positive integer, this function flips the x
-            components along the dimension number <literal>dim</literal> of
-            <literal>x</literal> (<literal>x</literal> and <literal>y</literal> have
-            the same size)
+            two positive integers <literal>dim</literal> and <literal>sb</literal>,
+            this function flips the x components by blocks of size <literal>sb</literal>
+            along the dimension number <literal>dim</literal> of <literal>x</literal>
+            (<literal>x</literal> and <literal>y</literal> have the same size).
+        </para>
+        <para>
+            The optional parameter <literal>sb</literal> (for Size Block) allows flipping
+            <literal>x</literal> by blocks of size <literal>sb*size(x,2)</literal>
+            (<literal>dim=1</literal>) or <literal>size(x,1)*sb</literal> (<literal>dim=2</literal>).
         </para>
     </refsection>
     <refsection>
@@ -93,7 +102,19 @@ n = inv_coeff(grand(3, 9, "uin", 0, 3), 2);
 d = inv_coeff(grand(3, 9, "uin", 0, 3), 2);
 r = n./d
 flipdim(r, 2)
+ ]]></programlisting>
+        <para>
+            Examples using <literal>sb</literal>:
+        </para>
+        <programlisting role="example"><![CDATA[
+X = [0 1 2 3 4 5 6 7 8 9 10 11];
+flipdim(X, 2, 2) // => [10 11   8 9   6 7   4 5   2 3   0 1] // Block size = 2.
+flipdim(X, 2, 3) // => [9 10 11   6 7 8   3 4 5   0 1 2]
+flipdim(X, 2, 4) // => [8 9 10 11   4 5 6 7   0 1 2 3]
+flipdim(X, 2, 6) // => [6 7 8 9 10 11   0 1 2 3 4 5]
 
+// Error if sb does not divide the targeted dimension of x.
+y = flipdim(x, 2, 5)", refMsg); // size(X) = [1 12] and sb=5 does not divide 12.
  ]]></programlisting>
     </refsection>
     <refsection>
@@ -101,7 +122,10 @@ flipdim(r, 2)
         <revhistory>
             <revision>
                 <revnumber>5.5.0</revnumber>
-                <revremark>Extension from decimals to any type: booleans, integer-encoded numbers, texts, polynomials and rationals.</revremark>
+                <revremark>
+                    Extension from decimals to any type: booleans, integers, strings, polynomials and rationals.
+                    New input argument <literal>sb</literal> to flip <literal>x</literal> blockwise.
+                </revremark>
             </revision>
         </revhistory>
     </refsection>
index a59ef30..b9bfd5b 100644 (file)
@@ -2,7 +2,7 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- * Copyright (C) 2013 - Samuel GOUGEON : restriction to decimal numbers removed. Examples added for booleans, integer-encoded numbers, text, polynomials, rationals
+ * Copyright (C) 2013 - Samuel GOUGEON: restriction to decimal numbers removed. Examples added for booleans, integer-encoded numbers, text, polynomials, rationals
  *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
     <refnamediv>
         <refname>flipdim</refname>
         <refpurpose>
-            retourne les éléments de <literal>x</literal> selon une dimension
+            retourne par blocs les éléments de <literal>x</literal> selon une dimension
         </refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Séquence d'appel</title>
-        <synopsis>y = flipdim(x, dim)</synopsis>
+        <synopsis>y = flipdim(x, dim [,sb])</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Paramètres</title>
                     <para>entier positif</para>
                 </listitem>
             </varlistentry>
+            <varlistentry>
+                <term>sb</term>
+                <listitem>
+                    <para>entier positif, la taille des blocs à permuter</para>
+                </listitem>
+            </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
         <title>Description</title>
         <para>
             A partir de <literal>x</literal>, un scalaire/vecteur/matrice/hypermatrice de n'importe quel type et
-            <literal>dim</literal> un entier positif, cette fonction retourne les éléments de
-            <literal>x</literal> selon le nombre dimension <literal>dim</literal> de
+            deux entiers positifs <literal>dim</literal> et <literal>sb</literal>, cette fonction retourne les éléments de
+            <literal>x</literal> par blocs de taille <literal>sb</literal> selon le nombre dimension <literal>dim</literal> de
             <literal>x</literal> (<literal>x</literal> et <literal>y</literal> ont la même taille)
         </para>
+        <para>
+            Le paramètre optionnel <literal>sb</literal> (pour Size Block) permet de retounerles éléments de <literal>x</literal>
+            par blocs de taille <literal>sb*size(x,2)</literal> (<literal>dim=1</literal>)
+            ou <literal>size(x,1)*sb</literal> (<literal>dim=2</literal>).
+        </para>
     </refsection>
     <refsection>
         <title>Exemples</title>
@@ -93,13 +104,29 @@ d = inv_coeff(grand(3, 9, "uin", 0, 3), 2);
 r = n./d
 flipdim(r, 2)
  ]]></programlisting>
+        <para>
+            Exemples utilisant <literal>sb</literal> :
+        </para>
+        <programlisting role="example"><![CDATA[
+X = [0 1 2 3 4 5 6 7 8 9 10 11];
+flipdim(X, 2, 2) // => [10 11   8 9   6 7   4 5   2 3   0 1] // Taille du bloc = 2.
+flipdim(X, 2, 3) // => [9 10 11   6 7 8   3 4 5   0 1 2]
+flipdim(X, 2, 4) // => [8 9 10 11   4 5 6 7   0 1 2 3]
+flipdim(X, 2, 6) // => [6 7 8 9 10 11   0 1 2 3 4 5]
+
+// Erreur si sb ne divise pas la dimension ciblée de x.
+y = flipdim(x, 2, 5)", refMsg); // size(X) = [1 12] et sb=5 ne divise pas 12.
+ ]]></programlisting>
     </refsection>
     <refsection>
         <title>History</title>
         <revhistory>
             <revision>
                 <revnumber>5.5.0</revnumber>
-                <revremark>Extension de décimaux à tout type: booléens, entiers, textes, polynômes et rationals.</revremark>
+                <revremark>
+                    Extension de décimaux à tout type : booléens, entiers, chaines de caractères, polynômes et fractions rationnelles.
+                    Nouveau paramètre optionnel <literal>sb</literal> pour retourner <literal>x</literal> par blocs.
+                </revremark>
             </revision>
         </revhistory>
     </refsection>
index c07e0ea..2c46e6d 100644 (file)
@@ -9,22 +9,23 @@
 // are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 
-function y = flipdim(x, dim)
+function y = flipdim(x, dim, sb)
 
     // FLIPDIM function
-    // Given x, a scalar/vector/matrix of any type and an integer dim, this function flips the x components  along the dimension number dim (x and y have the same size).
+    // Given x, a scalar/vector/matrix of any type, two integers dim and sb, this function flips the x components by blocks along the dimension number dim (x and y have the same size).
     // -Inputs :
     //  x : a scalar/vector/array
     //  dim : a positive integer
+    //  sb : size of the block to permute
     // -Output :
     //  y : a scalar/vector/array
     //
     // F.Belahcene
 
     rhs = argn(2);
-    if rhs <> 2 then
-        msg = _("%s: Wrong number of input argument(s): %d expected.\n");
-        error(msprintf(msg, "flipdim", 2));
+    if rhs < 2 then
+        msg = _("%s: Wrong number of input argument(s): %d to %d expected.\n");
+        error(msprintf(msg, "flipdim", 2, 3));
     end
 
     if size(dim, "*") <> 1 then
@@ -34,23 +35,58 @@ function y = flipdim(x, dim)
         msg = _("%s: Wrong type for input argument #%d: A positive integer expected.\n");
         error(msprintf(msg, "flipdim", 2));
     end
+    if rhs >= 3 then
+        if size(sb, "*") <> 1 then
+            msg = _("%s: Wrong size for input argument #%d: A positive integer expected.\n")
+            error(msprintf(msg, "flipdim", 3));
+        elseif type(sb) <> 8 & (type(sb) <> 1 | sb < 1 ) then
+            msg = _("%s: Wrong type for input argument #%d: A positive integer expected.\n");
+            error(msprintf(msg, "flipdim", 3));
+        elseif dim > 2 then
+            msg = _("%s: Cannot flip hypermatrix blockwise. %d input arguments expected.\n");
+            error(msprintf(msg, "flipdim", 2));
+        end
+    else
+        sb = 1;
+    end
 
     dim = floor(dim);
+    sb  = floor(sb);
 
     if dim > ndims(x)
         y = x;
         return
-    end
+    elseif or(dim == [1 2]) & ndims(x) < 3 then  // flipdim(x, dim, sb) is not suited for hypermatrices.
 
-    l = list();
-    for k = 1:dim - 1
-        l(k) = eye();
-    end
-    l(dim) = $:-1:1;
-    for k = dim + 1:ndims(x)
-        l(k) = eye();
-    end
+        if dim == 1 then
+            x = x.';
+        end
 
-    y = x(l(:));
+        nC = size(x, 2);
+        if modulo(nC, sb) ~= 0 then
+            msg = _("%s: Wrong value for input argument #%d: A divisor of the selected dimension size expected.\n");
+            error(msprintf(msg, "flipdim", 3));
+        end
+        nb = nC/sb; // Number of blocks.
+        c2 = ((nb:-1:1).*.ones(1,sb))*sb + ones(1, nb).*.(1-sb:0);
+        y = x(:, c2);
+        if dim == 1 then
+            y = y.';
+        end
+
+    else
+
+        l = list();
+        for k = 1:dim - 1
+            l(k) = eye();
+        end
+        l(dim) = $:-1:1;
+        for k = dim + 1:ndims(x)
+            l(k) = eye();
+        end
+
+        y = x(l(:));
+
+    end
 
 endfunction
diff --git a/scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.dia.ref b/scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.dia.ref
new file mode 100644 (file)
index 0000000..65863ab
--- /dev/null
@@ -0,0 +1,61 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2014 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 12130 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12130
+//
+// <-- Short Description -->
+// flipdim can now flip an input by blocks.
+//---------------------------------------------------
+// Double type
+x = [0 1 2 3 4 5 6 7 8 9 10 11];
+x = [x ; x];
+ref1 = [11 10 9 8 7 6 5 4 3 2 1 0];
+ref1 = [ref1 ; ref1];
+y = flipdim(x, 2); // Present action.
+assert_checkequal(y, ref1);
+ref2 = [10 11   8 9   6 7   4 5   2 3   0 1];
+ref2 = [ref2 ; ref2];
+y = flipdim(x, 2, 2); // Block size = 2.
+assert_checkequal(y, ref2);
+ref3 = [9 10 11   6 7 8   3 4 5   0 1 2];
+ref3 = [ref3 ; ref3];
+y = flipdim(x, 2, 3);
+assert_checkequal(y, ref3);
+ref4 = [8 9 10 11   4 5 6 7   0 1 2 3];
+ref4 = [ref4 ; ref4];
+y = flipdim(x, 2, 4);
+assert_checkequal(y, ref4);
+ref5 = [6 7 8 9 10 11   0 1 2 3 4 5];
+ref5 = [ref5 ; ref5];
+y = flipdim(x, 2, 6);
+assert_checkequal(y, ref5);
+//---------------------------------------------------
+// String type
+x = string(x);
+ref1 = string(ref1);
+y = flipdim(x, 2); // Present action.
+assert_checkequal(y, ref1);
+ref2 = string(ref2);
+y = flipdim(x, 2, 2); // Block size = 2.
+assert_checkequal(y, ref2);
+ref3 = string(ref3);
+y = flipdim(x, 2, 3);
+assert_checkequal(y, ref3);
+ref4 = string(ref4);
+y = flipdim(x, 2, 4);
+assert_checkequal(y, ref4);
+ref5 = string(ref5);
+y = flipdim(x, 2, 6);
+assert_checkequal(y, ref5);
+// Error checks
+refMsg = msprintf(_("%s: Wrong value for input argument #%d: A divisor of the selected dimension size expected.\n"), "flipdim", 3);
+assert_checkerror("y = flipdim(x, 2, 5)", refMsg); // size(x) = [2 12] and 5 does not divide 12.
diff --git a/scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.tst b/scilab/modules/elementary_functions/tests/nonreg_tests/bug_12130.tst
new file mode 100644 (file)
index 0000000..e59a612
--- /dev/null
@@ -0,0 +1,77 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2014 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 12130 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12130
+//
+// <-- Short Description -->
+// flipdim can now flip an input by blocks.
+
+
+//---------------------------------------------------
+// Double type
+x = [0 1 2 3 4 5 6 7 8 9 10 11];
+x = [x ; x];
+
+ref1 = [11 10 9 8 7 6 5 4 3 2 1 0];
+ref1 = [ref1 ; ref1];
+y = flipdim(x, 2); // Present action.
+assert_checkequal(y, ref1);
+
+ref2 = [10 11   8 9   6 7   4 5   2 3   0 1];
+ref2 = [ref2 ; ref2];
+y = flipdim(x, 2, 2); // Block size = 2.
+assert_checkequal(y, ref2);
+
+ref3 = [9 10 11   6 7 8   3 4 5   0 1 2];
+ref3 = [ref3 ; ref3];
+y = flipdim(x, 2, 3);
+assert_checkequal(y, ref3);
+
+ref4 = [8 9 10 11   4 5 6 7   0 1 2 3];
+ref4 = [ref4 ; ref4];
+y = flipdim(x, 2, 4);
+assert_checkequal(y, ref4);
+
+ref5 = [6 7 8 9 10 11   0 1 2 3 4 5];
+ref5 = [ref5 ; ref5];
+y = flipdim(x, 2, 6);
+assert_checkequal(y, ref5);
+
+
+//---------------------------------------------------
+// String type
+x = string(x);
+ref1 = string(ref1);
+y = flipdim(x, 2); // Present action.
+assert_checkequal(y, ref1);
+
+ref2 = string(ref2);
+y = flipdim(x, 2, 2); // Block size = 2.
+assert_checkequal(y, ref2);
+
+ref3 = string(ref3);
+y = flipdim(x, 2, 3);
+assert_checkequal(y, ref3);
+
+ref4 = string(ref4);
+y = flipdim(x, 2, 4);
+assert_checkequal(y, ref4);
+
+ref5 = string(ref5);
+y = flipdim(x, 2, 6);
+assert_checkequal(y, ref5);
+
+
+// Error checks
+refMsg = msprintf(_("%s: Wrong value for input argument #%d: A divisor of the selected dimension size expected.\n"), "flipdim", 3);
+assert_checkerror("y = flipdim(x, 2, 5)", refMsg); // size(x) = [2 12] and 5 does not divide 12.
+
index 1d4209a..050c05b 100644 (file)
@@ -2,7 +2,9 @@
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="jmat">
     <refnamediv>
         <refname>jmat</refname>
-        <refpurpose>row or column block permutation</refpurpose>
+        <refpurpose>row or column block permutation.
+            <emphasis role="bold">This function is obsolete.</emphasis>
+        </refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
@@ -46,5 +48,21 @@ m*j   // permutes blocks of columns
 j'*m  // permutes blocks of rows
  ]]></programlisting>
     </refsection>
-    
+    <refsection role="see also">
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="flipdim">flipdim</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revremark>Tagged as obsolete, please use flipdim instead. Will be removed in Scilab 5.5.1.</revremark>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index f29cb56..b509917 100644 (file)
@@ -14,6 +14,7 @@ function [j]=jmat(n,m)
     //   n   : number of block rows or block columns of the matrix
     //   m   : size of the (square) blocks
     //!
+    warnobsolete("flipdim", "5.5.1");
     j=[];
     for k=1:n,j(k,n-k+1)=1;end;
     j=j.*.eye(m,m);
index 82d1563..fb262e4 100644 (file)
@@ -77,9 +77,13 @@ function [la,lb]=lattn(n,p,mat_cov)
     end
 
 
-    a=eye(d);b=eye(d);
-    z=poly(0,"z");la=list();lb=list();
-    no=p-n-1;cv=mat_cov;
+    a=eye(d);
+    b=eye(d);
+    z=poly(0,"z");
+    la=list();
+    lb=list();
+    no=p-n-1;
+    cv=mat_cov;
 
     if -no >= floor(l/d) then
         error(msprintf(gettext("%s: Wrong values for input arguments #%d and #%d.\n"),"lattn",1, 2))
@@ -95,16 +99,18 @@ function [la,lb]=lattn(n,p,mat_cov)
     end
 
     for j=0:n-1,
-        jd=jmat(j+1,d);
-        r1=jd*cv((p+1)*d+1:(p+2+j)*d,:);
-        r2=jd*cv(p*d+1:(p+1+j)*d,:);
-        r3=jd*cv((p-1-j)*d+1:p*d,:);
-        r4=jd*cv((p-j)*d+1:(p+1)*d,:);
-        c1=coeff(a);c2=coeff(b);
+        r1=flipdim(cv((p+1)*d+1:(p+2+j)*d,:), 1, d);
+        r2=flipdim(cv(p*d+1:(p+1+j)*d,:), 1, d);
+        r3=flipdim(cv((p-1-j)*d+1:p*d,:), 1, d);
+        r4=flipdim(cv((p-j)*d+1:(p+1)*d,:), 1, d);
+        c1=coeff(a);
+        c2=coeff(b);
         k1=(c1*r1)*inv(c2*r2);
         k2=(c2*r3)*inv(c1*r4);
         a1=a-k1*z*b;
-        b=-k2*a+z*b;a=a1;
-        la(j+1)=a;lb(j+1)=b;
+        b=-k2*a+z*b;
+        a=a1;
+        la(j+1)=a;
+        lb(j+1)=b;
     end;
 endfunction
index d62b83f..790ba38 100644 (file)
@@ -12,13 +12,17 @@ function [la,lb]=lattp(n,p,mat_cov)
     // See lattn for more information
     [l,d]=size(mat_cov);
     id=eye(d);
-    [a,b]=lattn(n,0,mat_cov);a=a(n);b=b(n);
-    z=poly(0,"z");la=list();lb=list();
-    jd=jmat(n+1,d);
+    [a,b]=lattn(n,0,mat_cov);
+    a=a(n);
+    b=b(n);
+    z=poly(0,"z");
+    la=list();
+    lb=list();
     for j=0:p-1,
-        r1=jd*mat_cov((j+1)*d+1:(j+n+2)*d,:);
-        r2=jd*mat_cov(j*d+1:(j+n+1)*d,:);
-        c1=coeff(a);c2=coeff(b);
+        r1=flipdim(mat_cov((j+1)*d+1:(j+n+2)*d,:), 1, d);
+        r2=flipdim(mat_cov(j*d+1:(j+n+1)*d,:), 1, d);
+        c1=coeff(a);
+        c2=coeff(b);
         k=(c1*r1)*inv(c2*r2);
         hst=-inv(c1(:,n*d+1:(n+1)*d));
         r=k*hst;
index 754c382..f2f2eab 100644 (file)
@@ -67,23 +67,23 @@ function [la, sig, lb] = levin(n, Cov)
     end
     for j=0:n-1
         //
-        //   Bloc permutation matrix
-        //
-        jd = jmat(j+1, d);
-        //
         //   Levinson algorithm
         //
-        r1 = jd*cv((p+1)*d+1:(p+2+j)*d, :);
-        r2 = jd*cv(p*d+1:(p+1+j)*d, :);
-        r3 = jd*cv((p-1-j)*d+1:p*d, :);
-        r4 = jd*cv((p-j)*d+1:(p+1)*d, :);
-        c1 = coeff(a); c2 = coeff(b);
-        sig1 = c1*r4; gam1 = c2*r2;
+        r1 = flipdim(cv((p+1)*d+1:(p+2+j)*d, :), 1, d);
+        r2 = flipdim(cv(p*d+1:(p+1+j)*d, :),     1, d);
+        r3 = flipdim(cv((p-1-j)*d+1:p*d, :),     1, d);
+        r4 = flipdim(cv((p-j)*d+1:(p+1)*d, :),   1, d);
+        c1 = coeff(a);
+        c2 = coeff(b);
+        sig1 = c1*r4;
+        gam1 = c2*r2;
         k1 = (c1*r1)*inv(gam1);
         k2 = (c2*r3)*inv(sig1);
         a1 = a-k1*z*b;
         b = -k2*a+z*b;
         a = a1;
-        la(j+1) = a; lb(j+1) = b; sig(j+1) = sig1;
+        la(j+1) = a;
+        lb(j+1) = b;
+        sig(j+1) = sig1;
     end
 endfunction