Re: [R] writing spdiags function for R
, (Adata(i), i=1,mn)
do 10 i=1,mn
Bdata(i) = Adata(i)
10 continue
kk1 = (k-1)*k
do 15 i=1,kk1
Adata(i)=0.0
15 continue
js = kk1+1
je = kk1+mn
do 20 i=js,je
Adata(i) = Bdata(i-js+1)
20 continue
js = je+1
je = je+kk1
do 25 i=js,je
Adata(i) = 0.0
25 continue
C print *, Augmented Adata with , je, elements
C print 1000,(Adata(i), i=1,je)
1000 FORMAT(1H ,25f4.0)
do 100 i=1,q
qx = k*(i-1)+1
C print *,Top element is no ,qx,=,Adata(qx)
tv(1) = Adata(qx)
do 30 j=1,(k-1)
tv(j+1) = 0.0
30 continue
qx = min(q-i, k-1)
C print *,qx=,qx
if (qx .gt. 0) then
C qx will be 0 when we are at last superdiagonal,i.e., i == q
do 35 j=1,qx
tv(j+1) = Adata(k*(i+j-1)+j+1)
35 continue
endif
not0 = .FALSE.
do 40 j=1,k
if (tv(j) .ne. 0.0) not0 = .TRUE.
40 continue
C if (.NOT. not0) print *, Zero column for i=,i
if (not0) then
C print *,Nonzero column for i=,i
jb = jb+1
d(jb) = (i-k)
C print *, New jb=,jb, d element =,d(jb)
C print 1000,(tv(j), j=1,k)
do 45 j=1,k
Bdata(kend+j) = tv(j)
45 continue
kend = kend + k
C save the diagonal as column of B in vector form
endif
100 continue
C print *,Bdata , jb, columns
C qx = k*jb
C print 1000,(Bdata(j), j=1,qx)
return
end
#
R Wrapper
#
spDIAGS - function(A) { # A is a matrix
m - dim(A)[[1]]
n - dim(A)[[2]]
k - min(m, n) # length of diagonals
if (m n) { # tall matrix
Adata - as.vector(t(A)) # convert to vector BY ROWS
} else { # fat or square matrix (m = n)
Adata - as.vector(A) # convert to vector BY COLUMNS
}
# print(Adata)
la-length(Adata)
na-m*n+2*k*(k-1)
Adata - c(Adata, rep(0,(na-la)))
nb-(m+n)*k
nd-m+n-1
jb-0
Bdata-rep(0,nb)
d-rep(0,nd)
tv-rep(0,k)
tres-.Fortran(jspd,m=as.integer(m), n=as.integer(n), k=as.integer(k),
Adata=as.double(Adata), jb=as.integer(jb),
Bdata=as.double(Bdata), d=as.integer(d), tv=as.double(tv),
na=as.integer(na), nb=as.integer(nb), nd=as.integer(nd) )
# print(str(tres))
jb-tres$jb
d-tres$d[1:jb]
Bdata-tres$Bdata[1:(jb*k)]
# print(Bdata)
if (m n) d - -d # reset index
B - matrix(Bdata, nrow=k, byrow=FALSE) # convert to matrix form
result-list(B=B, d=d)
}
#
--
View this message in context:
http://r.789695.n4.nabble.com/writing-spdiags-function-for-R-tp4552626p4675010.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
Dear R-list,
l have been working on a translation of a matlab library into R, it took me
a while but I am almost there to submit it to CRAN...however, for this
library to be computationally competitive I need to solve an issue with the
home-made R version of the spdiags.m function (an old issue coming back).
What I have now working is:
# A = rbind(c(1,3,1,0), c(1,1,0,2), c(2,1,0,0) )
# A = replicate(100, rnorm(50))
A = replicate(1000, rnorm(1000))
exdiag - function(mat, off) {mat[row(mat)+off == col(mat)]}
spdiags - function(A){
require(Matrix)
indx = which(A != 0, arr.ind = T)
i = indx[,1]; j = indx[,2]
d = sort(j-i);
d = d[which(diff(c(-Inf, d) ) != 0)]
m = nrow(A); n = ncol(A);
p = length(d)
empty = vector()
A = as.matrix(A) ## make A logical, otherwise exdiag throws horrible
warnings
## might have to condition the above command on the type of input matrix
B = Matrix(0, nrow = min(c(m,n)), ncol = p, sparse = TRUE);
system.time(
for (k in 1:p){
print(k)
if (m = n){
i = max(c(1, 1 + d[k])):min(c(n, m + d[k]) )
} else {
i = max(c(1, 1 - d[k])):min(c(m, n-d[k]) )
}
if (length(i) != 0){
B[i, k] = exdiag(A, d[k])
# print(B)
# print(i)
# print(k)
}
}
)
return (list( B = B, d = d) )
}
the advantages are: 1) that it works with matrices of any dimension
(previous versions worked only with squared matrices). 2) it does not need
other functions, beside exdiag to work.
However, I run into serious issues with computational efficiency. As I am
extracting non-zero diagonals and filling B in a for loop, the time it takes
to complete each iteration highly depends on the dimensions of the matrix.
Try the different starting A matrices to test this.
Unfortunately, I really need to optimize this aspect, as the people who I am
expecting to use this library would have really long time-series, i.e. it
would make .m win over .R
So, as I am wondering if there is a more efficient solution to do this
(perhaps using apply?) that I cannot currently see. Perhaps using sparse
matrices could help, but exdiag throws a warning: sparse[ logic ] :
.M.sub.i.logical() maybe inefficient (noted also in previous responses to
this post of mine).
Using the function ?band has the drawback of padding the matrix with 0
(except the diagonal considered, e.g., band(A,1,1) ). Often, however, you
want to keep the 0s along the diagonal of interest, distinguished from the
0s of the other diagonals. Thus, doing something like: which(band(A,1,1) !=
0, arr.ind = TRUE) to get the non-zero element would also discard the 0s
along the diagonal.
The nice getband,R posted by Ben above, makes a combined use of diag and
band, together with a system of indeces, and tries to work around the issue
of the zero padding. However, I did not manage to use/adapt this function
to work with matrices that are not squared. And, as I am not really 100% a
programmer, I did not find a solution to it...
Moreover, I believe the problem of computational efficiency would still be
there, even fixing the dimensionality issue.
I really hope somebody out there can help me to figure this out.
thanks a lot for your precious time and attention,
Moreno
--
View this message in context:
http://r.789695.n4.nabble.com/writing-spdiags-function-for-R-tp4552626p4674540.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
Hi Ben, A couple of quick comments: the semicolons at the end of each line are unnecessary (they're harmless, but considered bad style -- most common in code of C and MATLAB coders yes, I know, just forgot to clean it up after translating :) You don't need to source(getband.R) inside the body of spdiags, probably better to do it outside (so it doesn't get run every time the function is called) if (length(which(bd) == T) != 0) can be shortened to if (any(bd!=0)) I haven't looked too carefully but can't all the i/j/d/p stuff be replaced by n - nrow(A) p - 2*n-1 d - seq(-(n-1),(n-1)) yes, all good points, I have acted the changes you suggested, and it works smooth. thanks, Moreno -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
Hi Ben,
thank you soo much:) your code worked, and we
managed to halven computing time.
I am copying below the final version of the function,
if anybody is willing to use/improve it.
thanks again,
Moreno
###
spdiags = function(A){
require(Matrix)
source(getband.R)
# Find all nonzero diagonals
i = rep(seq(1, nrow(A),1),nrow(A));
j = sort(i);
d = sort(j-i);
d = unique(d);
p = length(d); ##the nr. col of the new matrix
m = nrow(A); n = ncol(A);
B = Matrix(FALSE, nrow = min(c(m,n)), ncol = p, sparse = TRUE);
count = 0
for (k in 1:p){
if (m = n){
i = max(1, 1+d[k]):min(n, m+d[k]);
} else { i = max(1, 1-d[k]):min(m,n-d[k]); }
if (length(i) != 0){
bd = getband(A, d[k]);
if (length(which(bd) == T) != 0){
count = count + 1;
# print( paste( non-zero diagonal, count) )
B[i,count] = bd
} else {
# print (paste (ncol, ncol(B), diag, k ) )
B = B[, -ncol(B)] }
}
}
return (list( B = B, d = d) )
}
Quoting Ben Bolker bbol...@gmail.com on Fri, 13 Apr 2012 20:09:52 +:
Ben Bolker bbolker at gmail.com writes:
I'm not quite sure how to do it, but I think you should look
at the ?band function in Matrix. In combination with diag() of a
suitably truncated matrix, you should be able to extract bands
of sparse matrices efficiently ...
getband - function(A,k) {
n - nrow(A)
if (abs(k)(n-1)) stop(bad band requested)
if (k0) {
v - seq(n-k) ## -seq((n-k+1),n)
w - seq(k+1,n) ## -seq(n-k-1)
} else if (k0) {
v - seq(-k+1,n)
w - seq(n+k)
} else return(diag(A))
diag(band(A,k,k)[v,w,drop=FALSE])
}
PS: I think this should extract the k^th off-diagonal
band in a way that should (?) work reasonably efficiently
with sparse matrices. I have not tested it carefully,
nor benchmarked it.
Ben Bolker
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
On 12-04-16 09:32 PM, Moreno I. Coco wrote:
Hi Ben,
thank you soo much:) your code worked, and we
managed to halven computing time.
I am copying below the final version of the function,
if anybody is willing to use/improve it.
thanks again,
Moreno
A couple of quick comments:
the semicolons at the end of each line are unnecessary (they're
harmless, but considered bad style -- most common in code of C and
MATLAB coders
You don't need to source(getband.R) inside the body of spdiags,
probably better to do it outside (so it doesn't get run every time the
function is called)
if (length(which(bd) == T) != 0)
can be shortened to
if (any(bd!=0))
I haven't looked too carefully but can't all the i/j/d/p stuff be
replaced by
n - nrow(A)
p - 2*n-1
d - seq(-(n-1),(n-1))
?
###
spdiags = function(A){
require(Matrix)
source(getband.R)
# Find all nonzero diagonals
i = rep(seq(1, nrow(A),1),nrow(A));
j = sort(i)
d = sort(j-i);
d = unique(d);
p = length(d); ##the nr. col of the new matrix
m = nrow(A); n = ncol(A);
B = Matrix(FALSE, nrow = min(c(m,n)), ncol = p, sparse = TRUE);
count = 0
for (k in 1:p){
if (m = n){
i = max(1, 1+d[k]):min(n, m+d[k]);
} else { i = max(1, 1-d[k]):min(m,n-d[k]); }
if (length(i) != 0){
bd = getband(A, d[k]);
if (length(which(bd) == T) != 0){
count = count + 1;
# print( paste( non-zero diagonal, count) )
B[i,count] = bd
} else {
# print (paste (ncol, ncol(B), diag, k ) )
B = B[, -ncol(B)] }
}
}
return (list( B = B, d = d) )
}
Quoting Ben Bolker bbol...@gmail.com on Fri, 13 Apr 2012 20:09:52 +:
Ben Bolker bbolker at gmail.com writes:
I'm not quite sure how to do it, but I think you should look
at the ?band function in Matrix. In combination with diag() of a
suitably truncated matrix, you should be able to extract bands
of sparse matrices efficiently ...
getband - function(A,k) {
n - nrow(A)
if (abs(k)(n-1)) stop(bad band requested)
if (k0) {
v - seq(n-k) ## -seq((n-k+1),n)
w - seq(k+1,n) ## -seq(n-k-1)
} else if (k0) {
v - seq(-k+1,n)
w - seq(n+k)
} else return(diag(A))
diag(band(A,k,k)[v,w,drop=FALSE])
}
PS: I think this should extract the k^th off-diagonal
band in a way that should (?) work reasonably efficiently
with sparse matrices. I have not tested it carefully,
nor benchmarked it.
Ben Bolker
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
Moreno I. Coco M.I.Coco at sms.ed.ac.uk writes: [snip snip snip] So, I have written my own spdiags function (below); following also a suggestion in an old, and perhaps unique post, about this issue. It works only for square matrices (that's my need), however I have a couple of issues, mainly related to computational efficiency: 1) if I use it with a sparseMatrix, it throws a tedious warning sparse[ logic ] : .M.sub.i.logical() maybe inefficient; can I suppress this warning somehow, this is slowing the computation very radically; quick answer: use suppressMessages() 2) I can go around this problem by translating a sparseMatrix back into a logical matrix before I run spdiags on it. However, the loop gets very slow for large matrices (e.g., 2000x2000), which is the kind of matrices I have to handle. If you look in the code, I have placed a system.time() where the code is slowing down, and it takes about: I'm not quite sure how to do it, but I think you should look at the ?band function in Matrix. In combination with diag() of a suitably truncated matrix, you should be able to extract bands of sparse matrices efficiently ... __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing spdiags function for R
Ben Bolker bbolker at gmail.com writes:
I'm not quite sure how to do it, but I think you should look
at the ?band function in Matrix. In combination with diag() of a
suitably truncated matrix, you should be able to extract bands
of sparse matrices efficiently ...
getband - function(A,k) {
n - nrow(A)
if (abs(k)(n-1)) stop(bad band requested)
if (k0) {
v - seq(n-k) ## -seq((n-k+1),n)
w - seq(k+1,n) ## -seq(n-k-1)
} else if (k0) {
v - seq(-k+1,n)
w - seq(n+k)
} else return(diag(A))
diag(band(A,k,k)[v,w,drop=FALSE])
}
PS: I think this should extract the k^th off-diagonal
band in a way that should (?) work reasonably efficiently
with sparse matrices. I have not tested it carefully,
nor benchmarked it.
Ben Bolker
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] writing spdiags function for R
Dear R-list,
I am in the process of translating a long function written in Matlab
into R (mainly because I am a big of fan of R, and folks will not
have to pay to use it :). In the translation of this function
I got stack because they use spdiags, which, as far as I can tell
it is not available in R. I have explored the Matrix package, from
which I borrowed some of the functions (e.g., sparseMatrix), but
I could not actually find an equivalent to spdiags (please, let
me know if it is there somewhere).
So, I have written my own spdiags function (below); following
also a suggestion in an old, and perhaps unique post, about
this issue.
It works only for square matrices (that's my need), however I
have a couple of issues, mainly related to computational
efficiency:
1) if I use it with a sparseMatrix, it throws a tedious warning
sparse[ logic ] : .M.sub.i.logical() maybe inefficient;
can I suppress this warning somehow, this is slowing the computation
very radically;
2) I can go around this problem by translating a sparseMatrix back
into a logical matrix before I run spdiags on it. However, the loop
gets very slow for large matrices (e.g., 2000x2000), which is the
kind of matrices I have to handle. If you look in the code,
I have placed a system.time() where the code is slowing down, and
it takes about:
user system elapsed
0.280.050.33
to complete an iteration...thus, I was wondering whether there is
a more efficient way to do what I am doing...also, if you spot
other places where the function could be optimized I would be
very glad to hear it!
thank you very much in advance for your kind help,
Best,
Moreno
###
## it works only for square matrices
## it could work with sparse matrices but it spits a tedious warning
## it is definitely inefficient compared to the original matlab code
## choose below different matrices to test the function.
# r = c(2,3,5,5); c = c(2,1,4,5)
# A = sparseMatrix(r, c)
# A = replicate(1000, rnorm(1000) )
# A = rbind(c(1,2,3),c(2,3,4),c(3,4,5))
spdiags = function(A){
# Find all nonzero diagonals
i = rep(seq(1, nrow(A),1),nrow(A));
j = sort(i);
d = sort(j-i);
# d = d(find(diff([-inf; d]))); ## from Matlab ...
# d = c(d[which(diff(d) == 1)], d[length(d)] ) ## this emulate
above but needs to stick in last element
d = unique(d); ##this should work just fine and it is simpler
p = length(d); ##the nr. col of the new matrix
m = nrow(A); n = ncol(A);
B = matrix(0, nrow = min(c(m,n)), ncol = p);
for (k in 1:p){
# print(k)
cl = vector();
if (m = n){
i = max(1, 1+d[k]):min(n, m+d[k]);
} else { i = max(1, 1-d[k]):min(m,n-d[k]); }
system.time(
if (length(i) != 0){
B[i,k] = A[ col(A) == row (A) - d[k]]
} )
}
return (list( B = B, d = d) )
}
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
