* Bug 15428 fixed: repmat() was slow: rewritten => > 7x faster 82/19782/16
Samuel GOUGEON [Fri, 16 Feb 2018 06:35:47 +0000 (07:35 +0100)]
  http://bugzilla.scilab.org/15428

  http://mailinglists.scilab.org/Scilab-users-More-rapid-calculation-tp4037596p4037608.html

Change-Id: I599bbfa3988dbf865f2ff32eb360eabddd53227b

scilab/CHANGES.md
scilab/modules/elementary_functions/macros/repmat.sci
scilab/modules/elementary_functions/tests/benchmarks/bench_repmat.tst [new file with mode: 0644]
scilab/modules/elementary_functions/tests/unit_tests/repmat.tst

index 20e0c96..b922a1d 100644 (file)
@@ -145,6 +145,7 @@ Feature changes and additions
 * `intersect()` now supports complex numbers.
 * `setdiff()` now supports complex numbers.
 * `twinkle` can now blink together several hierarchically independent objects, like a curve and its labels, etc.
+* `repmat()` has been rewritten. It is 7 times faster now.
 
 
 Help pages:
@@ -272,6 +273,7 @@ Bug Fixes
 * [#15392](http://bugzilla.scilab.org/show_bug.cgi?id=15392): `comet` and `comet3d` did not allow specifying colors with colors names.
 * [#15393](http://bugzilla.scilab.org/show_bug.cgi?id=15393): In a new figure, `nicholschart` plotted nothing. The default frame color was a flashy cyan. The position of gain labels could be puzzling. It was not possible to specify colors by their names. Postprocessing the frames and the set of labels was not easy.
 * [#15425](http://bugzilla.scilab.org/show_bug.cgi?id=15425): The Kronecker product `a.*.b` failed when `a` or `b` or both are hypermatrices, with one or both being polynomials or rationals.
+* [#15428](http://bugzilla.scilab.org/show_bug.cgi?id=15428): `repmat()` was slow. Its code did not use `kron()` and was complex.
 * [#15431](http://bugzilla.scilab.org/show_bug.cgi?id=15431): The empty matrix `[]` and its non-convertibility were poorly documented.
 * [#15451](http://bugzilla.scilab.org/show_bug.cgi?id=15451): The code was not adapted to use `lucene 4.10` in Debian.
 * [#15514](http://bugzilla.scilab.org/show_bug.cgi?id=15514): The `set()` documentation page needed to be overhauled.
index a6d0bda..fb3fe1c 100644 (file)
@@ -1,6 +1,7 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) - 2011 -  INRIA, Serge Steer
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) - 2018 - Samuel GOUGEON
 //
 // This file is hereby licensed under the terms of the GNU GPL v2.0,
 // pursuant to article 5.3.4 of the CeCILL v.2.1.
 // and continues to be available under such terms.
 // For more information, see the COPYING file which you should have received
 // along with this program.
+
 function B = repmat(A,varargin)
-    //ajouter test sur les type des elements de varargin
+
+    // CHECKING INPUT ARGUMENTS
+    // ------------------------
     rhs = argn(2);
     if rhs < 2 then
-        error(msprintf(_("%s: Wrong number of input arguments: at least %d expected.\n"), "repmat", 2))
+        msg = _("%s: Wrong number of input arguments: at least %d expected.\n")
+        error(msprintf(msg, "repmat", 2))
+    end
+
+    if A == [] then
+        B = A
+        return
     end
 
-    narg=size(varargin);
+    narg = size(varargin);
 
-    // Test varargin
-    if narg==1 then
+    // Reading & setting replication numbers
+    if narg == 1 then
         // Scalar of vector needed
-        if typeof(varargin(1)) <> "constant" then
-            error(msprintf(_("%s: Wrong type for input argument #%d: A real scalar or vector expected.\n"), "repmat", 2))
+        v = varargin(1);
+        if type(v) <> 1 | ~isreal(v, 0) | ~and(v == int(v)) | v == [] | v <0 then
+            msg = _("%s: Argument #%d: Non-negative integers expected.\n")
+            error(msprintf(msg, "repmat", 2))
+        else
+            v = real(v);    // in case of complex encoding with null imag
+        end
+        if ~isscalar(v) & ~isvector(v) then
+            msg = _("%s: Wrong size for input argument #%d: A scalar or vector expected.\n")
+            error(msprintf(msg, "repmat", 2))
         end
-        if size(varargin(1),"*")<>1 & isempty(find(size(varargin(1))==1)) then
-            error(msprintf(_("%s: Wrong size for input argument #%d: A real scalar or vector expected.\n"), "repmat", 2))
+        sizes = v;
+        if length(v) == 1 then
+            sizes = [v,v];
         end
     else
-        for i=1:narg
-            if typeof(varargin(i)) <> "constant" then
-                error(msprintf(_("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "repmat", i+1))
+        sizes = [];
+        for i = 1:narg
+            v = varargin(i);
+            if type(v) <> 1 | ~isreal(v, 0) | ~and(v == int(v)) | v == [] | v < 0 then
+                msg = _("%s: Argument #%d: Non-negative integers expected.\n")
+                error(msprintf(msg, "repmat", i+1))
             end
-            if size(varargin(i),"*")<>1 then
-                error(msprintf(_("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "repmat", i+1))
+            if length(v)<>1 then
+                msg = _("%s: Argument #%d: Scalar (1 element) expected.\n")
+                error(msprintf(msg, "repmat", i+1))
             end
+            sizes = [sizes v];
         end
     end
 
-    if type(A)>10 then
-        if typeof(A)=="rational" then
-            B=rlist(repmat(A.num,varargin(:)),repmat(A.den,varargin(:)),A.dt)
-            return
+    // Preprocessing sizes = replication factors
+    // -----------------------------------------
+    // Trimming ending ones in sizes
+    for s = length(sizes):-1:3
+        if sizes(s) == 1 then
+            sizes(s) = []
         else
-            execstr("B=%"+typeof(A)+"_repmat(A,varargin(:))")
-            return
+            break
         end
     end
-
-    if narg==1 then
-        if size(varargin(1),"*")==1 then
-            siz=list(varargin(1),varargin(1))
-        else //convert array into list
-            tmp=varargin(1)
-            siz=list();
-            for i=1:size(tmp,"*"),siz(i)=tmp(i); end
-        end
+    // B expected sizes
+    sA = size(A)
+    ndimsA = length(sA)
+    ndimsR = length(sizes)
+    if ndimsA >= ndimsR then
+        Bsizes = size(A).*[sizes ones(1,(ndimsA-ndimsR))]
     else
-        siz=list();
-        for i=1:narg
-            siz(i)=varargin(i)
-        end
+        Bsizes = [sA ones(1,(ndimsR-ndimsA))] .* sizes
     end
-
-    nd=size(siz)
-    if or(type(A)==[5 6]) then //sparse matrices
-        if nd>2 then
-            error(msprintf(_("%s: Wrong number of output matrix dimensions required: %d expected for sparse matrices.\n"), "repmat", 2))
-        end
+    // Sparse matrices are only 2D (limited output dimensions):
+    if or(type(A) == [5 6]) & length(Bsize) > 2 then
+        msg = _("%s: Wrong dimensions for input arguments: The sparse result must be 2D, not %dD.\n");
+        error(msprintf(msg, "repmat", length(Bsizes)));
     end
-
-    for i=size(siz):-1:3
-        if siz(i)>1 then break,end
-        nd=nd-1
+    // Replication factors
+    Rfactors = list();
+    for s = sizes
+        Rfactors($+1) = s
     end
-    sizA=size(A)
-    nda=size(sizA,"*")
 
-    if and(sizA==1) then //scalar case
+    // PROCESSING
+    // ----------
+    if find(sizes == 0) then
+        B = []
+        return
+    end
 
-        //this case can also be handled by the general one but in a less
-        //efficient way
-        if nd<=2 then
-            B=A(ones(siz(1:nd)))
+    // Special case where the input A is scalar: special efficient processing
+    if size(A,"*") == 1 then
+        i = uint8(0);
+        i(Rfactors(:)) = 0;
+        i = i + 1;
+        if typeof(A) <> "rational" then
+            B = matrix(A(i), sizes);
         else
-            s=1;for k=1:nd;s=s*siz(k),end
-            B=matrix(A(ones(s,1)),siz(1:nd))
-        end
-    else //general case
-        if nda<nd then
-            sizA(nda+1:nd)=1;
-        elseif  nda>nd then
-            for k=nd+1:nda
-                siz(k)=1
-            end
-            nd=nda
-        end
-        I=list();
-        for i=1:nd
-            ind=matrix(1:sizA(i),-1,1);
-            ind=ind(:,ones(1,siz(i)));
-            I(i)=ind;
+            B = matrix(A(i,0), sizes);
         end
+        return
+    end
 
-        if nda > 2 | (size(varargin(1),"*") <> 1 & size(varargin(1)) <3) then // Works if A is hypermat but not for int8,int16 matrix
-            B=A(I(:));
-        else // Works for int8, int16... matrix but not for hypermat
-            if rhs ==2 then
-                if size(varargin(1),"*") <> 1 then // case repmat(A,size)
-                    varargin_temp=varargin;
-                    if size(varargin(1),1) <> 1 & size(varargin(1),2) == 1 then
-                        for i=1:size(varargin(1),1)
-                            varargin(i)=varargin_temp(1)(i);
-                        end
-                    elseif size(varargin(1),2) <> 1 & size(varargin(1),1) == 1 then
-                        for i=1:size(varargin(1),2)
-                            varargin(i)=varargin_temp(1)(i);
-                        end
-                    else
-                        error(msprintf(_("%s: Wrong size for input argument #%d: a vector expected.\n"),"repmat",2));
-                    end
-                end
-            end
-            res=A(I(1),I(2));
-            A_base=matrix(res,size(res,1)*size(res,2),1);
-            dims=[size(A,1)*varargin(1)];
-            for i=2:size(varargin) //compute dims
-                dims_1=[size(A,i)*varargin(i)];
-                dims=[dims, dims_1];
-            end
-            nb=1;
-            flag=0;
-            for i=3:size(varargin)
-                nb=nb*varargin(i);
-                if varargin(i)~=1 then // dims[2 2 1 1]=>dims[2 2]
-                    flag=1;
-                end
-            end
-            if size(varargin) ==1 then // dims[2] => dims[2 2]
-                varargin(2) = varargin(1);
-            end
-            if flag == 0 then
-                B=matrix(A_base,varargin(1)*size(A,1), varargin(2)*size(A,2));
-            else
-                J1=[1:size(A_base,1)];
-                J=J1;
-                for k=1:nb-1
-                    J=[J,J1];
-                end
-                J=J';
-                A_final=A_base(J);
-                B=matrix(A_final, dims);
-            end
-        end
+    // Numerical data: double, integers, polynomials, sparse : kron() is used
+    if or(type(A) == [1 2 5 7 8]) then
+        B = ones(Rfactors(:)) .*. A;
+        return
+    end
+
+    // Other regular types : booleans, textes..:
+    // => we replicate and use the matrix of linearized indices of A components
+    if or(type(A) == [4 6 9 10]) then
+        // Processing (without intermediate variable, to save memory)
+        B = matrix(A(ones(Rfactors(:)) .*. matrix(1:size(A,"*"), size(A))), Bsizes);
+        return
     end
-endfunction
 
+    // Rationals, other mlists, overloads
+    if typeof(A) == "rational" then
+        B = rlist(repmat(A.num,sizes), repmat(A.den,sizes), A.dt)
+    else
+        execstr("B=%" + typeof(A) + "_repmat(A, sizes)")
+    end
+endfunction
diff --git a/scilab/modules/elementary_functions/tests/benchmarks/bench_repmat.tst b/scilab/modules/elementary_functions/tests/benchmarks/bench_repmat.tst
new file mode 100644 (file)
index 0000000..c00f70a
--- /dev/null
@@ -0,0 +1,22 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2018 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+//==============================================================================
+// Benchmark for repmat() function
+//==============================================================================
+
+// <-- BENCH NB RUN : 10 -->
+
+r = rand(1,10000);
+
+// <-- BENCH START -->
+
+repmat(r, 100,60);
+
+// <-- BENCH END -->
+
+clear ans
index b4c05a4..c9d1d17 100644 (file)
@@ -1,6 +1,7 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2011 - INRIA - Serge Steer
+// Copyright (C) 2018 - Samuel GOUGEON
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
@@ -31,3 +32,36 @@ assert_checkequal(repmat(matrix(1:2, [1 2]),2),[1 2 1 2; 1 2 1 2]);
 assert_checkequal(repmat([1,2;3,4],[2,3]),[1,2,1,2,1,2;3,4,3,4,3,4;1,2,1,2,1,2;3,4,3,4,3,4]);
 assert_checkequal(repmat(int8([1,2;3,4]),[2,3]),int8([1,2,1,2,1,2;3,4,3,4,3,4;1,2,1,2,1,2;3,4,3,4,3,4]));
 assert_checkequal(repmat(matrix(1:2, 1, 2),[2,3]),[1,2,1,2,1,2;1,2,1,2,1,2]);
+
+// With other hypermatrices
+    // Polynomials
+H = repmat([%z %z^2], [1 2 2]);
+res = cat(3, [%z %z^2 %z %z^2], [%z %z^2 %z %z^2]);
+assert_checkequal(H,res);
+H = repmat(cat(3,[%z %z^2],[1 %z]), 1, 2);
+res = cat(3, [%z %z^2 %z %z^2], [1 %z 1 %z]);
+assert_checkequal(H,res);
+    // Rationals
+H = repmat(1./[%z %z^2], [1 2 2]);
+res = cat(3, 1./[%z %z^2 %z %z^2], 1./[%z %z^2 %z %z^2]);
+assert_checkequal(H,res);
+H = repmat(1./cat(3,[%z %z^2],[1 %z]), 1, 2);
+res = cat(3, 1./[%z %z^2 %z %z^2], 1./[1 %z 1 %z]);
+assert_checkequal(H,res);
+    // Booleans
+H = repmat([%t %f ; %f %t], [1 2 2]);
+res = cat(3, [%t %f %t %f; %f %t %f %t], [%t %f %t %f; %f %t %f %t]);
+assert_checkequal(H,res);
+H = repmat(cat(3,[%f %t %t],[%t %f %t]), 2, 1);
+res = cat(3, [%f %t %t ; %f %t %t], [%t %f %t; %t %f %t]);
+assert_checkequal(H,res);
+    // Texts
+H = repmat(["a" "b" ; "c" "d"], [2 1 2]);
+res = cat(3, ["a" "b" ; "c" "d" ; "a" "b" ; "c" "d"], ..
+             ["a" "b" ; "c" "d" ; "a" "b" ; "c" "d"]);
+assert_checkequal(H, res);
+H = repmat(cat(3,"a", "b", "c"), 2,3,2);
+res = cat(3,["a" "a" "a"; "a" "a" "a"],["b" "b" "b"; "b" "b" "b"],["c" "c" "c"; "c" "c" "c"]);
+res = cat(3, res, res);
+assert_checkequal(H, res);
+