From Roger.Bivand at nhh.no  Fri Nov 19 13:56:55 2004
From: Roger.Bivand at nhh.no (Roger Bivand)
Date: Fri Nov 19 13:57:08 2004
Subject: [Rd] Re: Method dispatch S3/S4 through optimize()
In-Reply-To: <D007130A-39C3-11D9-964C-000A95A7E3AA@uiuc.edu>
Message-ID: <Pine.LNX.4.44.0411191335550.17908-100000@reclus.nhh.no>

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@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@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@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@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@nhh.no


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