* Bug 15873 fixed: sprand() yielded wrong density
[scilab.git] / scilab / modules / sparse / tests / unit_tests / sprand.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
4 // Copyright (C) 2018 - Samuel GOUGEON
5 //
6 //  This file is distributed under the same license as the Scilab package.
7 // =============================================================================
8 // <-- CLI SHELL MODE -->
9 // <-- NO CHECK REF -->
10
11 // Unitary tests of sprand()
12 // =========================
13
14 // Check default "typ" = "uniform"
15 grand("setsd",0);
16 nrows = 1000;
17 ncols = 2000;
18 density = 1/100;
19 s=sprand(nrows,ncols,density);
20 assert_checkequal ( size(s) , [nrows,ncols] );
21 nnzs=nnz(s);
22 [ij,values]=spget(s);
23 assert_checkequal ( min(values) >= 0 , %t );
24 assert_checkequal ( max(values) <= 1 , %t );
25 assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
26
27 // Get empty matrix
28 s=sprand(0,0,0.01);
29 assert_checkequal ( s , sparse([],[],[0,0]) );
30
31 // Test the scientific part
32 // In the following script, we check that the entries of the matrix have the expected distribution.
33 // We use the spget function in order to get the nonzero entries.
34 // Then we compute the min, mean and max of the entries and compare them with the limit values.
35
36 grand("setsd",0);
37 typ = "normal";
38 nrows = 1000;
39 ncols = 2000;
40 density = 1/100;
41 s=sprand(nrows,ncols,density,typ);
42 assert_checkequal ( size(s) , [nrows,ncols] );
43 nnzs=nnz(s);
44 [ij,values]=spget(s);
45 assert_checkequal ( nnzs > 10000 , %t );
46 assert_checkequal ( min(values) < -3 , %t );
47 assert_checkalmostequal ( mean(values) , 0 , [] , 1.e-2 );
48 assert_checkequal ( max(values) > 3 , %t );
49 assert_checkalmostequal ( variance(values) , 1 , 2.2e-2 );
50
51 grand("setsd",0);
52 typ = "uniform";
53 nrows = 1000;
54 ncols = 2000;
55 density = 1/100;
56 s=sprand(nrows,ncols,density,typ);
57 assert_checkequal ( size(s) , [nrows,ncols] );
58 nnzs=nnz(s);
59 [ij,values]=spget(s);
60 assert_checkequal ( nnzs > 10000 , %t );
61 assert_checkalmostequal ( min(values) , 0 , [] , 1.e-2 );
62 assert_checkalmostequal ( mean(values) , 0.5 , 1.e-2 );
63 assert_checkalmostequal ( max(values) , 1 , 1.e-2 );
64 assert_checkalmostequal ( variance(values) , 1/12 , 1.e-2 );
65
66 // In the following script, we check that the entry indices,
67 // which are also chosen at random, have the correct distribution.
68 // We generate kmax sparse random matrices with uniform distribution.
69 // For each matrix, we consider the indices of the nonzero entries
70 // which were generated, i.e. we see if the event Aij = {the entry (i,j) is nonzero}
71 // occurred for each i and j, for i=1,2,...,nrows and j=1,2,...,ncols.
72 // The matrix C(i,j) stores the number of times that the event Aij occurred.
73 // The matrix R(k) stores the actual density of the try number k, where k=1,2,...,kmax.
74
75 grand("setsd",0);
76 kmax = 1000;
77 ncols=100;
78 nrows=100;
79 density=0.01;
80 typ="uniform";
81 C=zeros(nrows,ncols);
82 R=[];
83 for k=1:kmax
84   M=sprand(nrows,ncols,density,typ);
85   NZ=find(M<>0);
86   NZratio = size(NZ,"*")/(nrows*ncols);
87   R=[R NZratio];
88   C(NZ)=C(NZ)+1;
89 end
90 // Now that this algorithm has been performed (which may require some time),
91 // we can compute elementary statistics to check that the algorithm performed well.
92
93 // The average number should be close to the expectation.
94 assert_checkalmostequal ( density*kmax , mean(C) , 1.e-2 );
95 // The density should be close to expected density
96 assert_checkalmostequal ( density , mean(R) , 1.e-2 );
97
98 // More deeper tests should involve the particular distribution of
99 // C, which follows a binomial law.
100 // May be a chi-square test should be used for this.
101
102 // CHECK THE ACTUAL DENSITY
103 // ------------------------
104 clc
105 grand("setsd",0);
106 for n = [100 1000 1e4 1e5 1e6]
107     for d = [1e-5 3e-5 1e-4 3e-4 1e-3 0.003 0.01 0.03 0.1 0.3 1]
108         clear m
109         for i = 1:min(5e6/n, 50)
110             s = sprand(1,n,d);
111             m(i) = nnz(s);
112         end
113         assert_checkalmostequal(mean(m), d*n, 1e-3, 2);
114     end
115 end