227ab0af4027c31ee6e87b862333da09d2505f47
[scilab.git] / scilab / modules / sparse / macros / sprand.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 // Copyright (C) 2010 - DIGITEO - Michael Baudin
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 // Copyright (C) 2018 - Samuel GOUGEON
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 function a = sprand(m, n, density, law)
15
16     // CHECK ARGUMENTS
17     // ---------------
18     rhs = argn(2)
19     if rhs < 3 | rhs > 4 then
20         error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"), "sprand" , 3 , 4 ));
21     end
22     if rhs < 4 then
23         law = "unf";
24     end
25     if m==0 | n==0 then
26         a = sparse([], [], [0 0])
27         return
28     end
29     // law
30     if type(law)<>10 then
31         error(msprintf(gettext("%s: Wrong type for input argument #%d: string expected.\n"),"sprand",4));
32     end
33     if and(law<>["u" "unf" "uniform" "def" "n" "nor" "normal"]) then
34         error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or ''%s'' expected.\n"),"sprand",4,"uniform","normal"));
35     end
36     if or(law == ["u" "uniform" "def"]) then
37         law = "unf";
38     elseif or(law == ["n" "normal"]) then
39         law = "nor";
40     end
41
42     // INITIALIZATION
43     // --------------
44     N = m*n
45     density = max(min(density,1),0);
46     nel = round(N*density); // the objective number of non zero elements
47     if nel==0 then
48         a = sparse([],[], [m,n])
49         return
50     end
51
52     // GENERATION OF RANDOM UNIQUE INDICES OF NON-ZERO VALUES
53     // ------------------------------------------------------
54     // generate a sequence of increments greater than 1 ------------------------
55     di = N / nel
56     mdist = 1/density //the mean distance between two consecutive non-zero values
57     ij    = grand(nel*1.01, 1, "exp", (mdist-1))
58
59     s = (nel + sum(ij(1:nel)))/(N-di/2)  // rescaling factor
60     ij = 1 + ij/s              // rescales intervals & guaranties that indices are unique
61     ij = floor(cumsum(ij));    // Intervals between indices => indices
62     ij(find(ij > N)) = [];     // Garanties that linear indices are <= m*n
63     nel = length(ij)
64     ij = ind2sub([m,n], ij);   // => [I J] conversion
65
66     // RESULT
67     // ------
68     a = sparse(ij, grand(nel, 1, law, 0, 1),[m,n]);
69 endfunction
70