Mercurial > octave-nkf
annotate scripts/sparse/private/__sprand_impl__.m @ 20266:83792dd9bcc1
Use in-place operators in m-files where possible.
* scripts/audio/@audioplayer/set.m, scripts/audio/@audiorecorder/set.m,
scripts/audio/mu2lin.m, scripts/elfun/cosd.m, scripts/general/del2.m,
scripts/general/profexplore.m, scripts/general/quadl.m, scripts/general/rat.m,
scripts/general/rotdim.m, scripts/help/get_first_help_sentence.m,
scripts/help/private/__strip_html_tags__.m, scripts/image/cubehelix.m,
scripts/io/textread.m, scripts/linear-algebra/duplication_matrix.m,
scripts/linear-algebra/housh.m, scripts/linear-algebra/krylov.m,
scripts/linear-algebra/logm.m, scripts/linear-algebra/normest.m,
scripts/linear-algebra/onenormest.m, scripts/optimization/fminsearch.m,
scripts/optimization/lsqnonneg.m, scripts/optimization/qp.m,
scripts/plot/appearance/annotation.m, scripts/plot/appearance/axis.m,
scripts/plot/appearance/legend.m, scripts/plot/appearance/specular.m,
scripts/plot/draw/colorbar.m, scripts/plot/draw/hist.m,
scripts/plot/draw/plotmatrix.m, scripts/plot/draw/private/__stem__.m,
scripts/plot/util/__actual_axis_position__.m,
scripts/plot/util/__gnuplot_drawnow__.m, scripts/plot/util/findobj.m,
scripts/plot/util/print.m, scripts/plot/util/private/__go_draw_axes__.m,
scripts/plot/util/private/__print_parse_opts__.m, scripts/plot/util/rotate.m,
scripts/polynomial/pchip.m, scripts/polynomial/polyaffine.m,
scripts/polynomial/polyder.m, scripts/polynomial/private/__splinefit__.m,
scripts/polynomial/residue.m, scripts/signal/arch_fit.m,
scripts/signal/arch_rnd.m, scripts/signal/bartlett.m,
scripts/signal/blackman.m, scripts/signal/freqz.m, scripts/signal/hamming.m,
scripts/signal/hanning.m, scripts/signal/spectral_adf.m,
scripts/signal/spectral_xdf.m, scripts/signal/stft.m,
scripts/sparse/bicgstab.m, scripts/sparse/cgs.m,
scripts/sparse/private/__sprand_impl__.m, scripts/sparse/qmr.m,
scripts/sparse/sprandsym.m, scripts/sparse/svds.m, scripts/specfun/legendre.m,
scripts/special-matrix/gallery.m, scripts/statistics/base/gls.m,
scripts/statistics/models/logistic_regression.m,
scripts/statistics/tests/kruskal_wallis_test.m,
scripts/statistics/tests/manova.m, scripts/statistics/tests/wilcoxon_test.m,
scripts/time/datevec.m:
Use in-place operators in m-files where possible.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 26 May 2015 21:07:42 -0700 |
parents | 4197fc428c7d |
children |
rev | line source |
---|---|
19731
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19630
diff
changeset
|
1 ## Copyright (C) 2004-2015 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; | |
20266
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
138 S *= V; |
18624 | 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 |