Mercurial > octave-nkf
annotate scripts/sparse/private/__sprand_impl__.m @ 19630:0e1f5a750d00
maint: Periodic merge of gui-release to default.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 20 Jan 2015 10:24:46 -0500 |
| parents | 54a1e95365e1 446c46af4b42 |
| children | 4197fc428c7d |
| rev | line source |
|---|---|
|
17744
d63878346099
maint: Update copyright notices for release.
John W. Eaton <jwe@octave.org>
parents:
17712
diff
changeset
|
1 ## Copyright (C) 2004-2013 Paul Kienzle |
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13929
diff
changeset
|
2 ## Copyright (C) 2012 Jordi Gutiérrez Hermoso |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
3 ## |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
4 ## This file is part of Octave. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
5 ## |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
7 ## under the terms of the GNU General Public License as published by |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
9 ## your option) any later version. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
10 ## |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
14 ## General Public License for more details. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
15 ## |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
16 ## You should have received a copy of the GNU General Public License |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
17 ## along with Octave; see the file COPYING. If not, see |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
18 ## <http://www.gnu.org/licenses/>. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
19 ## |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
20 ## Original version by Paul Kienzle distributed as free software in the |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
21 ## public domain. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
22 |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
23 ## -*- texinfo -*- |
|
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13202
diff
changeset
|
24 ## @deftypefn {Function File} {} __sprand_impl__ (@var{s}, @var{randfun}) |
| 18624 | 25 ## @deftypefnx {Function File} {} __sprand_impl__ (@var{m}, @var{n}, @var{d}, @var{fcnname}, @var{randfun}) |
| 26 ## @deftypefnx {Function File} {} __sprand_impl__ (@var{m}, @var{n}, @var{d}, @var{rc}, @var{fcnname}, @var{randfun}) | |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
27 ## Undocumented internal function. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
28 ## @end deftypefn |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
29 |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
30 ## Actual implementation of sprand and sprandn happens here. |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
31 |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
32 function S = __sprand_impl__ (varargin) |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
33 |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
34 if (nargin == 2) |
| 18624 | 35 [m, randfun] = deal (varargin{1:2}); |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
36 [i, j] = find (m); |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
37 [nr, nc] = size (m); |
|
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
38 S = sparse (i, j, randfun (size (i)), nr, nc); |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
39 else |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
40 if (nargin == 5) |
| 18624 | 41 [m, n, d, fcnname, randfun] = deal (varargin{:}); |
| 42 else | |
| 43 [m, n, d, rc, fcnname, randfun] = deal (varargin{:}); | |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
44 endif |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
45 |
| 18624 | 46 if (! (isscalar (m) && m == fix (m) && m > 0)) |
| 47 error ("%s: M must be an integer greater than 0", fcnname); | |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
48 endif |
| 18624 | 49 if (! (isscalar (n) && n == fix (n) && n > 0)) |
| 50 error ("%s: N must be an integer greater than 0", fcnname); | |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
51 endif |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
52 if (d < 0 || d > 1) |
| 18624 | 53 error ("%s: density D must be between 0 and 1", fcnname); |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
54 endif |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
55 |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
56 if (nargin == 5) |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
57 mn = m*n; |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
58 k = round (d*mn); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
59 if (mn > sizemax ()) |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
60 ## randperm will overflow, so use alternative methods |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
61 |
| 18624 | 62 idx = unique (fix (rand (1.01*k, 1) * mn)) + 1; |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
63 |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
64 ## idx contains random numbers in [1,mn] |
| 18624 | 65 ## Generate 1% more random values than necessary in order to reduce the |
| 66 ## probability that there are less than k distinct values; maybe a | |
| 67 ## better strategy could be used but I don't think it's worth the price. | |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
68 |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
69 ## actual number of entries in S |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
70 k = min (length (idx), k); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
71 j = floor ((idx(1:k) - 1) / m); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
72 i = idx(1:k) - j * m; |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
73 j++; |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
74 else |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
75 idx = randperm (mn, k); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
76 [i, j] = ind2sub ([m, n], idx); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
77 endif |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
78 |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
79 S = sparse (i, j, randfun (k, 1), m, n); |
| 18624 | 80 |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
81 elseif (nargin == 6) |
| 18624 | 82 ## Create a matrix with specified reciprocal condition number. |
| 83 | |
| 84 if (! isscalar (rc) && ! isvector (rc)) | |
| 85 error ("%s: RC must be a scalar or vector", fcnname); | |
| 86 endif | |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
87 |
| 18624 | 88 ## We want to reverse singular valued decomposition A=U*S*V'. |
| 89 ## First, first S is constructed and then U = U1*U2*..Un and | |
| 90 ## V' = V1*V2*..Vn are seen as Jacobi rotation matrices with angles and | |
| 91 ## planes of rotation randomized. Repeatedly apply rotations until the | |
| 92 ## required density for A is achieved. | |
|
16779
8fce0ed4894a
Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents:
14868
diff
changeset
|
93 |
| 18624 | 94 if (isscalar (rc)) |
| 95 if (rc < 0 || rc > 1) | |
| 96 error ("%s: reciprocal condition number RC must be between 0 and 1", fcnname); | |
| 97 endif | |
| 98 ## Reciprocal condition number is ratio of smallest SV to largest SV | |
| 99 ## Generate singular values randomly and sort them to build S | |
| 100 ## Random singular values in range [rc, 1]. | |
| 101 v = rand (1, min (m,n)) * (1 - rc) + rc; | |
| 102 v(1) = 1; | |
| 103 v(end) = rc; | |
| 104 v = sort (v, "descend"); | |
| 105 S = sparse (diag (v, m, n)); | |
| 106 else | |
| 107 ## Only the min (m, n) greater singular values from rc vector are used. | |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
108 if (length (rc) > min (m,n)) |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
109 rc = rc(1:min(m, n)); |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
110 endif |
| 18624 | 111 S = sparse (diag (sort (rc, "descend"), m, n)); |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
112 endif |
|
13197
6db186dfdeaa
Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff
changeset
|
113 |
| 18624 | 114 Uinit = speye (m); |
| 115 Vinit = speye (n); | |
| 116 k = round (d*m*n); | |
| 117 while (nnz (S) < k) | |
| 118 if (m > 1) | |
| 119 ## Construct U randomized rotation matrix | |
| 120 rot_angleu = 2 * pi * rand (); | |
| 121 cu = cos (rot_angleu); su = sin (rot_angleu); | |
| 122 rndtmp = randperm (m, 2); | |
| 123 i = rndtmp(1); j = rndtmp(2); | |
| 124 U = Uinit; | |
| 125 U(i, i) = cu; U(i, j) = -su; | |
| 126 U(j, i) = su; U(j, j) = cu; | |
| 127 S = U * S; | |
| 128 endif | |
| 129 if (n > 1) | |
| 130 ## Construct V' randomized rotation matrix | |
| 131 rot_anglev = 2 * pi * rand (); | |
| 132 cv = cos (rot_anglev); sv = sin (rot_anglev); | |
| 133 rndtmp = randperm (n, 2); | |
| 134 i = rndtmp(1); j = rndtmp(2); | |
| 135 V = Vinit; | |
| 136 V(i, i) = cv; V(i, j) = sv; | |
| 137 V(j, i) = -sv; V(j, j) = cv; | |
| 138 S = S * V; | |
| 139 endif | |
|
18623
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
140 endwhile |
|
35a5e7740a6d
Added implementation for 4th argument of sprand/sprandn (bug #41839).
Eduardo Ramos (edu159) <eduradical951@gmail.com>
parents:
17744
diff
changeset
|
141 endif |
|
16779
8fce0ed4894a
Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents:
14868
diff
changeset
|
142 endif |
|
8fce0ed4894a
Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents:
14868
diff
changeset
|
143 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
144 endfunction |
|
17338
1c89599167a6
maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents:
17176
diff
changeset
|
145 |