* Bug 15873 fixed: sprand() yielded wrong density 12/20612/5
Samuel GOUGEON [Sat, 10 Nov 2018 20:32:29 +0000 (21:32 +0100)]
   http://bugzilla.scilab.org/15873

   Actual density relative accuracy more than 100 times better.
   Processing time slightly smaller (~ -20%)

   test_run sparse // passes, except spisp (as without this patch)

Change-Id: I1c303f272bf3e247660da2cf4f401860cbcd49bb

scilab/CHANGES.md
scilab/modules/sparse/macros/sprand.sci
scilab/modules/sparse/tests/unit_tests/sprand.dia.ref [deleted file]
scilab/modules/sparse/tests/unit_tests/sprand.tst

index 0f5dd32..92989b5 100644 (file)
@@ -541,7 +541,7 @@ Known issues
 * [#15284](http://bugzilla.scilab.org/show_bug.cgi?id=15284): Port names are not set to the corresponding I/O block labels.
 * [#15310](http://bugzilla.scilab.org/show_bug.cgi?id=15310): `isdef` considered void arguments as defined (regression)
 * [#15346](http://bugzilla.scilab.org/show_bug.cgi?id=15346): In an array of rationals, there was no way to address components with their linearized indices.
-* [#15403](http://bugzilla.scilab.org/show_bug.cgi?id=15403): `bar(..,"stacked")` could start from y<>0. Input arguments and possible syntaxes were poorly tested. 
+* [#15403](http://bugzilla.scilab.org/show_bug.cgi?id=15403): `bar(..,"stacked")` could start from y<>0. Input arguments and possible syntaxes were poorly tested.
 * [#15404](http://bugzilla.scilab.org/show_bug.cgi?id=15404): `surf()` and `mesh()` did not allow to specify `foreground`, `facecolor`, `markforeground` and `markbackground` global properties colors as a predefined named color out of a list of the 9 main color names. Colors specifications as "#RRGGBB" hexa code or Colors indices in the color map were nor allowed.
 * [#15422](http://bugzilla.scilab.org/show_bug.cgi?id=15422): `strsubst("ab", "", "cd")` crashed Scilab.
 * [#15423](http://bugzilla.scilab.org/show_bug.cgi?id=15423): `tbx_make(myModule,sections)` executed the existing builder (if any), instead of targeting only selected module sections. Otherwise, tbx_make(myModule, "help"|"macros") yielded an error, and tbx_make(myModule,"localization") never built it.
@@ -686,6 +686,7 @@ Known issues
 * [#15860](http://bugzilla.scilab.org/show_bug.cgi?id=15860): Horizontal concatenations with cblock tables had troubles.
 * [#15866](http://bugzilla.scilab.org/show_bug.cgi?id=15866): Unlike all other set functions, `setdiff` did not accept any "r" or "c" option.
 * [#15867](http://bugzilla.scilab.org/show_bug.cgi?id=15867): For input encoded integers, `setdiff` returned selected elements in decreasing order.
+* [#15873](http://bugzilla.scilab.org/show_bug.cgi?id=15873): In average, `sprand(100, 100, 0.8)` yielded ~8800 non-zero values instead of ~8000.
 * [#15878](http://bugzilla.scilab.org/show_bug.cgi?id=15878): `sgrid` and `evans` were broken.
 * [#15880](http://bugzilla.scilab.org/show_bug.cgi?id=15880): `sgrid` needed some improvements: Labeling was sometimes ambiguous ; large circles were not labeled ; data_bounds did not always take the input wn into account ; named and #RRGGBB colors specifications could not be used. `evans` needed some improvements: the block of legends hid data ; asymptotes were too visible.
 * [#15881](http://bugzilla.scilab.org/show_bug.cgi?id=15881): Demos about stems, bars and histograms were spread in several subsections.
index 2ad65c6..227ab0a 100644 (file)
@@ -1,8 +1,8 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA
 // Copyright (C) 2010 - DIGITEO - Michael Baudin
-//
 // 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.
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function a=sprand(m,n,density,typ)
+function a = sprand(m, n, density, law)
 
-    //---- Check arguments----------------------------------------------
+    // CHECK ARGUMENTS
+    // ---------------
     rhs = argn(2)
-
-    if ( rhs < 3 | rhs > 4 ) then
+    if rhs < 3 | rhs > 4 then
         error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"), "sprand" , 3 , 4 ));
     end
-
-    if ( rhs < 4 ) then
-        typ="def";
+    if rhs < 4 then
+        law = "unf";
     end
-
-    if type(typ)<>10 then
+    if m==0 | n==0 then
+        a = sparse([], [], [0 0])
+        return
+    end
+    // law
+    if type(law)<>10 then
         error(msprintf(gettext("%s: Wrong type for input argument #%d: string expected.\n"),"sprand",4));
     end
-
-    if and(typ<>["u";"n";"uniform";"normal";"def";"nor"]) then
+    if and(law<>["u" "unf" "uniform" "def" "n" "nor" "normal"]) then
         error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or ''%s'' expected.\n"),"sprand",4,"uniform","normal"));
     end
-    if typ == "u" | typ == "uniform" then //"uniform" is the syntax for uniform distribution with rand, the equivalent with grand is "def"
-        typ = "def";
-    elseif typ == "n" | typ == "normal" then //"normal" is the syntax for normal distribution with rand, the equivalent with grand is "nor"
-        typ = "nor";
+    if or(law == ["u" "uniform" "def"]) then
+        law = "unf";
+    elseif or(law == ["n" "normal"]) then
+        law = "nor";
     end
 
-    density=max(min(density,1),0);
-
-    nel=m*n*density; //the objective number of non zero elements
+    // INITIALIZATION
+    // --------------
+    N = m*n
+    density = max(min(density,1),0);
+    nel = round(N*density); // the objective number of non zero elements
     if nel==0 then
-        a=sparse([],[],[m,n])
-        return
-    end
-
-    //---- generate a sequence of increments----------------------------
-    mdist = 1/density //the mean distance between to consecutive index
-    nel1  = (2.2-density)*nel; //generate more increments than requested nnz elements
-    if nel1<1
-        ij = 1;
-    else
-        ij    = round(1+grand(nel1,1,"exp",(mdist-1)))
-    end
-
-    //---- sum the increments to get the index--------------------------
-    ij=cumsum(ij);
-
-    //---- eliminate the index with exceed the maximum matrix index
-    ij(find(ij>m*n))=[];
-    nel1=size(ij,"*");
-    if nel1==0 then
-        a=sparse([],[],[m,n]);
+        a = sparse([],[], [m,n])
         return
     end
-    //---- generate the row and column indices from the 1D index--------
-
-    ij=ind2sub([m,n],ij);
-
-    //----  generates the random non zeros elements --------------------
-    //according to the requested law and create the sparse matrix
-    if typ == "nor" then // Because of the syntax of grand, we have two cases, one with "def" and one with "nor"
-        a=sparse(ij,grand(nel1,1,typ,0,1),[m,n]);
-    elseif typ == "def" then
-        a=sparse(ij,grand(nel1,1,typ),[m,n]);
-    end
-
 
+    // GENERATION OF RANDOM UNIQUE INDICES OF NON-ZERO VALUES
+    // ------------------------------------------------------
+    // generate a sequence of increments greater than 1 ------------------------
+    di = N / nel
+    mdist = 1/density //the mean distance between two consecutive non-zero values
+    ij    = grand(nel*1.01, 1, "exp", (mdist-1))
+
+    s = (nel + sum(ij(1:nel)))/(N-di/2)  // rescaling factor
+    ij = 1 + ij/s              // rescales intervals & guaranties that indices are unique
+    ij = floor(cumsum(ij));    // Intervals between indices => indices
+    ij(find(ij > N)) = [];     // Garanties that linear indices are <= m*n
+    nel = length(ij)
+    ij = ind2sub([m,n], ij);   // => [I J] conversion
+
+    // RESULT
+    // ------
+    a = sparse(ij, grand(nel, 1, law, 0, 1),[m,n]);
 endfunction
 
diff --git a/scilab/modules/sparse/tests/unit_tests/sprand.dia.ref b/scilab/modules/sparse/tests/unit_tests/sprand.dia.ref
deleted file mode 100644 (file)
index 97da5c6..0000000
+++ /dev/null
@@ -1,90 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// Check default "typ" = "uniform"
-rand("seed",0);
-grand("setsd",0);
-nrows = 1000;
-ncols = 2000;
-density = 1/100;
-s=sprand(nrows,ncols,density);
-assert_checkequal ( size(s) , [nrows,ncols] );
-nnzs=nnz(s);
-[ij,values]=spget(s);
-assert_checkequal ( min(values) >= 0 , %t );
-assert_checkequal ( max(values) <= 1 , %t );
-assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
-// Get empty matrix
-s=sprand(0,0,0.01);
-assert_checkequal ( s , [] );
-// Test the scientific part
-// In the following script, we check that the entries of the matrix have the expected distribution. 
-// We use the spget function in order to get the nonzero entries. 
-// Then we compute the min, mean and max of the entries and compare them with the limit values.
-rand("seed",0);
-grand("setsd",0);
-typ = "normal";
-nrows = 1000;
-ncols = 2000;
-density = 1/100;
-s=sprand(nrows,ncols,density,typ);
-assert_checkequal ( size(s) , [nrows,ncols] );
-nnzs=nnz(s);
-[ij,values]=spget(s);
-assert_checkequal ( nnzs > 10000 , %t );
-assert_checkequal ( min(values) < -3 , %t );
-assert_checkalmostequal ( mean(values) , 0 , [] , 1.e-2 );
-assert_checkequal ( max(values) > 3 , %t );
-assert_checkalmostequal ( variance(values) , 1 , 1.e-2 );
-rand("seed",0);
-grand("setsd",0);
-typ = "uniform";
-nrows = 1000;
-ncols = 2000;
-density = 1/100;
-s=sprand(nrows,ncols,density,typ);
-assert_checkequal ( size(s) , [nrows,ncols] );
-nnzs=nnz(s);
-[ij,values]=spget(s);
-assert_checkequal ( nnzs > 10000 , %t );
-assert_checkalmostequal ( min(values) , 0 , [] , 1.e-2 );
-assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
-assert_checkalmostequal ( max(values) , 1 , 1.e-2 );
-assert_checkalmostequal ( variance(values) , 1/12 , 1.e-2 );
-// In the following script, we check that the entry indices, 
-// which are also chosen at random, have the correct distribution. 
-// We generate kmax sparse random matrices with uniform distribution. 
-// For each matrix, we consider the indices of the nonzero entries 
-// which were generated, i.e. we see if the event Aij = {the entry (i,j) is nonzero}
-// occurred for each i and j, for i=1,2,...,nrows and j=1,2,...,ncols. 
-// The matrix C(i,j) stores the number of times that the event Aij occurred. 
-// The matrix R(k) stores the actual density of the try number k, where k=1,2,...,kmax.
-rand("seed",0);
-grand("setsd",0);
-kmax = 1000;
-ncols=100;
-nrows=100;
-density=0.01;
-typ="uniform";
-C=zeros(nrows,ncols);
-R=[];
-for k=1:kmax
-  M=sprand(nrows,ncols,density,typ);
-  NZ=find(M<>0);
-  NZratio = size(NZ,"*")/(nrows*ncols);
-  R=[R NZratio];
-  C(NZ)=C(NZ)+1;
-end
-// Now that this algorithm has been performed (which may require some time), 
-// we can compute elementary statistics to check that the algorithm performed well.
-// The average number should be close to the expectation.
-assert_checkalmostequal ( density*kmax , mean(C) , 1.e-2 );
-// The density should be close to expected density
-assert_checkalmostequal ( density , mean(R) , 1.e-2 );
-// More deeper tests should involve the particular distribution of 
-// C, which follows a binomial law.
-// May be a chi-square test should be used for this.
index 2ded8e7..dd66075 100644 (file)
@@ -1,15 +1,17 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+// Copyright (C) 2018 - Samuel GOUGEON
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 // <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
 
-
+// Unitary tests of sprand()
+// =========================
 
 // Check default "typ" = "uniform"
-rand("seed",0);
 grand("setsd",0);
 nrows = 1000;
 ncols = 2000;
@@ -24,14 +26,13 @@ assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
 
 // Get empty matrix
 s=sprand(0,0,0.01);
-assert_checkequal ( s , [] );
+assert_checkequal ( s , sparse([],[],[0,0]) );
 
 // Test the scientific part
-// In the following script, we check that the entries of the matrix have the expected distribution. 
-// We use the spget function in order to get the nonzero entries. 
+// In the following script, we check that the entries of the matrix have the expected distribution.
+// We use the spget function in order to get the nonzero entries.
 // Then we compute the min, mean and max of the entries and compare them with the limit values.
 
-rand("seed",0);
 grand("setsd",0);
 typ = "normal";
 nrows = 1000;
@@ -45,9 +46,8 @@ assert_checkequal ( nnzs > 10000 , %t );
 assert_checkequal ( min(values) < -3 , %t );
 assert_checkalmostequal ( mean(values) , 0 , [] , 1.e-2 );
 assert_checkequal ( max(values) > 3 , %t );
-assert_checkalmostequal ( variance(values) , 1 , 1.e-2 );
+assert_checkalmostequal ( variance(values) , 1 , 2.2e-2 );
 
-rand("seed",0);
 grand("setsd",0);
 typ = "uniform";
 nrows = 1000;
@@ -63,16 +63,15 @@ assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
 assert_checkalmostequal ( max(values) , 1 , 1.e-2 );
 assert_checkalmostequal ( variance(values) , 1/12 , 1.e-2 );
 
-// In the following script, we check that the entry indices, 
-// which are also chosen at random, have the correct distribution. 
-// We generate kmax sparse random matrices with uniform distribution. 
-// For each matrix, we consider the indices of the nonzero entries 
+// In the following script, we check that the entry indices,
+// which are also chosen at random, have the correct distribution.
+// We generate kmax sparse random matrices with uniform distribution.
+// For each matrix, we consider the indices of the nonzero entries
 // which were generated, i.e. we see if the event Aij = {the entry (i,j) is nonzero}
-// occurred for each i and j, for i=1,2,...,nrows and j=1,2,...,ncols. 
-// The matrix C(i,j) stores the number of times that the event Aij occurred. 
+// occurred for each i and j, for i=1,2,...,nrows and j=1,2,...,ncols.
+// The matrix C(i,j) stores the number of times that the event Aij occurred.
 // The matrix R(k) stores the actual density of the try number k, where k=1,2,...,kmax.
 
-rand("seed",0);
 grand("setsd",0);
 kmax = 1000;
 ncols=100;
@@ -88,7 +87,7 @@ for k=1:kmax
   R=[R NZratio];
   C(NZ)=C(NZ)+1;
 end
-// Now that this algorithm has been performed (which may require some time), 
+// Now that this algorithm has been performed (which may require some time),
 // we can compute elementary statistics to check that the algorithm performed well.
 
 // The average number should be close to the expectation.
@@ -96,8 +95,21 @@ assert_checkalmostequal ( density*kmax , mean(C) , 1.e-2 );
 // The density should be close to expected density
 assert_checkalmostequal ( density , mean(R) , 1.e-2 );
 
-// More deeper tests should involve the particular distribution of 
+// More deeper tests should involve the particular distribution of
 // C, which follows a binomial law.
 // May be a chi-square test should be used for this.
 
-
+// CHECK THE ACTUAL DENSITY
+// ------------------------
+clc
+grand("setsd",0);
+for n = [100 1000 1e4 1e5 1e6]
+    for d = [1e-5 3e-5 1e-4 3e-4 1e-3 0.003 0.01 0.03 0.1 0.3 1]
+        clear m
+        for i = 1:min(5e6/n, 50)
+            s = sprand(1,n,d);
+            m(i) = nnz(s);
+        end
+        assert_checkalmostequal(mean(m), d*n, 1e-3, 2);
+    end
+end