[Rd] Re: Method dispatch S3/S4 through optimize()

[Roger Bivand]

After looking at this with Roger Koenker off-list, the problem can be 
replicated running R CMD check tryout on:

http://reclus.nhh.no/R/Devel/tryout_0.1-1.tar.gz

With a NAMESPACE file in tryout, the S4 dispatch breaks down for det(), 
without a NAMESPACE file, it runs without error (on R-devel with SparseM 
0.52 for me). Roger pointed me to 

http://tolstoy.newcastle.edu.au/R/devel/04a/0042.html

which looks similar, but was said to be resolved. A practical but too 
invasive correction is to change the SparseM code to:

setGeneric("Det", function(x, ...) standardGeneric("Det"))
setMethod("Det","matrix",get("det", pos=NULL, mode= "function"))
setMethod("Det","matrix.csr", function(x, ...) Det(chol(x))^2)
setMethod("Det","matrix.csr.chol", function(x, ...) x at det)

rather than original below, the alternative is to drop NAMESPACE in 
spdep, which is not desirable either. 

Roger


On Fri, 19 Nov 2004, Roger Bivand wrote:

> On Thu, 18 Nov 2004, roger koenker wrote:
> 
> > I am far from an expert on S4 method dispatch, so at the risk of
> > embarrassing myself I will just indicate how the det() function is
> > written in SparseM now in the hope that someone else will see the source
> > of the problem:
> > 
> > setGeneric("det")
> > setMethod("det","matrix",get("det", pos=NULL, mode= "function"))
> > setMethod("det","matrix.csr", function(x, ...) det(chol(x))^2)
> > setMethod("det","matrix.csr.chol", function(x, ...) x at det)
> > 
> > The Cholesky function chol in SparseM was recently extended slightly to
> > compute this determinant as an additional component, so if the argument
> > is already of the matrix.csr.chol class det() just extracts this
> > component, otherwise it tries to do the cholesky, or defaults to the
> > base function.
> 
> As I read the SparseM code, it is following the specification, and if
> there was something wrong, the simple example 1) would not have worked
> either. But it does, so something changes in the context of where a
> function like det() and chol() is called from, here through optimize(),
> that sends it to the wrong place based on assuming that the class is 
> unknown. 
> 
> I'll try to make an example that doesn't depend on spdep, just on SparseM 
> and optimize().
> 
> This is frustrating, because with Ord's trick of using matrices similar to
> symmetric, the SparseM route will be very useful for most situations, 
> even though it requires strictly symmetric matrices. So your provision of 
> det() for "matrix.csr" is very useful, once we sort out this strange 
> behaviour.
> 
> Roger
> 
> > 
> > I hope that this might expose some fly in the soup.
> > 
> > 
> > url:	www.econ.uiuc.edu/~roger        	Roger Koenker
> > email	rkoenker at uiuc.edu			Department of Economics
> > vox: 	217-333-4558				University of Illinois
> > fax:   	217-244-6678				Champaign, IL 61820
> > 
> > On Nov 18, 2004, at 3:34 PM, Roger Bivand wrote:
> > 
> > > I have been running into difficulties with dispatching on an S4 class
> > > defined in the SparseM package, when the method calls are inside a
> > > function passed as the f= argument to optimize() in functions in the 
> > > spdep
> > > package. The S4 methods are typically defined as:
> > >
> > > setMethod("det","matrix.csr", function(x, ...) det(chol(x))^2)
> > >
> > > that is within setMethod() rather than by name before the setMethod().
> > >
> > > When called from within functions passed as the f= argument to 
> > > optimize,
> > > the S3 generics for det() and chol() get picked up, not the S4 generics
> > > for the S4 SparseM classes. This looks for instance like (from
> > > example(boston)):
> > >
> > >> gp2 <- lagsarlm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + 
> > >> I(RM^2)
> > > +  +  AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT),
> > > +  data=boston.c, nb2listw(boston.soi), method="SparseM")
> > > matrix.csr
> > > Error in chol(tmp1) : non-numeric argument to chol
> > >
> > > (this is the error message in chol.R in base). tmp1 is of class
> > > "matrix.csr", but is being sent to the S3 class (error in function
> > > sar.lag.mix.f.sM() in R/lag.spsarlm.R).
> > >
> > > sar.lag.mix.f.sM <- function(rho, W, e.a, e.b, e.c, n, quiet)
> > > {
> > > 	SSE <- e.a - 2*rho*e.b + rho*rho*e.c
> > > 	s2 <- SSE/n
> > > 	tmp1 <- asMatrixCsrIrW(W, rho)
> > > cat(class(tmp1), "\n")
> > > 	tmp2 <- chol(tmp1)
> > > #	tmp2 <- .cholMatrixCsr(tmp1)
> > > cat(class(tmp2), "\n")
> > > 	tmp3 <- (tmp2 at det)^2
> > > cat(tmp3, "\n")
> > > 	ret <- (log(tmp3)
> > > 		- ((n/2)*log(2*pi)) - (n/2)*log(s2) - (1/(2*s2))*SSE)
> > > 	if (!quiet)
> > > 	    cat("(SparseM) rho:\t", rho, "\tfunction value:\t", ret, "\n")
> > > 	ret
> > > }
> > >
> > >
> > > Three further observations:
> > >
> > > 1) In a simpler case:
> > >
> > >> library(spdep)
> > > Loading required package: tripack
> > > Loading required package: maptools
> > > Loading required package: foreign
> > > Loading required package: SparseM
> > > [1] "SparseM library loaded"
> > >> example(similar.listw)
> > > ...
> > > smlr.l> log(det(asMatrixCsrIrW(asMatrixCsrListw(similar.listw(COL.W)), 
> > > 0.5)))
> > > [1] -1.627660
> > >
> > > just works, and the appropriate det() and chol() are found; det()
> > > dispatches on class "matrix.csr" and finds the right chol() for that
> > > class - this is the same code context as found in the function passed 
> > > to
> > > optimize();
> > >
> > > 2) When the R source file containing the lagsarlm() function, and the
> > > function sent as f= through optimize is sourced(), the appropriate 
> > > det()
> > > and chol() are found; and
> > >
> > > 3) When the SparseM R code is modified, and the function defined prior 
> > > to
> > > its setMethod (so visible as a function independently of setMethod) and
> > > called by name rather than as a method, everything works, as it should.
> > > But this isn't a solution, setMethod() defined methods should dispatch 
> > > on
> > > class irrespective of context, I think.
> > >
> > > SparseM does not have a namespace, spdep does, but I don't think this 
> > > is
> > > the issue. I'm pretty sure this isn't an issue with SparseM either,
> > > because of 1).
> > >
> > > I've put a copy of the problematic spdep on:
> > >
> > > http://reclus.nhh.no/R/spdep/spdep_0.2-24.tar.gz
> > >
> > > should that be of any use, the issue is present in both R-2.0.1 and
> > > R-devel (2004-11-16), SparseM is 0.52.
> > >
> > > I hope that I'm missing something fairly obvious.
> > >
> > > Roger
> > >
> > > -- 
> > > Roger Bivand
> > > Economic Geography Section, Department of Economics, Norwegian School 
> > > of
> > > Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
> > > Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
> > > e-mail: Roger.Bivand at nhh.no
> > >
> > >
> > 
> > 
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no