Re: [Rd] Compiling R with ATLAS

[Douglas Bates]

There is a Debian/Ubuntu specific list, r-sig-deb...@r-project.org,
which I am cc:'ing on this reply, that is a better location for this
discussion.

It appears that you may be going about things the hard way.  There are
Ubuntu packages for atlas and for R that can handle all of this for
you.  Take a look at http://cran.r-project.org/bin/linux/ubuntu for
information on how to get started.  It will probably be much easier to
install the pre-compiled packages first and worry about compilation
for customization later.

On Thu, May 13, 2010 at 11:16 PM, Avraham Adler  wrote:
> Hello. I know almost nothing about Linux, so I apologize if the answer to
> this is obvious.
>
> I am trying to compile R with ATLAS (successfully compiled 3.9.24 on Ubuntu
> 10.04 64 bit). I configured with:
>
> ./configure --enable-R-shlib --enable-BLAS-shlib
> --with-blas="-L/usr/local/atlas/lib -lf77blas -latlas"
> --with-lapack="-L/usr/local/atlaslib -llapack -lcblas" --with-x
>
> and I get this error (from config.log):
>
> configure:28566: checking for dgemm_ in -L/usr/local/atlas/lib -lf77blas
> -latlas
> configure:28587: gcc -std=gnu99 -o conftest -g -O2 -fpic
> -I/usr/local/include  -L/usr/local/lib64 conftest.c -L/usr/local/atlas/lib
> -lf77blas -latlas  -lgfortran -lm -ldl -lm  >&5
> conftest.c: In function 'main':
> conftest.c:193: warning: implicit declaration of function 'dgemm_'
> /usr/local/atlas/lib/libf77blas.so: undefined reference to
> `atl_f77wrap_zger2u_'
> /usr/local/atlas/lib/libf77blas.so: undefined reference to
> `atl_f77wrap_cger2c_'
> /usr/local/atlas/lib/libf77blas.so: undefined reference to
> `atl_f77wrap_zger2c_'
> /usr/local/atlas/lib/libf77blas.so: undefined reference to
> `atl_f77wrap_cger2u_'
> collect2: ld returned 1 exit status
> configure:28587: $? = 1
>
> I found the following e-mail in the archives <
> http://www.mail-archive.com/r-devel@r-project.org/msg16853.html>, but it has
> a different error, and the suggestion of switching out the BLAS later does
> not seem to apply to ATLAS which requires two files, I think.
>
> Any advice would be very appreciated.
>
> Thank you.
>
>        [[alternative HTML version deleted]]
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Failure to load the recommended package Matrix (Was: [R] Can one get a list of recommended packages?)

[Douglas Bates]

On Sun, Jun 13, 2010 at 8:28 AM, Uwe Ligges
 wrote:
>
>
> On 13.06.2010 01:09, Dr. David Kirkby wrote:
>>
>> On 06/12/10 05:27 PM, Douglas Bates wrote:
>>>
>>> On Sat, Jun 12, 2010 at 8:37 AM, Dr. David Kirkby
>>>  wrote:
>>>>
>>>> R 2.10.1 is used in the Sage maths project. Several recommended packages
>>>> (Matrix, class, mgcv, nnet, rpart, spatial, and survival) are failing to
>>>> build on Solaris 10 (SPARC).
>>>
>>> Have you checked the dependencies for those packages? Some require GNU
>>> make.
>>
>> We used GNU make.
>>
>>>> We would like to be able to get a list of the recommended packages for R
>>>> 2.10.1, but ideally via a call to R, so it is not necessary to update
>>>> that
>>>> list every time a new version of R is released. We do not want to
>>>> access the
>>>> Internet to get this information.
>>>
>>>> Is there a way in R to list the recommended packages?
>>>
>>> I'm not sure I understand the logic of this. If you are going to
>>> build R then presumably you have the tar.gz file which contains the
>>> sources for the recommended packages in the subdirectory
>>> src/library/Recommended/. Why not get the list from there?
>>
>> The reason is when the version of R gets updated in Sage, then someone
>> will have to check that list again, and more than likely fail to do so,
>> with the result tests will fail since packages do not exist, or worst
>> still we will be unaware they have failed to build properly.
>>
>> Therefore, being able to get them from a command would be useful, but
>> can understand if that is not possible.
>>
>>> $ cd ~/src/R-devel/src/library/Recommended/
>>> $ ls *.tgz
>>> boot.tgz codetools.tgz lattice.tgz mgcv.tgz rpart.tgz
>>> class.tgz foreign.tgz MASS.tgz nlme.tgz spatial.tgz
>>> cluster.tgz KernSmooth.tgz Matrix.tgz nnet.tgz survival.tgz
>>
>> OK, thank you for that list.
>>
>>>> Better still, is there a way to list the recommended packages which
>>>> have not
>>>> been installed, so getting a list of any failures?
>>>
>>> Again, this seems to be a rather convoluted approach. Why not check
>>> why the packages don't install properly?
>>
>> R had built, and the failure of the packages to build was not very
>> obvious, since it did not cause make to exit with a non-zero exit code.
>> Nobody had noticed until very recently that there was a problem.
>>
>> Therefore I proposed to make a test of the packages that should have
>> been installed, and ensure they actually all had.
>>
>> You need to be aware that R is just one part of Sage. Building the whole
>> of Sage takes a long time (>24 hours on some computers) so needless to
>> say, people will not view every line of error messages. The fact that
>> 'make' succeeded left us a false sense of security, when later it was
>> realsed there were problems when R run its self-tests.
>>
>> Dave
>
>
> But if you really want to sense some security, you should really run
> make check-all
> after the installation, particularly since you are on a platform that is not
> really mainstream any more.
>
> Uwe Ligges

And if you look at the other R-help message posted by David Kirby you
will find a link to the trouble ticket report in Sage as
http://trac.sagemath.org/sage_trac/ticket/9201
which, fortunately, contains a link to the compilation log file
http://sage.math.washington.edu/home/mpatel/trac/8306/r-2.10.1.p2.log

The log file correctly points out the problem, which is that the
compiled shared object for the Matrix package can't be loaded because
it contains references to objects in another library from the compiler
and that library is not on the LD_LIBRARY_PATH during the dyn.load
from within R.  So the original poster needs to look at how the
DLLpath is being configured.

Exactly how this ended up as a request to somehow determine a list of
recommended packages on the fly from within R is a wonderful example
of misguided speculation but has nothing whatsoever to do with the
problem.  Even the idea that there are a host of other recommended
packages that failed to compile because they depend on the Matrix
package is pure speculation and, of course, hopelessly wrong.  The
installation of the Matrix package failed because the check on whether
the package could be loaded failed.  The compilation of the
recommended packages stopped at that point.  The other packages
weren't compiled because they came after the Matrix package in the
compilation sequence.

Problems in compilation of R or its recommended packages should be
reported on the R-devel mailing list and I have moved this thread from
R-help to R-devel.

Admittedly it was the original poster of the Sage trouble ticket who
decided that the issue to focus on was missing recommended packages.
That's a red herring.  The issue is the visibility of the some of the
libraries from the gcc toolchain when R is loading a package.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Linking to lapack

[Douglas Bates]

On Fri, Oct 22, 2010 at 4:30 AM, Nick Sabbe  wrote:
> Hello all.

> I'm developing a package for R holding a Gibbs sampler, which tends to have
> better performance when written in C than in R.
> During each iteration in the Gibbs sampler, I need the inverse of a
> symmetric matrix.

You may want to check that.  The first rule of numerical linear
algebra is that you hardly ever need to invert a matrix.  If you need
to solve a linear system, you use a factorization, as below and then
use one of the solvers based on that factorization.  If the matrix is
small calculating the inverse is not a big problem.  If it is large
and you do so within each iteration of your sampler then you are
probably wasting time.

> For this, I wish to use lapack, as is concisely suggested in "Writing R
> extensions", since this will have better performance than I could ever write
> myself.
>
> After some twiddling I have got my code to compile by including
> "R_ext/Lapack.h" and using "F77_CALL(dpotrf)", but unfortunately, I don't
> get this to link properly.
> I get this message: " testc.o:testc.c:(.text+0x255): undefined reference to
> `dpotrf_'" which seems logical to me as far as my understanding of C
> reaches, but I don't know how to resolve it. I'm quite sure I need some
> extra parameters in my makefile, but as I come from a world where all these
> complexities are happily abstracted away for me by an IDE, I have no actual
> clue on how to surmount this.

Go back to section 1.2.1 of "Writing R Extensions" that deals with
"Using Makevars".

> However: when I'm done with all my code, I wish to build a package for
> publication on CRAN, so I want to be sure that not only I can build it on my
> system, but it will also work well when distributed to other computers (if I
> understand the package process well, source files are compiled and linked
> during installation of the package), so I would also like to know how to do
> this.
>
> It should not be relevant, but either way: I'm doing all this on a Windows 7
> machine, though the package will probably be used on Linux-based servers
> eventually.
>
> Finally: I have found no comprehensive list of the functions available to an
> R package developer, nor, strangely, questions about that. Does such a thing
> exist, or are we up to hoping we find what we are looking for in the header
> files? If it does not exist already, I would surely be willing to work on
> it.

Do you mean other than the section on "The R API" in the "Writing R
Extensions" manual?
> Thanks for any input.
>
> Nick Sabbe
> --
> ping: nick.sa...@ugent.be
> link: http://biomath.ugent.be
> wink: A1.056, Coupure Links 653, 9000 Gent
> ring: 09/264.59.36
>
> -- Do Not Disapprove
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [Rcpp-devel] must .Call C functions return SEXP?

[Douglas Bates]

On Thu, Oct 28, 2010 at 1:44 PM, Dominick Samperi  wrote:
> See comments on Rcpp below.
>
> On Thu, Oct 28, 2010 at 11:28 AM, William Dunlap  wrote:
>>
>> > -Original Message-
>> > From: r-devel-boun...@r-project.org
>> > [mailto:r-devel-boun...@r-project.org] On Behalf Of Andrew Piskorski
>> > Sent: Thursday, October 28, 2010 6:48 AM
>> > To: Simon Urbanek
>> > Cc: r-devel@r-project.org
>> > Subject: Re: [Rd] must .Call C functions return SEXP?
>> >
>> > On Thu, Oct 28, 2010 at 12:15:56AM -0400, Simon Urbanek wrote:
>> >
>> > > > Reason I ask, is I've written some R code which allocates two long
>> > > > lists, and then calls a C function with .Call.  My C code
>> > writes to
>> > > > those two pre-allocated lists,
>> >
>> > > That's bad! All arguments are essentially read-only so you should
>> > > never write into them!
>> >
>> > I don't see how.  (So, what am I missing?)  The R docs themselves
>> > state that the main point of using .Call rather than .C is that .Call
>> > does not do any extra copying and gives one direct access to the R
>> > objects.  (This is indeed very useful, e.g. to reorder a large matrix
>> > in seconds rather than hours.)
>> >
>> > I could allocate the two lists in my C code, but so far it was more
>> > convenient to so in R.  What possible difference in behavior can there
>> > be between the two approaches?
>>
>> Here is an example of how you break the rule that R-language functions
>> do not change their arguments if you use .Call in the way that you
>> describe.  The C code is in alter_argument.c:
>>
>> #include 
>> #include 
>>
>> SEXP alter_argument(SEXP arg)
>> {
>>    SEXP dim ;
>>    PROTECT(dim = allocVector(INTSXP, 2));
>>    INTEGER(dim)[0] = 1 ;
>>    INTEGER(dim)[1] = LENGTH(arg) ;
>>    setAttrib(arg, R_DimSymbol, dim);
>>    UNPROTECT(1) ;
>>    return dim ;
>> }
>>
>> Make a shared library out of this.  E.g., on Linux do
>>    R CMD SHLIB -o Ralter_argument.so alter_argument.so
>> and load it into R with
>>    dyn.open("./Ralter_argument.so")
>> (Or, on any platform, put it into a package along with
>> the following R code and build it.)
>>
>> The associated R code is
>>     myDim <- function(v).Call("alter_argument", v)
>>     f <- function(z) myDim(z)[2]
>> Now try using it:
>>     > myData <- 6:10
>>     > myData
>>     [1]  6  7  8  9 10
>>     > f(myData)
>>     [1] 5
>>     > myData
>>          [,1] [,2] [,3] [,4] [,5]
>>     [1,]    6    7    8    9   10
>> The argument to f was changed!  This should never happen in R.
>>
>> If you are very careful you might be able ensure that
>> no part of the argument to be altered can come from
>> outside the function calling .Call().  It can be tricky
>> to ensure that, especially when the argument is more complicated
>> than an atomic vector.
>>
>> "If you live outside the law you must be honest" - Bob Dylan.
>
> This thread seems to suggest (following Bob Dylan) that one needs
> to be very careful when using C/C++ to modify R's memory
> directly, because you may modify other R variables that point
> to the same memory (due to R's copy-by-value semantics and
> optimizations).
>
> What are the implications for the Rcpp package where R
> objects are exposed to the C++ side in precisely this way,
> permitting unrestricted modifications? (In the original
> or "classic" version of this package direct writes to R's
> memory were done only for performance reasons.)
>
> Seems like extra precautions need to be taken to
> avoid the aliasing problem.

The current Rcpp facilities has the same benefits and dangers as the C
macros used in .Call.  You get access to the memory of the R object
passed as an argument, saving a copy step.  You shouldn't modify that
memory.  If you do, bad things can happen and they will be your fault.
 If you want to get a read-write copy you clone the argument (in Rcpp
terminology).

To Bill:  I seem to remember the Dylan quote as "To live outside the
law you must be honest."


>
> Dominick
>
>> In R, .Call() does not copy its arguments but the C code
>> writer is expected to do so if they will be altered.
>> In S+ (and S), .Call() copies the arguments if altering
>> them would make a user-visible change in the environment,
>> unless you specify that the C code will not be altering them.
>>
>> Bill Dunlap
>> Spotfire, TIBCO Software
>> wdunlap tibco.com
>>
>> > > R has pass-by-value(!) semantics, so semantically you code has
>> > > nothing to do with the result.1 and result.2 variables since only
>> > > their *values* are guaranteed to be passed (possibly a copy).
>> >
>> > Clearly C code called from .Call must be allowed to construct R
>> > objects, as that's how much of R itself is implemented, and further
>> > down, it's what you recommend I should do instead.
>> >
>> > But why does it follow that C code must never modify an object
>> > initially allocated by R code?  Are you saying there is some special
>> > magic difference in the state of an object allocated by R's C code
>> > vs. one all

Re: [Rd] [R] DBLEPR?

[Douglas Bates]

Reply redirected to the R-devel list.

On Tue, Nov 16, 2010 at 1:02 PM, Prof. John C Nash  wrote:
> Ravi Varadhan and I have been looking at UCMINF to try to identify why it 
> gives occasional
> (but not reproducible) errors, seemingly on Windows only. There is some 
> suspicion that its
> use of DBLEPR for finessing the Fortran WRITE() statements may be to blame. 
> While I can
> find DBLEPR in Venables and Ripley, it doesn't get much mention after about 
> 2000 in the
> archives, though it is in the R FAQ and Brian R. mentions they are in libR in 
> a 2009 post.
> Are the xxxPR routines now deprecated (particularly for 64 bit systems) or 
> still OK to
> use?  If OK, can anyone point to documentation and examples?

Can you give more context?  What is UNCMINF?  Is it a package?  Is it a file?

There is a file R_SRC/src/appl/uncmin.c but it doesn't use DBLEPR
because it is C source code.  As stated in the comments

/* ../appl/uncmin.f
   -- translated by f2c (version of 1 June 1993 23:00:00).
   -- and hand edited by Saikat DebRoy
   */

/*--- The Dennis + Schnabel Minimizer -- used by R's  nlm() ---*/

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] DBLEPR?

[Douglas Bates]

On Tue, Nov 16, 2010 at 2:35 PM, Prof. John C Nash  wrote:
> I normally see digest once per day, but got msg from Doug Bates so responding 
> with context.

> UCMINF is a package on CRAN that implements a variable metric minimizer.

A pedant might point out that the package is called "ucminf".

> It does quite
> well on unconstrained problems. Stig Mortensen packaged the Fortran version 
> for R, but is
> not at moment responding to emails. There's also a Matlab version. We have it 
> in optimx
> and get some occasions where it just stops if we set trace>0. Other times we 
> can't get it
> to fail. My guess is something like an undefined variable or a bad 
> declaration of a
> pointer, but C and C++ are not languages I've much experience with.

Well, it doesn't work well for me because my version of gfortran (GNU
Fortran (Ubuntu/Linaro 4.4.4-14ubuntu5) 4.4.5) objects to the format
strings in some of the Fortran WRITE statements.

The recommended approach is to avoid all Fortran I/O including writing
to a Fortran character array.   As there are only 3 such WRITE
statements in the Fortran code it would be very simple to replace them
with calls to C functions that in turn call Rprintf.  However, it
would be best if Stig could take ownership of the modifications.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Competing with one's own work

[Douglas Bates]

On Fri, Dec 3, 2010 at 11:01 AM, Ravi Varadhan  wrote:
> Dear Duncan,

> What constitutes a convincing argument for making significant changes?
> Taking the example of optimization algorithms (say, for smooth objective
> functions), how does one make a convincing argument that a particular class
> of algorithms are "better" than another class? This can be a difficult task,
> but quite doable with good benchmarking practices.

> Supposing for the moment that such a convincing argument has been made, is
> that sufficient to get the R-core to act upon it?  Are there compelling
> factors other than just "algorithm A is better than algorithm B"?

> I'd think that the argument is relatively easy if the need for the change is
> driven by consumer demand. But, even here I am not sure how to make an
> argument to the R-core to consider the big changes.  For example, there is a
> reasonable demand for constrained (smooth) optimization algorithms in R
> (based on R-help queries).  Currently, there are only 3 packages that can
> handle this.  However, in the base distribution only `constrOptim' function
> is provided, which cannot handle anything more than linear, inequality
> constraints.  I think that the base distribution needs to have a package for
> constrained optimization that can handle linear/nonlinear and
> equality/inequality constraints.

constrOptim is in the stats package, not the base package.  Functions
that are already in the required packages are maintained by R core.
If you know of bugs in such functions you should report them.  Because
there is a heavy burden in maintaining the large corpus of software in
R and its required packages, additions are viewed skeptically,
Adopting new capabilities and new code in a required package like
stats means that some member of R core has to be willing to maintain
it.  If the capabilities can be incorporated in a contributed package
then that is the preferred method of extending R. The burden of
maintaining the code, fixing bugs or other infelicities, etc. is on
the package maintainer.

I don't see anything in what you are proposing that could not be
incorporated in a contributed package.

> John, thanks for raising an important issue.
>
> Thanks & Best,
> Ravi.
>
> ---
> Ravi Varadhan, Ph.D.
> Assistant Professor,
> Division of Geriatric Medicine and Gerontology School of Medicine Johns
> Hopkins University
>
> Ph. (410) 502-2619
> email: rvarad...@jhmi.edu
>
>
> -Original Message-
> From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org]
> On Behalf Of Duncan Murdoch
> Sent: Friday, December 03, 2010 11:13 AM
> To: nas...@uottawa.ca
> Cc: r-devel@r-project.org
> Subject: Re: [Rd] Competing with one's own work
>
> On 03/12/2010 10:57 AM, Prof. John C Nash wrote:
>> No, this is not about Rcpp, but a comment in that overly long discussion
> raised a question
>> that has been in my mind for a while.
>>
>> This is that one may have work that is used in R in the base functionality
> and there are
>> improvements that should be incorporated.
>>
>> For me, this concerns the BFGS, Nelder-Mead and CG options of optim(),
> which are based on
>> the 1990 edition (Pascal codes) of my 1979 book "Compact numerical
> methods...", which were
>> themselves derived from other people's work. By the time Brian Ripley took
> that work (with
>> permission, even though not strictly required. Thanks!) there were already
> some
>> improvements to these same algorithms (mainly bounds and masks) in the
> BASIC codes of the
>> 1987 book by Mary Walker-Smith and I. However, BASIC to R is not something
> I'd wish on
>> anyone.
>>
>> Now there are some R packages, including some I've been working on, that
> do offer
>> improvements on the optim() offerings. I would not say mine are yet fully
> ready for
>> incorporation into the base, but they are pretty close. Equally I think
> some of the tools
>> in the base should be deprecated and users encouraged to try other
> routines. It is also
>> getting more and more important that novice users be provided with
> sensible guidance and
>> robust default settings and choices. In many areas, users are faced with
> more choice than
>> is efficient for the majority of problems.
>>
>> My question is: How should such changes be suggested / assisted? It seems
> to me that this
>> is beyond a simple feature request. Some discussion on pros and cons would
> be appropriate,
>> and those like myself who are familiar with particular tools can and
> should offer help.
>>
>> Alternatively, is there a document available in the style "Writing R
> Extensions" that has
>> a title like "How the R Base Packages are Updated"? A brief search was
> negative.
>>
>> I'm happy to compete with my own prior work to provide improvements. It
> would be nice to
>> see some of those improvements become the benchmark for further progress.
>
>
> There are answers at many different levels to your question

Re: [Rd] Calling C++ from R

[Douglas Bates]

It is more effective to send such inquiries to the Rcpp-devel mailing
list which I am cc'ing on this reply.

On Thu, Jan 20, 2011 at 3:05 PM, mtck1982  wrote:
>
> Hi All,
>
> I am new to this area and use Rcpp to call C++ from R and try to build the
> package under Windows 7. I use Rtools and R 2.10.1 32bit. Everything works
> fine with me, except using R functions like "rnorm" or "runif" in the C++
> code. When I use "R CMD check" the package, it always return error
>
> ** libs
>  making DLL ...
> g++ -I"c:/PROGRA~2/R/R-210~1.1/include"
> -I"c:/PROGRA~2/R/R-210~1.1/library/Rcpp/include"     -O2 -Wall  -c func.cpp
> -o func.o
> func.cpp: In function 'SEXPREC* myfunction(SEXPREC*, SEXPREC*)':
> func.cpp:220: error: 'rgamma' was not declared in this scope
> func.cpp:225: error: 'rnorm' was not declared in this scope
> func.cpp:244: error: 'runif' was not declared in this scope
> func.cpp:274: error: 'rbeta' was not declared in this scope
> make: *** [func.o] Error 1
>  ... don

It is not clear if you are trying to use the Rcpp "sugar"
constructions, which can apply to entire vectors, of if you are using
the functions in the R API.  If the latter then you will need to
preface the name with Rf_, as in Rf_runif. Those are the actual names
of the functions.  The Rinternals.h file defines a number of aliases,
such as runif, but in Rcpp those aliases are turned off, so as to
avoid name clashes.

You should note that when calling random number generator functions
the programmer is responsible for getting and restoring the seed
structure.  I can't remember the details right now and I am on a
Windows system without the sources so I will rely on someone else to
fill in the details.

> although I already put ,  and  in the header file.
> //
> The func.cpp file has the following structure
>
> #include "func.h"
> SEXP myfunction(SEXP a, SEXP b) {
> }
> //
> The  header file "func.h" has the following content
>
> #include 
> #include 
> #include 
> #include 
> #include 
> RcppExport SEXP myfunction(SEXP a, SEXP b);
> /
> What could the problem be?
>
> Many thanks,
> Xiaochun
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Calling-C-from-R-tp3228490p3228490.html
> Sent from the R devel mailing list archive at Nabble.com.
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Problem using F77_CALL(dgemm) in a package

[Douglas Bates]

On Sun, Feb 20, 2011 at 4:39 PM, Jason Rudy  wrote:
> I've never used C++ before, so for this project I think I will stick
> with just using the BLAS and LAPACK routines directly.  Another issue
> is that I will need to do some sparse matrix computations, for which I
> am planning to use CSPARSE, at least to begin with.  I am interested
> by RcppArmadillo, and would consider it for future projects.  If you
> don't mind, what in your opinion are the major pros and cons of an
> RcppArmadillo solution compared to simply using the BLAS or LAPACK
> routines through the .C interface?

You may want to consider the API exported by the Matrix package that
allows access to CHOLMOD functions in that package's library.  The
entire CSPARSE library is also included in the Matrix package but most
of it is not exported because the CHOLMOD functions are generally more
effective.  (Both CHOLMOD and CSPARSE are written by Tim Davis.
CSPARSE is good code but it was written more for instructional purpose
than as an "industrial strength" package.)

> On Sun, Feb 20, 2011 at 8:11 AM, Dirk Eddelbuettel  wrote:
>>
>> On 20 February 2011 at 09:50, Dirk Eddelbuettel wrote:
>> | There is of course merit in working through the barebones API but in case 
>> you
>> | would consider a higher-level alternative, consider these few lines based 
>> on
>> | RcppArmadillo (which end up calling dgemm() for you via R's linkage to the 
>> BLAS)
>>
>> PS I always forget that we have direct support in Rcpp::as<> for Armadillo
>> matrices. The examples reduces to three lines in C++, and you never need to
>> worry about row or column dimension, or memory allocation or deallocation:
>>
>>  R> suppressMessages(library(inline))
>>  R> txt <- '
>>  +    arma::mat Am = Rcpp::as< arma::mat >(A);
>>  +    arma::mat Bm = Rcpp::as< arma::mat >(B);
>>  +    return Rcpp::wrap( Am * Bm );
>>  +    '
>>  R> mmult <- cxxfunction(signature(A="numeric", B="numeric"),
>>  +                      body=txt,
>>  +                      plugin="RcppArmadillo")
>>  R> A <- matrix(1:9, 3, 3)
>>  R> B <- matrix(9:1, 3, 3)
>>  R> C <- mmult(A, B)
>>  R> print(C)
>>       [,1] [,2] [,3]
>>  [1,]   90   54   18
>>  [2,]  114   69   24
>>  [3,]  138   84   30
>>  R>
>>
>> Matrices A and B from directly initialise Armadillo matrices, and the result
>> can be returned directly.
>>
>> Hth, Dirk
>>
>> --
>> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
>>
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] core Matrix package segfaulted on R CMD check --use-gct

[Douglas Bates]

Can you provide the output from

sessionInfo()

so we can know the platform?  Also, did you configure R with
--enable-strict-barrier or set the C compilation flag
-DTESTING_WRITE_BARRIER?  I think that run-time error message can only
be thrown under those circumstances (not that it isn't an error, it's
just not checked for in other circumstances).

On Sat, Mar 26, 2011 at 5:21 PM, Hin-Tak Leung  wrote:
> Current core/Recommended Matrix package (0.999375-48) has been segfaulting 
> against R 2.13-alpha/2.14-trunk for the last week or so (since R-2.13 was 
> branched, when I started trying) when "run with R CMD check --use-gct":
>
> --
>> pkgname <- "Matrix"
>> source(file.path(R.home("share"), "R", "examples-header.R"))
>> gctorture(TRUE)
>> options(warn = 1)
>> library('Matrix')
> Loading required package: lattice
> Error : .onLoad failed in loadNamespace() for 'Matrix', details:
>  call: fun(...)
>  error: unprotected object (0x2768b18) encountered (was REALSXP)
> Error: package/namespace load failed for 'Matrix'
> Execution halted
> ---
>
> I traced to this because "R CMD check --use-gct snpStats" (both 1.1.13 and 
> 1.1.12) segfaults with the same message, and before that, the snpMatrix 
> 1.15.8.4 which includes some of David's newly written ld() ( which depends on 
> Matrix.)
>
> If the Matrix package segfaults, David's new ld() isn't useable.
>
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] core Matrix package segfaulted on R CMD check --use-gct

[Douglas Bates]

On Tue, Mar 29, 2011 at 5:34 AM, Hin-Tak Leung
 wrote:
> Martin Maechler wrote:
>>>>>>>
>>>>>>> Douglas Bates 
>>>>>>>    on Mon, 28 Mar 2011 09:24:39 -0500 writes:
>>
>>    > Can you provide the output from sessionInfo()
>>
>>    > so we can know the platform?  Also, did you configure R
>>    > with --enable-strict-barrier or set the C compilation flag
>>    > -DTESTING_WRITE_BARRIER?  I think that run-time error
>>    > message can only be thrown under those circumstances (not
>>    > that it isn't an error, it's just not checked for in other
>>    > circumstances).
>>
>> interesting.
>>
>> In the mean time, I *did* run --- for several hours! ---
>> your code example below,
>> and it did *not* segfault for me (64-bit, Linux Fedora 13).
>>
>> Martin
>
> 64-bit fedora 14. For building R svn (and checking soon-to-be-released R
> packages, rather than daily R-related work), I also have these, and indeed
> have "--enable-strict-barrier":
>
> export DEFS='-DUSE_TYPE_CHECKING_STRICT -DR_MEMORY_PROFILING' \
> ./configure --enable-memory-profiling --enable-strict-barrier
> --enable-byte-compiled-packages --with-valgrind-instrumentation=2
>
>> sessionInfo()
> R version 2.14.0 Under development (unstable) (--)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8
>  [5] LC_MONETARY=C              LC_MESSAGES=en_GB.UTF-8
>  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base

Thanks for the information.  I can replicate the problem on a Red Hat
EL 5 64-bit system and will start to debug now.

>>>>>>> Douglas Bates 
>>>>>>>    on Mon, 28 Mar 2011 09:24:39 -0500 writes:
>>
>>    > Can you provide the output from
>>    > sessionInfo()
>>
>>    > so we can know the platform?  Also, did you configure R with
>>    > --enable-strict-barrier or set the C compilation flag
>>    > -DTESTING_WRITE_BARRIER?  I think that run-time error message can
>> only
>>    > be thrown under those circumstances (not that it isn't an error, it's
>>    > just not checked for in other circumstances).
>>
>>    > On Sat, Mar 26, 2011 at 5:21 PM, Hin-Tak Leung
>>  wrote:
>>    >> Current core/Recommended Matrix package (0.999375-48) has been
>> segfaulting against R 2.13-alpha/2.14-trunk for the last week or so (since
>> R-2.13 was branched, when I started trying) when "run with R CMD check
>> --use-gct":
>>    >>     >> --
>>    >>> pkgname <- "Matrix"
>>    >>> source(file.path(R.home("share"), "R", "examples-header.R"))
>>    >>> gctorture(TRUE)
>>    >>> options(warn = 1)
>>    >>> library('Matrix')
>>    >> Loading required package: lattice
>>    >> Error : .onLoad failed in loadNamespace() for 'Matrix', details:
>>    >>  call: fun(...)
>>    >>  error: unprotected object (0x2768b18) encountered (was REALSXP)
>>    >> Error: package/namespace load failed for 'Matrix'
>>    >> Execution halted
>>    >> ---
>>    >>     >> I traced to this because "R CMD check --use-gct snpStats"
>> (both 1.1.13 and 1.1.12) segfaults with the same message, and before that,
>> the snpMatrix 1.15.8.4 which includes some of David's newly written ld() (
>> which depends on Matrix.)
>>    >>     >> If the Matrix package segfaults, David's new ld() isn't
>> useable.
>>    >>     >>
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] core Matrix package segfaulted on R CMD check --use-gct

[Douglas Bates]

On Tue, Mar 29, 2011 at 2:17 PM, Douglas Bates  wrote:
> On Tue, Mar 29, 2011 at 5:34 AM, Hin-Tak Leung
>  wrote:
>> Martin Maechler wrote:
>>>>>>>>
>>>>>>>> Douglas Bates 
>>>>>>>>    on Mon, 28 Mar 2011 09:24:39 -0500 writes:
>>>
>>>    > Can you provide the output from sessionInfo()
>>>
>>>    > so we can know the platform?  Also, did you configure R
>>>    > with --enable-strict-barrier or set the C compilation flag
>>>    > -DTESTING_WRITE_BARRIER?  I think that run-time error
>>>    > message can only be thrown under those circumstances (not
>>>    > that it isn't an error, it's just not checked for in other
>>>    > circumstances).
>>>
>>> interesting.
>>>
>>> In the mean time, I *did* run --- for several hours! ---
>>> your code example below,
>>> and it did *not* segfault for me (64-bit, Linux Fedora 13).
>>>
>>> Martin
>>
>> 64-bit fedora 14. For building R svn (and checking soon-to-be-released R
>> packages, rather than daily R-related work), I also have these, and indeed
>> have "--enable-strict-barrier":
>>
>> export DEFS='-DUSE_TYPE_CHECKING_STRICT -DR_MEMORY_PROFILING' \
>> ./configure --enable-memory-profiling --enable-strict-barrier
>> --enable-byte-compiled-packages --with-valgrind-instrumentation=2
>>
>>> sessionInfo()
>> R version 2.14.0 Under development (unstable) (--)
>> Platform: x86_64-unknown-linux-gnu (64-bit)
>>
>> locale:
>>  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C
>>  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8
>>  [5] LC_MONETARY=C              LC_MESSAGES=en_GB.UTF-8
>>  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C
>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> Thanks for the information.  I can replicate the problem on a Red Hat
> EL 5 64-bit system and will start to debug now.

I isolated the problem and tested then committed a fix. I am going to
ask Martin to upload the new release as I have gotten out of sync with
some of his recent changes and he will, I hope, reconcile the branches
and trunk.  If you need the fixed version immediately, say for
testing, the changes are in the file Matrix/src/chm_common.c  You can
visit the SVN archive for the project,
https://r-forge.r-project.org/scm/?group_id=61, click on the link to
Browse Subversion Repository and go to pkg/Matrix/src/chm_common.c

I made the mistake that I chastised others for making, performing an
operation on the result of a call to install and another value that
may have required allocation of memory.  The first time install is
called on a particular string it allocates and SEXPREC, which may
cause a garbage collection.
>>>>>>>> Douglas Bates 
>>>>>>>>    on Mon, 28 Mar 2011 09:24:39 -0500 writes:
>>>
>>>    > Can you provide the output from
>>>    > sessionInfo()
>>>
>>>    > so we can know the platform?  Also, did you configure R with
>>>    > --enable-strict-barrier or set the C compilation flag
>>>    > -DTESTING_WRITE_BARRIER?  I think that run-time error message can
>>> only
>>>    > be thrown under those circumstances (not that it isn't an error, it's
>>>    > just not checked for in other circumstances).
>>>
>>>    > On Sat, Mar 26, 2011 at 5:21 PM, Hin-Tak Leung
>>>  wrote:
>>>    >> Current core/Recommended Matrix package (0.999375-48) has been
>>> segfaulting against R 2.13-alpha/2.14-trunk for the last week or so (since
>>> R-2.13 was branched, when I started trying) when "run with R CMD check
>>> --use-gct":
>>>    >>     >> --
>>>    >>> pkgname <- "Matrix"
>>>    >>> source(file.path(R.home("share"), "R", "examples-header.R"))
>>>    >>> gctorture(TRUE)
>>>    >>> options(warn = 1)
>>>    >>> library('Matrix')
>>>    >> Loading required package: lattice
>>>    >> Error : .onLoad failed in loadNamespace() for 'Matrix', details:
>>>    >>  call: fun(...)
>>>    >>  error: unprotected object (0x2768b18) encountered (was REALSXP)
>>>    >> Error: package/namespace load failed for 'Matrix'
>>>    >> Execution halted
>>>    >> ---
>>>    >>     >> I traced to this because "R CMD check --use-gct snpStats"
>>> (both 1.1.13 and 1.1.12) segfaults with the same message, and before that,
>>> the snpMatrix 1.15.8.4 which includes some of David's newly written ld() (
>>> which depends on Matrix.)
>>>    >>     >> If the Matrix package segfaults, David's new ld() isn't
>>> useable.
>>>    >>     >>
>>
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] C-Side: Applying a function (given as param) to data (given as param)

[Douglas Bates]

On Fri, Jun 3, 2011 at 5:17 AM, oliver  wrote:
> Hello,
>
> I'm implementing a package (C-extension),
> where one function gets data and a function
> that needs to be applied to the data.
>
> I want to apply the function to (parts of)
> the data on the C-side.
>
> 1) how do I apply a function (given via SEXP) to data
> 2) how do I select parts of the data (also provided via SEXP)?

Not to be facetious but you begin by reading the "Writing R Extensions" manual.

An alternative is to read the vignette Rcpp-Introduction available as
http://dirk.eddelbuettel.com/code/rcpp/Rcpp-introduction.pdf and soon
to be in The R Journal.  They show an explicit example of apply in
compiled code (C++ using the Rcpp structures, in their case).

> Comment on 2: I want to implement a moving/rolling/running
> apply on the data, so I want to select a data window.
> How can this be done, and how can it be done efficiently?
>
> For example on page 101 of "Writing R extensions" there was a trick
> on how to use REAL on the result, before using it in a loop.
> Can this be used for efficiency reasons in general?
>
>
>  Ciao,
>    Oliver
>
>
> P.S.: What do these macros (like REAL()) do? Do they convert data
>      or rather cast datatypes?
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Error using RcppGSL

[Douglas Bates]

Questions like this would get a faster response on the Rcpp-devel list, to
which I am copying this reply.
On Jun 21, 2011 6:35 AM, "oyvfos"  wrote:
> Hi, I get an error using RcppGSL: fatal error: gsl/gsl_vector.h:No such
file
> or directory. What is the best way to install these files as they seem to
> be missing?
> Thanks,
> Oyvind
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/Error-using-RcppGSL-tp3613535p3613535.html
> Sent from the R devel mailing list archive at Nabble.com.
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Performance of .C and .Call functions vs. native R code

[Douglas Bates]

On Thu, Jul 14, 2011 at 10:21 AM, Alireza Mahani
 wrote:
> (I am using a LINUX machine)
>
> Jeff,
>
> In creating reproducible results, I 'partially' answered my question. I have
> attached two scripts, 'mvMultiply.r' and 'mvMultiply.cc'. Please copy both
> files into your chosen directory, then run 'Rscript mvMultiply.r' in that
> directory while changing the two boolean parameters 'INCLUDE_DATAPREP' and
> 'ROWMAJOR' to all four permutations. (The variable 'diffVec' is there to
> verify that the two methods produce the same output vector.)
>
> Below are the results that I get, along with discussion (tR and tCall are in
> sec):
>
> INCLUDE_DATAPREP,ROWMAJOR,tR,tCall
> F,F,13.536,13.875
> F,T,13.824,14.299
> T,F,13.688,18.167
> T,T,13.982,30.730
>
> Interpretation: The execution time for the .Call line is nearly identical to
> the call to R operator '%*%'. Two data preparation lines for the .Call
> method create the overhead:
>
> A <- t(A) (~12sec, or 12usec per call)
> dim(A) <- dim(A)[1] * dim(A)[2] (~4sec, or 4usec per call)
>
> While the first line can be avoided by providing options in c++ function (as
> is done in the BLAS API), I wonder if the second line can be avoided, aside
> from the obvious option of rewriting the R scripts to use vectors instead of
> matrices. But this defies one of the primary advantages of using R, which is
> succinct, high-level coding. In particular, if one has several matrices as
> input into a .Call function, then the overhead from matrix-to-vector
> transformations can add up. To summarize, my questions are:
>
> 1- Do the above results seem reasonable to you? Is there a similar penalty
> in R's '%*%' operator for transforming matrices to vectors before calling
> BLAS functions?
> 2- Are there techniques for reducing the above overhead for developers
> looking to augment their R code with compiled code?
>
> Regards,
> Alireza
>
> ---
> # mvMultiply.r
> # comparing performance of matrix multiplication in R (using '%*%' operator)
> vs. calling compiled code (using .Call function)
> # y [m x 1] = A [m x n] %*% x [n x 1]
>
> rm(list = ls())
> system("R CMD SHLIB mvMultiply.cc")
> dyn.load("mvMultiply.so")
>
> INCLUDE_DATAPREP <- F
> ROWMAJOR <- F #indicates whether the c++ function treats A in a row-major or
> column-major fashion
>
> m <- 100
> n <- 10
> N <- 100
>
> diffVec <- array(0, dim = N)
> tR <- 0.0
> tCall <- 0.0
> for (i in 1:N) {
>
>        A <- runif(m*n); dim(A) <- c(m,n)
>        x <- runif(n)
>
>        t1 <- proc.time()[3]
>        y1 <- A %*% x
>        tR <- tR + proc.time()[3] - t1
>
>        if (INCLUDE_DATAPREP) {
>                t1 <- proc.time()[3]
>        }
>        if (ROWMAJOR) { #since R will convert matrix to vector in a 
> column-major
> fashion, if the c++ function expects a row-major format, we need to
> transpose A before converting it to a vector
>                A <- t(A)
>        }
>        dim(A) <- dim(A)[1] * dim(A)[2]
>        if (!INCLUDE_DATAPREP) {
>                t1 <- proc.time()[3]
>        }
>        y2 <- .Call("matvecMultiply", as.double(A), as.double(x),
> as.logical(c(ROWMAJOR)))
>        tCall <- tCall + proc.time()[3] - t1
>
>        diffVec[i] <- max(abs(y2 - y1))
> }
> cat("Data prep time for '.Call' included: ", INCLUDE_DATAPREP, "\n")
> cat("C++ function expects row-major matrix: ", ROWMAJOR, "\n")
> cat("Time - Using '%*%' operator in R: ", tR, "sec\n")
> cat("Time - Using '.Call' function: ", tCall, "sec\n")
> cat("Maximum difference between methods: ", max(diffVec), "\n")
>
> dyn.unload("mvMultiply.so")
> ---
> # mvMultiply.cc
> #include 
> #include 
>
> extern "C" {
>
> SEXP matvecMultiply(SEXP A, SEXP x, SEXP rowmajor) {
>        double *rA = REAL(A), *rx = REAL(x), *ry;
>        int *rrm = LOGICAL(rowmajor);
>        int n = length(x);
>        int m = length(A) / n;
>        SEXP y;
>        PROTECT(y = allocVector(REALSXP, m));
>        ry = REAL(y);
>        for (int i = 0; i < m; i++) {
>                ry[i] = 0.0;
>                for (int j = 0; j < n; j++) {
>                        if (rrm[0] == 1) {
>                                ry[i] += rA[i * n + j] * rx[j];
>                        } else {
>                                ry[i] += rA[j * m + i] * rx[j];
>                        }
>                }
>        }
>        UNPROTECT(1);
>        return(y);
> }
>
> }
>

I realize that you are just beginning to use the .C and .Call
interfaces and your example is therefore a simple one.  However, if
you plan to continue with such development it is worthwhile learning
of some of the tools available.  I think one of the most important is
the "inline" package that can take a C or C++ code segment as a text
string and go through all the steps of creating and loading a
.Call'able compiled function.

Second, if you are going to use C++ (the code you show could be C code
as it doesn't use any C++ extensions) then you should look at th

Re: [Rd] Performance of .C and .Call functions vs. native R code

[Douglas Bates]

I just saw that I left a syntax error in the .R and the first
_Rout.txt files.  Notice that in the second _Rout.txt file the order
of the arguments in the constructors for the MMatrixXd and the
MVectorXd are in a different order than in the .R and the first
_Rout.txt files.  The correct order has the pointer first, then the
dimensions.  For the first _Rout.txt file this part of the code is not
used.

On Tue, Jul 19, 2011 at 10:00 AM, Douglas Bates  wrote:
> On Thu, Jul 14, 2011 at 10:21 AM, Alireza Mahani
>  wrote:
>> (I am using a LINUX machine)
>>
>> Jeff,
>>
>> In creating reproducible results, I 'partially' answered my question. I have
>> attached two scripts, 'mvMultiply.r' and 'mvMultiply.cc'. Please copy both
>> files into your chosen directory, then run 'Rscript mvMultiply.r' in that
>> directory while changing the two boolean parameters 'INCLUDE_DATAPREP' and
>> 'ROWMAJOR' to all four permutations. (The variable 'diffVec' is there to
>> verify that the two methods produce the same output vector.)
>>
>> Below are the results that I get, along with discussion (tR and tCall are in
>> sec):
>>
>> INCLUDE_DATAPREP,ROWMAJOR,tR,tCall
>> F,F,13.536,13.875
>> F,T,13.824,14.299
>> T,F,13.688,18.167
>> T,T,13.982,30.730
>>
>> Interpretation: The execution time for the .Call line is nearly identical to
>> the call to R operator '%*%'. Two data preparation lines for the .Call
>> method create the overhead:
>>
>> A <- t(A) (~12sec, or 12usec per call)
>> dim(A) <- dim(A)[1] * dim(A)[2] (~4sec, or 4usec per call)
>>
>> While the first line can be avoided by providing options in c++ function (as
>> is done in the BLAS API), I wonder if the second line can be avoided, aside
>> from the obvious option of rewriting the R scripts to use vectors instead of
>> matrices. But this defies one of the primary advantages of using R, which is
>> succinct, high-level coding. In particular, if one has several matrices as
>> input into a .Call function, then the overhead from matrix-to-vector
>> transformations can add up. To summarize, my questions are:
>>
>> 1- Do the above results seem reasonable to you? Is there a similar penalty
>> in R's '%*%' operator for transforming matrices to vectors before calling
>> BLAS functions?
>> 2- Are there techniques for reducing the above overhead for developers
>> looking to augment their R code with compiled code?
>>
>> Regards,
>> Alireza
>>
>> ---
>> # mvMultiply.r
>> # comparing performance of matrix multiplication in R (using '%*%' operator)
>> vs. calling compiled code (using .Call function)
>> # y [m x 1] = A [m x n] %*% x [n x 1]
>>
>> rm(list = ls())
>> system("R CMD SHLIB mvMultiply.cc")
>> dyn.load("mvMultiply.so")
>>
>> INCLUDE_DATAPREP <- F
>> ROWMAJOR <- F #indicates whether the c++ function treats A in a row-major or
>> column-major fashion
>>
>> m <- 100
>> n <- 10
>> N <- 100
>>
>> diffVec <- array(0, dim = N)
>> tR <- 0.0
>> tCall <- 0.0
>> for (i in 1:N) {
>>
>>        A <- runif(m*n); dim(A) <- c(m,n)
>>        x <- runif(n)
>>
>>        t1 <- proc.time()[3]
>>        y1 <- A %*% x
>>        tR <- tR + proc.time()[3] - t1
>>
>>        if (INCLUDE_DATAPREP) {
>>                t1 <- proc.time()[3]
>>        }
>>        if (ROWMAJOR) { #since R will convert matrix to vector in a 
>> column-major
>> fashion, if the c++ function expects a row-major format, we need to
>> transpose A before converting it to a vector
>>                A <- t(A)
>>        }
>>        dim(A) <- dim(A)[1] * dim(A)[2]
>>        if (!INCLUDE_DATAPREP) {
>>                t1 <- proc.time()[3]
>>        }
>>        y2 <- .Call("matvecMultiply", as.double(A), as.double(x),
>> as.logical(c(ROWMAJOR)))
>>        tCall <- tCall + proc.time()[3] - t1
>>
>>        diffVec[i] <- max(abs(y2 - y1))
>> }
>> cat("Data prep time for '.Call' included: ", INCLUDE_DATAPREP, "\n")
>> cat("C++ function expects row-major matrix: ", ROWMAJOR, "\n")
>> cat("Time - Using '%*%' operator in R: ", tR, "sec\n")
>> cat("Time - Using '.Call' function: ", tCall, "sec\n")
>> cat("Ma

Re: [Rd] Use of Matrix within packages in R-2.14.0

[Douglas Bates]

On Nov 11, 2011 6:13 PM, "Bryan W. Lewis"  wrote:
">
> Dear R-devel readers:
>
>  I am really stuck trying resolving an issue with the use of the
> Matrix in one of my packages, irlba, with R-2.14.0. When I use
> crossprod with at least one sparse argument in the packaged code I
> receive the error:
>
> Error in crossprod(x, y) :
>  requires numeric/complex matrix/vector arguments
You must create a NAMESPACE file and import the crossprod function from the
Matrix package.

> However, when I run the code outside of the context of the package it
> works fine.
>
> I have constructed the following trivial dummy package to replicate
> the problem, the contents of which follow:
>
> DESCRIPTION file:
> Package: dummy
> Type: Package
> Title: Dummy package for testing
> Version: 0.0.0
> Date: 2011-10-15
> Author: B. W. Lewis 
> Maintainer: B. W. Lewis 
> Description: A bogus test package
> Depends: Matrix
> License: GPL
>
> NAMESPACE file:
> export("test")
>
> R/test.R file:
> `test` <- function(A)
> {
>  x = crossprod(A)
>  return(0)
> }
>
>
> To replicate the problem, create a package called "dummy" from the
> above three files, install it, and run from R-2.14.0:
>
> library("dummy")
> A = Matrix(0,10,10)
> test(A)
>
> I have tested this with Ubuntu 10.10 GNU/Linux (64-bit), and 32-bit
> Windows versions of R-2.14.0.
>
> Any help in resolving this issue is greatly appreciated!
>
> Thanks,
>
> Bryan
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] RcppArmadillo compilation error: R CMD SHLIB returns status 1

[Douglas Bates]

On Dec 6, 2011 8:30 AM, "Duncan Murdoch"  wrote:
>
> On 05/12/2011 1:22 PM, Paul Viefers wrote:
>>
>> Dear all,
>>
>> running the example by D. Eddebuettel (
http://dirk.eddelbuettel.com/blog/2011/04/23/) I get an error message.
Specifically, the R code I was taking from the above example is
>>
>> ### BEGIN EXAMPLE ###
>>
>> suppressMessages(require(RcppArmadillo))
>> suppressMessages(require(Rcpp))
>> suppressMessages(require(inline))
>> code<- '
>>arma::mat coeff = Rcpp::as(a);
>>arma::mat errors = Rcpp::as(e);
>>int m = errors.n_rows; int n = errors.n_cols;
>>arma::mat simdata(m,n);
>>simdata.row(0) = arma::zeros(1,n);
>>for (int row=1; row>  simdata.row(row) = simdata.row(row-1)*trans(coeff)+errors.row(row);
>>}
>>return Rcpp::wrap(simdata);
>>  '
>> ## create the compiled function
>> rcppSim<- cxxfunction(signature(a="numeric",e="numeric"),
>> code,plugin="RcppArmadillo")
>>
>> ### END OF EXAMPLE ###
>>
>> Executing this inside R, returned the following:
>>
>> ERROR(s) during compilation: source code errors or compiler
configuration errors!
>>
>> Program source:
>>   1:
>>   2: // includes from the plugin
>>   3: #include
>>   4: #include
>>   5:
>>   6:
>>   7: #ifndef BEGIN_RCPP
>>   8: #define BEGIN_RCPP
>>   9: #endif
>>  10:
>>  11: #ifndef END_RCPP
>>  12: #define END_RCPP
>>  13: #endif
>>  14:
>>  15: using namespace Rcpp;
>>  16:
>>  17:
>>  18: // user includes
>>  19:
>>  20:
>>  21: // declarations
>>  22: extern "C" {
>>  23: SEXP file33765791( SEXP a, SEXP e) ;
>>  24: }
>>  25:
>>  26: // definition
>>  27:
>>  28: SEXP file33765791( SEXP a, SEXP e ){
>>  29: BEGIN_RCPP
>>  30:
>>  31:arma::mat coeff = Rcpp::as(a);
>>  32:arma::mat errors = Rcpp::as(e);
>>  33:int m = errors.n_rows; int n = errors.n_cols;
>>  34:arma::mat simdata(m,n);
>>  35:simdata.row(0) = arma::zeros(1,n);
>>  36:for (int row=1; row>  37:  simdata.row(row) =
simdata.row(row-1)*trans(coeff)+errors.row(row);
>>  38:}
>>  39:return Rcpp::wrap(simdata);
>>  40:
>>  41: END_RCPP
>>  42: }
>>  43:
>>  44:
>> Error in compileCode(f, code, language = language, verbose = verbose) :
>>   Compilation ERROR, function(s)/method(s) not created!
>> Executing command 'C:/PROGRA~1/R/R-214~1.0/bin/i386/R CMD SHLIB
file33765791.cpp 2>  file33765791.cpp.err.txt' returned status 1
>>
>> I am working under R 2.14.0 and as the pros among you might guess, I am
new to using the C++ interfaces within R. I think all I have to do is to
edit some settings on my Windows 7 machine here, but the error message is
too cryptic to me. Alas, I could also not find any thread or help topic
that deals with this online. I appreciate any direct reply or reference
where I can find a solution to this.
>> Please let me know in case I am leaving out some essential details here.
>
>
> If you put the program source into a file (e.g. fn.cpp) and in a Windows
cmd shell you run
>
> R CMD SHLIB fn.cpp
>
> what do you get?   I would guess you've got a problem with your setup of
the compiler or other tools, and this would likely show it.

I don't think that will work because you need the appropriate -I option to
get the headers from the RcppArmadillo package.  It may be easier to use
the RcppArmadillo.package.skeleton function to create a package.

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Cholesky update/downdate

[Douglas Bates]

On Thu, Dec 29, 2011 at 8:05 AM, Yves Deville  wrote:
> Dear R-devel members,

> I am looking for a fast Cholesky update/downdate. The matrix A being
> symmetric positive definite (n, n) and factorized as
> A = L %*% t(L), the goal is to factor the new matrix  A +- C %*% t(C) where
> C is (n, r). For instance, C is 1-column when adding/removing an observation
> in a linear regression. Of special interest is the case where A is sparse.

> Looking at the 'Matrix' package (help and source code), it seems that the
> CHOLMOD library shipped with 'Matrix' allows this,
> but is not (yet?) interfaced in 'Matrix', where the 'update' method for
> Cholesky decomposition objects seems limited to a new matrix A + m*I with a
> scalar (diagonal) modification.

The CHOLMOD library provides sparse matrix methods, especially the
Cholesky decomposition and modifications to that decomposition, which
is where the name comes from.  Do you expect to work with sparse
matrices?

I haven't seem too much code for update/downdate operations on the
Cholesky decomposition for dense matrices.  There were rank-1
update/downdate methods in Linpack but they didn't make it through to
Lapack.
> If this is true: are there plans to implement such up/downdates?
>
> Best,
>
> Yves
>
> Yves Deville, statistical consultant, France.
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Copying objects prior to .Call

[Douglas Bates]

On Wed, Jan 11, 2012 at 11:49 AM, Simon Urbanek
 wrote:
>
> On Jan 11, 2012, at 12:08 PM, Taylor Arnold wrote:
>
>> R-devel,
>>
>> I have noticed that making a copy of an object in R prior to using
>> .Call on the original object can
>> cause the C code to alter not only the object passed to it but also
>> the copy in R.

> Please see the docs - .Call does *NOT* have a DUP argument - you are 
> responsible for duplication at all times if you make modifications (e.g. 
> using duplicate()).

Except that duplicate will create a new SEXPREC and the original
poster wanted to modify the SEXPREC passed through a pointer in .Call.

Purposely changing the value of arguments to .Call is a bad design,
Taylor.  R is a functional language and this breaks the functional
semantics. It is the sort of thing that you do only when you can't
think of a better approach.  It may, perhaps, be justified if the R
objects you are passing happen to be fields in a reference class (see
?setRefClass) but otherwise it is just opening yourself up to errors.
At the level of C/C++ all R objects passed as arguments to .Call
should be regarded as

const SEXP


>> A simple example
>> is:
>>
>>> x <- 2
>>> y <- x
>>> .Call("addOne", x, DUP=TRUE) # Changing DUP does not alter output
>> NULL
>>> x
>> [1] 3
>>> y
>> [1] 3
>>
>> And corresponding simple C code:
>>
>> "test.c":
>> #include 
>> #include 
>> #include 
>>
>> SEXP addOne(SEXP input) {
>>   REAL(input)[0] = REAL(input)[0] + 1;
>>   return R_NilValue;
>> }
>>
>> I assume that this is simply a result of lazy loading in R, and well
>> documented. My question is, do
>> there exist functions to (1) force R to make a copy of an object
>> (force() does not work), and (2) to check
>> whether two objects are actually pointing to the same memory address.
>> For question 1, I have
>> found specific operations which force a copy of a given datatype, but
>> would prefer a more general
>> solution.
>>
>> Thank you,
>>
>> Taylor
>>
>>> sessionInfo()
>> R version 2.14.1 (2011-12-22)
>> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>>
>> locale:
>> [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> loaded via a namespace (and not attached):
>> [1] tools_2.14.1
>>
>>
>> --
>> Taylor B. Arnold
>> Department of Statistics
>> Yale University
>> 24 Hillhouse Avenue
>> New Haven, CT 06520
>>
>> e-mail: taylor.arn...@yale.edu
>>
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>>
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Julia

[Douglas Bates]

My purpose in mentioning the Julia language (julialang.org) here is
not to start a flame war.  I find it to be a very interesting
development and others who read this list may want to read about it
too.

It is still very much early days for this language - about the same
stage as R was in 1995 or 1996 when only a few people knew about it -
but Julia holds much potential.  There is a thread about "R and
statistical programming" on groups.google.com/group/julia-dev.  As
always happens, there is a certain amount of grumbling of the "R IS
S SLW" flavor but there is also some good discussion regarding
features of R (well, S actually) that are central to the language.
(Disclaimer: I am one of the participants discussing the importance of
data frames and formulas in R.)

If you want to know why Julia has attracted a lot of interest very
recently (like in the last 10 days), as a language it uses multiple
dispatch (like S4 methods) with methods being compiled on the fly
using the LLVM (http://llvm.org) infrastructure.  In some ways it
achieves the Holy Grail of languages like R, Matlab, NumPy, ... in
that it combines the speed of compiled languages with the flexibility
of the high-level interpreted language.

One of the developers, Jeff Bezanson, gave a seminar about the design
of the language at Stanford yesterday, and the video is archived at
http://www.stanford.edu/class/ee380/.  You don't see John Chambers on
camera but I am reasonably certain that a couple of the questions and
comments came from him.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Julia

[Douglas Bates]

On Thu, Mar 1, 2012 at 11:20 AM, Jeffrey Ryan  wrote:
> Doug,
>
> Agreed on the interesting point - looks like it has some real promise.
>  I think the spike in interest could be attributable to Mike
> Loukides's tweet on Feb 20. (editor at O'Reilly)
>
> https://twitter.com/#!/mikeloukides/status/171773229407551488
>
> That is exactly the moment I stumbled upon it.

I think Jeff Bezanson attributes the interest to a blog posting by
Viral Shah, another member of the development team, that hit Reddit.
He said that, with Viral now in India, it all happened overnight for
those in North America and he awoke the next day to find a firestorm
of interest.  I ran across Julia in the Release Notes of LLVM and
mentioned it to Dirk Eddelbuettel who posted about it on Google+ in
January.  (Dirk, being much younger than I, knows about these
new-fangled social media things and I don't.)

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Small changes to big objects (1)

[Douglas Bates]

Is there a difference in the copying behavior of

x@little <- other

and

x@little[] <- other

I was using the second form in (yet another!) modification of the internal
representation of mixed-effects models in the lme4 package in the hopes
that it would not trigger copying of the entire object.  The object
representing the model is quite large but the changes during iterations are
to small vectors representing parameters and coefficients.



On Thu, Jan 3, 2013 at 1:08 PM, John Chambers  wrote:

> Martin Morgan commented in email to me that a change to any slot of an
> object that has other, large slot(s) does substantial computation,
> presumably from copying the whole object.  Is there anything to be done?
>
> There are in fact two possible changes, one automatic but only partial,
> the other requiring some action on the programmer's part.  Herewith the
> first; I'll discuss the second in a later email.
>
> Some context:  The notion is that our object has some big data and some
> additional smaller things.  We need to change the small things but would
> rather not copy the big things all the time.  (With long vectors, this
> becomes even more relevant.)
>
> There are three likely scenarios: slots, attributes and named list
> components.  Suppose our object has "little" and "BIG" encoded in one of
> these.
>
> The three relevant computations are:
>
> x@little <- other
> attr(x, "little") <- other
> x$little <- other
>
> It turns out that these are all similar in behavior with one important
> exception--fixing that is the automatic change.
>
> I need to review what R does here. All these are replacement functions,
> `@<-`, `attr<-`, `$<-`.  The evaluator checks before calling any
> replacement whether the object needs to be duplicated (in a routine
> EnsureLocal()).  It does that by examining a special field that holds the
> reference status of the object.
>
> Some languages, such as Python (and S) keep reference counts for each
> object, de-allocating the object when the reference count drops back to
> zero.  R uses a different strategy. Its NAMED() field is 0, 1 or 2
> according to whether the object has been assigned never, once or more than
> once.  The field is not a reference count and is not decremented--relevant
> for this issue.  Objects are de-allocated only when garbage collection
> occurs and the object does not appear in any current frame or other context.
> (I did not write any of this code, so apologies if I'm misrepresenting it.)
>
> When any of these replacement operations first occurs for a particular
> object in a particular function call, it's very likely that the reference
> status will be 2 and EnsureLocal will duplicate it--all of it. Regardless
> of which of the three forms is used.
>
> Here the non-level-playing-field aspect comes in.  `@<-` is a normal R
> function (a "closure") but the other two are primitives in the main code
> for R.  Primitives have no frame in which arguments are stored.  As a
> result the new version of x is normally stored with status 1.
>
> If one does a second replacement in the same call (in a loop, e.g.) that
> should not normally copy again.  But the result of `@<-` will be an object
> from its frame and will have status 2 when saved, forcing a copy each time.
>
> So the change, naturally, is that R 3.0.0 will have a primitive
> implementation of `@<`.  This has been implemented in r-devel (rev. 61544).
>
> Please try it out _before_ we issue that version, especially if you own a
> package that does things related to this question.
>
> John
>
> PS:  Some may have noticed that I didn't mention a fourth approach: fields
> in a reference class object.  The assumption was that we wanted classical,
> functional behavior here.  Reference classes don't have the copy problem
> but don't behave functionally either.  But that is in fact the direction
> for the other approach.  I'll discuss that later, when the corresponding
> code is available.
>
> __**
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/**listinfo/r-devel
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] alpha, portable use

[Douglas Bates]

On 9/25/06, Paul Gilbert <[EMAIL PROTECTED]> wrote:
> I am still confused about this (and it is still happening with R-beta).
> Writing R Extensions suggests I need a Makefile or Makevars in my
> package,  but that has not been the case previously.  Is there now a
> requirement that all packages need Makefiles or Makevars if there is
> fortran to be compiled?  This is only happening on one of my systems.
> Building R and make check work fine on that system, but it seems that
> not all the information gets passed along to package compiles.
>
> (BTW, this is just a warning, but Kurt suggested we try to eliminate
> warnings.)

>From the warning it seems that you have a Makevars file in the src
directory for your package.  The change in R-2.4.0 is that packages
that did have a Makevars file of the form

PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS)

should change that to

PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

>
> Paul Gilbert
>
> Prof Brian Ripley wrote:
> > On Wed, 20 Sep 2006, Paul Gilbert wrote:
> >
> >> When I build one of my packages with alpha from yesterday I am getting
> >>
> >> * checking for portable use of $BLAS_LIBS ... WARNING
> >> apparently missing $(FLIBS) in 'PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS)'
> >>
> >> Is this something I should worry about?  (Possibly I got this before and
> >> didn't notice.)
> >
> > Yes, please do check Writing R Extensions 
> >
> 
>
> La version française suit le texte anglais.
>
> 
>
> This email may contain privileged and/or confidential inform...{{dropped}}
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] x86_64, acml-3.5.0-gfortran64 and lme4

[Douglas Bates]

I am encountering segfaults when checking the lme4 package on an
Athlon64 system if I use the acml blas.  R was built as a 64-bit
application using the GCC 4.0.3 compiler suite including gfortran.
The version of acml is 3.5.0 gfortran64.

I do not encounter the segfaults when I compile R with R's built-in
BLAS.  The errors occur in the first example in lme4 in a call to
lmer. It looks like they would occur in any call to lmer.  Running
under the debugger shows that the segfault occurs in a call to dtrsm
(a level-3 BLAS routine to solve a triangular system of equations)
that is called from within a cholmod (sparse matrix library) routine.

Has anyone succeeded in running R CMD check on the lme4 package with
accelerated BLAS?  I'm trying to pin down is this occurs only with
ACML or also with Atlas and/or Goto's BLAS.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] error in lme4 for R 2.4.0 (PR#9297)

[Douglas Bates]

On 10/14/06, [EMAIL PROTECTED] <[EMAIL PROTECTED]> wrote:
> Full_Name: Din Chen
> Version: 2.4.0
> OS: Windows XP
> Submission from: (NULL) (66.17.122.18)
>
>
> I just updated the R.2.4.0. and got the error message for random effect
> modelling, which was working on R.2.3.1.
>
> library(lme4)
> mmod <- lmer(bright ~ 1+(1|operator), pulp)
> summary(mmod)
>
> Then when I tried to extract the residuals and random effect using:
> resid(mmodr)
> ranef(mmodr)
>
> I got error message:
>
> Error in resid(mmod) : no slot of name "family" for this object of class 
> "lmer"

Insufficient testing by the author.  I need to modify the resid
methods for the new class definitions.  I'll upload lme4_0.9975-6 with
a fix for this later today.

The ranef method does work ok, I believe.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] lme4 install, was Re: [R] Which executable is associated with R CMD INSTALL?

[Douglas Bates]

On 31 Oct 2006 12:05:21 +0100, Peter Dalgaard <[EMAIL PROTECTED]> wrote:
>
> [move to r-devel, put maintainer in loop]
>
> Patrick Connolly <[EMAIL PROTECTED]> writes:
>
> > On Mon, 30-Oct-2006 at 04:44PM -0500, Duncan Murdoch wrote:
> >
> >
> > |> Try "R CMD printenv R_HOME" and you'll find which R home directory it is
> > |> using.  You can see a lot more with "R CMD printenv" or various options
> > |> to "R CMD config".
> >
> > Thanks for that information.  It knocks my theory on the head.  Pity
> > that, because I might have been able to do something about it if that
> > was the problem.  Now I'm at a loss to work out why lme4 installation
> > cannot find Matrix, and more strangely, why nothing else gave a
> > similar problem.  I think I've tried every version of lme4 and Matrix
> > that has emerged since R-2.4.0 has been released.
> >
> > The fact that no other Red hat user has a problem indicates the
> > problem is this end; but I'm short of ideas about where to look.
> > Looks like it's back to versions 6 months old -- like proprietary
> > software users have to put up with. :-(
>
> Hmmm, I can provoke this on SUSE too (I was going to ask you to try
> this!):
>
> mkdir foo
> MY_R=~/misc/r-release-branch/BUILD/bin/R
> R_LIBS=foo $MY_R --vanilla << EOF
>  options(repos=c(CRAN="http://cran.dk.r-project.org";))
>  install.packages("lme4", depend=TRUE)
> EOF
>
> This installs Matrix alright, then dies with
>
> Warning in install.packages("lme4", depend = TRUE) :
>  argument 'lib' is missing: using foo
> trying URL 'http://cran.dk.r-project.org/src/contrib/lme4_0.9975-8.tar.gz'
> Content type 'application/x-tar' length 235617 bytes
> opened URL
> ==
> downloaded 230Kb
>
> * Installing *source* package 'lme4' ...
> ** libs
>
> gcc -I/usr/local/src/pd/r-release-branch/BUILD/include
>   -I/usr/local/src/pd/r-release-branch/BUILD/include
>   -I/usr/local/include -fpic -g -O2 -std=gnu99 -c Wishart.c -o Wishart.o
>
> In file included from Wishart.h:4,
>  from Wishart.c:1:
> lme4_utils.h:9:20: Matrix.h: No such file or directory
> In file included from Wishart.h:4,
>  from Wishart.c:1:
> lme4_utils.h:25: error: syntax error before "c"
> lme4_utils.h:25: warning: type defaults to `int' in declaration of `c'
> lme4_utils.h:25: warning: data definition has no type or storage class
> lme4_utils.h:163: error: syntax error before "cholmod_factor"
> lme4_utils.h:177: error: syntax error before "cholmod_factor"
> make: *** [Wishart.o] Error 1
> ERROR: compilation failed for package 'lme4'
>
> And Matrix.h is sitting happily in foo/Matrix/include/ but obviously
> not getting included. Should there have been a copy inside lme4?

Not according to my understanding of how the "LinkingTo:"
specification in the DESCRIPTION file is to work.

The headers for the C functions exported by the Matrix package are
defined Matrix/inst/include so they can be changed in one place if the
API changes.  The "LinkingTo:" specification is designed to put these
headers on the -I path at the time of compilation of the lme4 package
source code.  For some reason it didn't do that in this case.

> What is not obvious to me is how this can work anywhere...
>
> I also tried unpacking the lme4 directory from its tarball, dropping
> all files from the Matrix installed include dir into lme4/src and then
>
> $MY_R CMD INSTALL -l foo lme4
>
> but no go
>
> gcc -shared -L/usr/local/lib64 -o lme4.so Matrix_stubs.o Wishart.o
> glmer.o init.o lme4_utils.o lmer.o local_stubs.o pedigree.o
> -L/usr/local/src/pd/r-release-branch/BUILD/lib -lRlapack
> -L/usr/local/src/pd/r-release-branch/BUILD/lib -lRblas -lg2c -lm
> -lgcc_s
>
> local_stubs.o(.text+0x0): In function `M_numeric_as_chm_dense':
> /home/bs/pd/tmp/lme4/src/Matrix_stubs.c:420: multiple definition of 
> `M_numeric_as_chm_dense'
> Matrix_stubs.o(.text+0x0):/home/bs/pd/tmp/lme4/src/Matrix_stubs.c:420: first 
> defined here

You don't want to compile both Matrix_stubs.c and local_stubs.c.  All
that local_stubs.c does is include Matrix_stubs.c, the reason being
that those wrapper functions are needed locally but it is much cleaner
if they are defined in only one place.  The file Matrix_stubs.c is
never compiled by itself.  It is just there to provide the definitions
of the functions.  Check Luke's documentation and examples for using
R_RegisterCCallable and R_GetCCallable.


> local_stubs.o(.text+0x70): In function `M_dpoMatrix_chol':
> /home/bs/pd/tmp/lme4/src/Matrix_stubs.c:412: multiple definition of 
> `M_dpoMatrix_chol'
> Matrix_stubs.o(.text+0x70):/home/bs/pd/tmp/lme4/src/Matrix_stubs.c:412: first 
> defined here
> ..etc..
>
> --
>O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
>   c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
>  (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
> ~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907
>

Re: [Rd] lme4 install, was Re: [R] Which executable is associated with R CMD INSTALL?

[Douglas Bates]

On 06 Nov 2006 15:41:11 +0100, Peter Dalgaard <[EMAIL PROTECTED]> wrote:
> "Douglas Bates" <[EMAIL PROTECTED]> writes:
>
> > On 31 Oct 2006 12:05:21 +0100, Peter Dalgaard <[EMAIL PROTECTED]> wrote:
> > >
> > > [move to r-devel, put maintainer in loop]
> > >
> > > Patrick Connolly <[EMAIL PROTECTED]> writes:
> > >
> > > > On Mon, 30-Oct-2006 at 04:44PM -0500, Duncan Murdoch wrote:
> > > >
> > > >
> > > > |> Try "R CMD printenv R_HOME" and you'll find which R home directory 
> > > > it is
> > > > |> using.  You can see a lot more with "R CMD printenv" or various 
> > > > options
> > > > |> to "R CMD config".
> > > >
> > > > Thanks for that information.  It knocks my theory on the head.  Pity
> > > > that, because I might have been able to do something about it if that
> > > > was the problem.  Now I'm at a loss to work out why lme4 installation
> > > > cannot find Matrix, and more strangely, why nothing else gave a
> > > > similar problem.  I think I've tried every version of lme4 and Matrix
> > > > that has emerged since R-2.4.0 has been released.
> > > >
> > > > The fact that no other Red hat user has a problem indicates the
> > > > problem is this end; but I'm short of ideas about where to look.
> > > > Looks like it's back to versions 6 months old -- like proprietary
> > > > software users have to put up with. :-(
> > >
> > > Hmmm, I can provoke this on SUSE too (I was going to ask you to try
> > > this!):
> > >
> > > mkdir foo
> > > MY_R=~/misc/r-release-branch/BUILD/bin/R
> > > R_LIBS=foo $MY_R --vanilla << EOF
> > >  options(repos=c(CRAN="http://cran.dk.r-project.org";))
> > >  install.packages("lme4", depend=TRUE)
> > > EOF
> > >
> > > This installs Matrix alright, then dies with
> 
> > > And Matrix.h is sitting happily in foo/Matrix/include/ but obviously
> > > not getting included. Should there have been a copy inside lme4?
> >
> > Not according to my understanding of how the "LinkingTo:"
> > specification in the DESCRIPTION file is to work.
>
> Kurt pointed out the real reason: Using a relative path for R_LIBS was
> asking for trouble.
>
> So it was my own fault, I just got blinded by seeing symptoms so
> close to those reported by Patrick.

Thanks for the clarification.

However, I think it exposes a problem with R CMD INSTALL.  If there is
a LinkingTo: directive in the DESCRIPTION file and the package's
include directory cannot be found I think that R CMD INSTALL should
abort.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Matrix*vector: coercing sparse to dense matrix for arithmetic

[Douglas Bates]

On 11/29/06, Tamas K Papp <[EMAIL PROTECTED]> wrote:

> I have a sparse Matrix (kronecker product of spline design matrices),
> and I need to multiply each row by a number to get another matrix.  If
> the matrix is A and the numbers are stored in a vector k, with plain
> vanilla matrices I would do

> A*k

> But when using the Matrix package (class of A is "dgCMatrix"), I get
> the warning "coercing sparse to dense matrix for arithmetic".  The
> error message is perfectly reasonable, I am looking for a way to do it
> right (keeping operations in the realm of sparse matrices).  Any help
> would be appreciated.

I would suggest premultiplying by a diagonal matrix.  However, when I
look at the method dispatch after an example it doesn't appear to be
making a good choice.  What I think should work well (and Martin or I
will make the necessary changes to ensure that it does) is

library(Matrix)
data(KNex); mm <- KNex$mm
system.time(Diagonal(x = rnorm(nrow(mm))) %*% mm)

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Carriage returns and Sweave output

[Douglas Bates]

On 3/20/07, Ernest Turro <[EMAIL PROTECTED]> wrote:
>
> On 20 Mar 2007, at 07:53, Martin Maechler wrote:
>
> >> "Wolfi" == Wolfgang Huber <[EMAIL PROTECTED]>
> >> on Mon, 19 Mar 2007 15:38:00 + writes:
> >
> >>> the problem with results=hide is that it suppresses everything. I
> >>> just
> >>> need Sweave to suppress strings ending in '\r'...
> >
> >
> > Wolfi> Dear Ernest,
> >
> > Wolfi> IMHO it is good practice to make the printing of these
> > progress reports
> > Wolfi> ("sweep 4 of 1024\r") optional and produce them only if
> > the user calls
> > Wolfi> your function with, say, "verbose=TRUE",
> >
> > I strongly agree.
> >
> > Wolfi> and furthermore set the default value of the
> > Wolfi> 'verbose' argument to "verbose=interactive()".
> >
> > or -- typically my choice -- to  'verbose = getOption("verbose")
> >
> > Martin Maechler, ETH Zurich
> >
> > Wolfi> Best wishes
> > Wolfi> Wolfgang
> >
> >  []
> >
>
> I agree that users should be free to choose the level of verbosity.
> Here, I want to show the verbose output and print it onto the tex
> file using Sweave to give users a good idea of what happens. What I
> don't want is countless lines being printed because \r is being
> interpreted as \n ...

In cases like this capture.output() is your friend.  Write one code
chunk with results=hide,echo=FALSE that uses capture.output to trap
the desired output as character strings then use string manipulation
functions to do the desired replacement.  A second code chunk with
eval=FALSE shows the code that apparently produces the output and a
third code chunk with echo=FALSE uses cat on the manipulated character
strings with quote=FALSE to show what apparently was produced.

>  Ernest Turro wrote:
> > Dear all,
> > I have a code chunk in my Rnw file that, when executed, outputs
> > carriage return characters ('\r') to inform on the progress (e.g.
> > "sweep 4 of 1024\r"). But Sweave interprets this as a newline
> > character, and therefore I get countless pages of output in my
> > vignette where I only really want one line. Any ideas?
> > Thanks
> > E
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [R] Use of 'defineVar' and 'install' in .Call

[Douglas Bates]

I did read your second message about the problem symptoms disappearing
but I thought that I might make a couple of suggestions about your
code anyway.

There are a number of helper functions declared in Rinternals.h such
as ScalarReal, which is equivalent to your mkans.  (Also
ScalarInteger, ScalarLogical, ...) The functions to convert scalars in
the other direction, i.e. taking an SEXP and returning an int or a
double or a ... , are called asInteger, asReal, ...

If you are writing a package you can define an initialization function
called R_init_ to perform initialization for the package.  If
you are going to use a symbol like x many times then you can save the
result of install("x") as a global, say, myPkg_xSymbol during
initialization and use the global variable instead of calling install
many times.  The overhead for calling install is small but if you can
avoid it I don't see a reason not to.

Finally, why do you want to use identifiers like u1, u2, ... when you
could pass a vector named "u" and use that.  In other words, wouldn't
it be easier to do the extraction of the components in the R code for
the function rather than generating a whole bunch of different names?

I suggest that this thread be moved to the r-devel list, which I have cc:d.

On 3/27/07, Daniel Berg <[EMAIL PROTECTED]> wrote:
> Dear all,
>
> [system and version information below]
>
> I am trying to modify a C function for finding the root of an
> expression. The function is to be called from R as .Call with input
> parameters:
>
> f: expression for which we will find the root
> guesses: interval for the solution
> stol: tolerance
> rho: environment
>
> The original functions I use are:
>
> SEXP mkans(double x) {
>   SEXP ans;
>   PROTECT(ans = allocVector(REALSXP, 1));
>   REAL(ans)[0] = x;
>   UNPROTECT(1);
>   return ans;
> }
> double feval(double x, SEXP f, SEXP rho) {
>   defineVar(install("x"), mkans(x), rho);
>   return(REAL(eval(f, rho))[0]);
> }
> SEXP zero(SEXP f, SEXP guesses, SEXP stol, SEXP rho) {
>   double x0 = REAL(guesses)[0], x1 = REAL(guesses)[1], tol = REAL(stol)[0];
>   double f0, f1, fc, xc;
>   if(tol <= 0.0) error("non-positive tol value");
>   f0 = feval(x0, f, rho); f1 = feval(x1, f, rho);
>   if(f0 == 0.0) return mkans(x0);
>   if(f1 == 0.0) return mkans(x1);
>   if(f0*f1 > 0.0) error("x[0] and x[1] have the same sign");
>   for(;;) {
> xc = 0.5*(x0+x1);
> if(fabs(x0-x1) < tol) return mkans(xc);
> fc = feval(xc, f, rho);
> if(fc == 0) return mkans(xc);
> if(f0*fc > 0.0) {
>   x0 = xc; f0 = fc;
> }
> else {
>   x1 = xc; f1 = fc;
> }
>   }
> }
>
>
> This works great. However, I wish to make it more general, by
> modifying 'feval'. Given that my problem involves a data set 'u', with
> dimension (i x j), I need to assign values to 'u1', 'u2', ..., 'ui'
> via defineVar(install(...)). I tried the following:
>
> double feval(double x, double *u, int d, double v, SEXP f, SEXP rho) {
>   int i;
>   char *str1="u", str2[1001], *str3;
>   defineVar(install("x"), mkans(x), rho);
>   defineVar(install("y"), mkans(v), rho);
>   for(i=0;i sprintf(str2,"%d",i+1);
> str3 = (char *)malloc((strlen(str1)+strlen(str2)+1)*sizeof(char));
> strcpy(str3,str1);
> strcat(str3,str2);
> defineVar(install(str3), mkans(u[i]), rho);
>   }
>   free(str3);
>   return(REAL(eval(f,rho))[0]);
> }
>
> My R-package still compiles without errors but R crashes due to the
> defineVar command.
>
> Any suggestions to how I can do the defineVar bit?
>
> Thanks in advance.
>
> Reagards,
> Daniel Berg
>
> 
> > R.Version()
> $platform
> [1] "i486-pc-linux-gnu"
> $arch
> [1] "i486"
> $os
> [1] "linux-gnu"
> $system
> [1] "i486, linux-gnu"
> $status
> [1] ""
> $major
> [1] "2"
> $minor
> [1] "3.1"
> $year
> [1] "2006"
> $month
> [1] "06"
> $day
> [1] "01"
> $`svn rev`
> [1] "38247"
> $language
> [1] "R"
> $version.string
> [1] "Version 2.3.1 (2006-06-01)"
>
> __
> R-help@stat.math.ethz.ch 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-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Developer work cycle

[Douglas Bates]

On 3/26/07, "José Luis Aznarte M." <[EMAIL PROTECTED]> wrote:
> Hi! I've been browsing through the last months' archive and I can't
> find an answer to my question, so here it is (let's hope it's not too
> obvious):
> I'm working on extensions of an R library, and I would be very
> surprised if everyone developing R packages is doing the following, as I do:
>
> 1.- Write down a modification of an R file
> 2.- Exit the current R session
> 3.- Install the package as root (sudo R CMD INSTALL...)
> 4.- Launch a new R session
> 5.- Test the change, if it does not work, go back to 1 or debug.
> 6.- Finish.
>
> Is this the proper (but quite awkward) way to proceed or there is an
> alternative to skip steps 2 to 4? I'm using emacs with ESS under linux.

John Chambers has provided an alternative approach of using

trace(fname, edit = TRUE)

where fname is the name of your function.  (Make sure that the server
for emacsclient has been started in your emacs session with M-x
server-start.)  This opens an emacs buffer containing the source for
the function which you can then edit.  After writing the file and
closing the client (C-x #) your ESS session has the new definition
installed in the package's namespace.

This will work even for objects hidden in the namespace.  The argument
"signature" allows you to edit S4 methods on the fly like this.  In my
experience you cannot edit registered S3 methods like this but it may
be that I am just not using trace correctly.

Of course you must copy the modified version of the source code to
your package sources when you are done.

As others have indicated, it is a good practice to install development
versions of packages in a private library so you do not need to use
sudo or worry about messing up system-wide directories.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] missing IntegerFromString()

[Douglas Bates]

On 6/6/07, Aniko Szabo <[EMAIL PROTECTED]> wrote:
> Thanks to everybody who responded to my question.
> asInteger(coerceVector(x,INTSXP)) indeed does what I need. I guess there
> is a lot I don't understand about type coercion, as I would not have
> expected it to work.

It is better to use

asInteger(x)

which will do the coercion if necessary.  When you do the coercion
yourself you should PROTECT the result then UNPROTECT it.  Calling
asInteger directly avoids this overhead without the risk of losing
data in a garbage collection.  asInteger can accomplish this because
only the first element of the SEXP x is converted to an integer.

>
> Aniko
>
> -Original Message-
> From: Seth Falcon [mailto:[EMAIL PROTECTED]
> Sent: Tuesday, June 05, 2007 11:24 PM
> To: Aniko Szabo
> Cc: r-devel@r-project.org
> Subject: Re: [Rd] missing IntegerFromString()
>
> Hi Aniko,
>
> "Aniko Szabo" <[EMAIL PROTECTED]> writes:
>
> > I have created a DLL not so long ago using C code. My code used the
> > IntegerFromString() function that used to be exported in the
> > Rinternals.h header file (and thus easily accessible). Recently I
> > upgraded to R 2.5.0 and my DLL stopped working. I see that the
> > IntegerFromString() function is not exported in any of the header
> files
> > in the RHOME\include directory. Is it possible for me to use it
> without
> > installing all R source files? I can see that the function is in
> > coerce.c, however it #includes other stuff that I don't have and I am
> > afraid to mess things about by doing things I don't understand. Or
> > perhaps there is another function that is intended to be used
> > instead?
>
> I think you want asInteger (which calls IntegerFromString).  This is
> in RHOME/include/Rinternals.h
>
> Best Wishes,
>
> + seth
>
> PS: Nice to see you again :-)
>
>
> --
> Seth Falcon | Computational Biology | Fred Hutchinson Cancer Research
> Center
> http://bioconductor.org
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] GPL 3 has been released

[Douglas Bates]

Version 3 of the GNU Public License (GPL) has been released.  A ZDNet
article about it can be found at

http://news.zdnet.com/2100-3513_22-6194139.html?part=rss&tag=feed&subj=zdnn

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Forthcoming change in the API of the Matrix package

[Douglas Bates]

Martin and I will soon release a new version of the Matrix package
with a modified API.  This will affect the authors of any packages
that use calls to the C function R_GetCCallable to directly access C
functions in the DLL or shared object object in the libs directory of
the Matrix package.  (If you didn't understand that last sentence,
relax - it means that you can ignore this message.)

We strongly suspect that I am the only such author (this mechanism is
used in the lme4 package) and, because I was the one who made the API
change, I do indeed know about it.  However, if others do use this
mechanism for, say, accessing functions in the CHOLMOD sparse matrix C
library, you should be aware of this.

The current version of the Matrix package is 0.99875-3.  This version
exports the C functions according to the old API.  The next version
will be 0.999375-0 using the new API.  I will soon upload version
0.99875-3 of the lme4 package that depends on

Matrix(<= 0.99875-3)

Version 0.999375-0 of the lme4 package will depend on

Matrix(>= 0.999375-0)

The changes in the API are in the functions as_cholmod_sparse,
as_cholmod_dense and as_cholmod_factor.  After the change the first
argument will be a pointer to a struct of the appropriate return type
(i.e. the first argument in as_cholmod_sparse is a cholmod_sparse *
and the second argument is an SEXP).  This allows the calling function
to handle both the allocation and the freeing of the storage for the
struct.

Also the new API provides several macros and typedefs for such
pointers to structs.

The development version of the Matrix package is available at

https://svn.R-project.org/R-packages/branches/Matrix-APIchange/

The corresponding  version of the lme4 package is at

https://svn.R-project.org/R-packages/branches/gappy-lmer/

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] [OT] help in setting up a doxygen configuration file

[Douglas Bates]

I would appreciate some pointers on how to set up a doxygen
configuration file for C source code.  In particular I would like to
be able to generate a call graph.  I tend to write a lot of short
utility functions and, by the time the final design reveals itself, it
is quite possible that some of these utilities are called in only one
place.  That's not a bad thing to have happen but I would like to know
about it when it does occur.

Doxygen seems to emphasize C++ classes and I can't manage to get it to
do much with my C functions.  The package sources are available at

https://svn.R-project.org/R-packages/branches/gappy-lmer/

The doxygen configuration file is gappy-lmer/inst/doc/Doxyfile

I have written all the Javadoc-style comments in the source files such
as gappy-lmer/src/lmer.c but I can't seem to get doxygen to notice
them.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Changes in the MEMSS package (data sets from Pinheiro and Bates, 2000, without the groupedData classes)

[Douglas Bates]

Some time ago Deepayan and I created a package called MEMSS for the
data sets from the nlme package  as data frames but not groupedData
objects.  Because of advances that Deepayan has made in lattice
graphics many of the specialized plotting methods for the groupedData
objects are no longer needed.  It is easier and less confusing to
store these data sets as data frames rather than as groupedData
objects.

I propose making two further changes in these objects.  Those objects
typically had factors converted to ordered factors in an artificial
way and I plan to revert them to factors.  Also, I will relabel any
factors with (26 or fewer) numeric levels with letters and to remove
the "ginfo" attribute.

Please tell me if this will cause you hardship because you have used
these data sets in your work and depend on the current formulation.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Changes in the MEMSS package (data sets from Pinheiro and Bates, 2000, without the groupedData classes)

[Douglas Bates]

On 7/24/07, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> Does this include datasets such as CO2 and ChickWeight
> which are in the datasets package?

> Could you post a list of the specific datasets you are referring to
> so there is no confusion what this is about.

I am only referring to the datasets in the MEMSS package.  There will
be CO2, ChickWeight, Theoph, etc. datasets in the MEMSS package
without the groupedData classes.  However, the versions in the
datasets package will stay as they are unless and until R-core decides
to change them.

Having datasets in the datasets package and in the MEMSS package with
the same name will produce a warning when you attach the MEMSS
package.  However, I generally use datasets with a call like

data(Rail, package = "MEMSS")

which doesn't cause a warning (it is assumed that if you are this
explicit about the data set then you probably want exactly that
version and haven't accidently overridden the other).

The bottom line is that unless you have been in the habit of attaching
the MEMSS package or getting data sets from it as above then you won't
notice any difference.

> On 7/24/07, Douglas Bates <[EMAIL PROTECTED]> wrote:
> > Some time ago Deepayan and I created a package called MEMSS for the
> > data sets from the nlme package  as data frames but not groupedData
> > objects.  Because of advances that Deepayan has made in lattice
> > graphics many of the specialized plotting methods for the groupedData
> > objects are no longer needed.  It is easier and less confusing to
> > store these data sets as data frames rather than as groupedData
> > objects.
> >
> > I propose making two further changes in these objects.  Those objects
> > typically had factors converted to ordered factors in an artificial
> > way and I plan to revert them to factors.  Also, I will relabel any
> > factors with (26 or fewer) numeric levels with letters and to remove
> > the "ginfo" attribute.
> >
> > Please tell me if this will cause you hardship because you have used
> > these data sets in your work and depend on the current formulation.
> >
> > __
> > R-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Compiling R for the Sony Playstation 3?

[Douglas Bates]

Has anyone installed Linux on a Sony Playstation 3 and compiled R for it?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Compiling R for the Sony Playstation 3?

[Douglas Bates]

On 8/1/07, Marc Schwartz <[EMAIL PROTECTED]> wrote:
> On Wed, 2007-08-01 at 11:13 -0500, Douglas Bates wrote:
> > Has anyone installed Linux on a Sony Playstation 3 and compiled R for it?

> Doug,
>
> I don't have any personal experience with both Linux and R on the PS3,
> but do know folks who have run Linux successfully on that platform.
>
> Here are some links that you might find helpful to at least cover the
> first part:
>
> http://www.playstation.com/ps3-openplatform/manual.html
>
> http://en.wikipedia.org/wiki/Linux_for_PlayStation_3
>
> http://www.engadget.com/2006/11/19/fedora-linux-up-and-running-on-playstation-3-with-video/
>
>
> There are additional links, as usual, on the Wikipedia page.

I did successfully install Ubuntu 7.04 ("feisty") on a Playstation 3.
Once you have Ubuntu installed you can use the standard package
management tools for Debian/Ubuntu to install the r-base-core,
r-recommended and r-base-dev packages.  Those packages are for R-2.4.1
but I had no trouble installing and checking the most recent versions
of R from the SVN site.

Essentially the Playstation 3 becomes a typical Ubuntu system with all
the tools from Ubuntu.  At present X11 uses the framebuffer device so
the graphics is poor by X11 standards and very poor by gamer's
standards.

Overall it is quite amazing given that the PS3 with a 60GB hard drive
is now being sold in the U.S. for $500 and that includes a Blu-ray DVD
drive and full 1080p HDMI interface.  At $500 it is one of the
cheapest Blu-ray DVD players, which is what I bought it for.  It's
great to be able to run R on my DVD player.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Compiling R for the Sony Playstation 3?

[Douglas Bates]

P.S. A site with instructions for installing Ubuntu on the PS3 is
psubuntu.com or google "PS3 Ubuntu" to get the community pages
documentation.

It is a good idea to have an Ethernet connection so that once you
install from the CD you can upgrade the packages.  The kernel on the
installation CD image has a known bug that causes it to flood the
system log with messages about querying the memory card readers.  You
want to replace that kernel with a later version.

The processor is apparently a 3.2 GHz dual-core PowerPC-64.  It
doesn't feel as fast as I expected but that could be because of memory
bandwidth (256 MB of memory plus 256 MB of graphics memory that isn't
used by Linux at present) or disk access.  However, it still runs R
faster than any other DVD player I have ever had.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Compiling R for the Sony Playstation 3?

[Douglas Bates]

On 8/3/07, Paul Gilbert <[EMAIL PROTECTED]> wrote:
> Doug

> Does it still play DVDs?

> (My wife keeps telling me that academics have to do all their serious
> work in the summer.)

The way that you install Linux on a PS3 it becomes a dual-boot
machine.  You have to reboot to the original operating system (called
XBM) to play DVD's at present.

Sony made provision for a second operating system to be installed on
the hard drive.  In the 'out-of-the-box' state there is a selection
under the 'System settings' that allows you to reformat the hard drive
with two partitions.  Unfortunately they only offer two possibilities
- either 10 Gb for XBM, 50 Gb for other or 50 Gb for XBM and 10 Gb for
other.  Reformatting does destroy any saved data but there is
provision for saving and restoring such data from a SD memory card or
a USB device.  Reformatting does not destroy the XBM operating system
(I believe it is stored in flash memory).

After reformatting the hard drive you insert the installation CD,
select "Install other OS"  from the System Settings and sit back and
watch the installation.

Once you install Linux the system boots to a boot loader called kboot.
 If you type boot-game-os at the kboot prompt you end up in XBM and
can play DVD's to your heart's content (I understand that you can also
play video games with it :-)

By the way, you can buy a remote control for the PS3 as a DVD player
but you don't need to have it.  You can use the game controller as a
DVD remote and I actually prefer it to having yet another remote
control with 85 buttons of which I use 2 and can never decide which 2
those are without turning on a light and getting reading glasses out.
The game remote is intended to be operated without looking at it and
all the feedback is on the screen.  It doesn't take long to become
proficient with it.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Sweave/ESS-like tools for HTML

[Douglas Bates]

My university provides me with a powerful course management system for
the courses that I teach.  Among other things I can create a wiki for
the course, which is very convenient for cross-linking different bits
of the course.

Naturally I use R extensively in my teaching and I want to incorporate
R code, output and graphics in such a wiki.  If I were producing LaTeX
sources instead of HTML sources I create .Rnw files for Sweave and I
would edit them using ESS in emacs.

What options do I have for producing HTML with embedded R content and
what is a good, preferably emacs-based, way of editing the source
code?

One basic problem is trying to present mathematical expressions in HTML (see
http://www.cs.tut.fi/~jkorpela/math/) but, aside from that, there are
questions of presenting input R expressions and the corresponding
output and of incorporating graphics files produced by R.  I could try
to use latex2html or texi2html but the output from latex2html at least
would be quite inconvenient to use because it generates so many linked
files.  Once they are uploaded it would be horrible trying to get all
the links straightened out.

In a sense there are already tools for this type of output from .Rd
files.  Would it be best to use those tools or to use texinfo tools or
...?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Sweave/ESS-like tools for HTML

[Douglas Bates]

On 10/18/07, Ben Bolker <[EMAIL PROTECTED]> wrote:

> Tom Short-2 wrote:

> > See this link for more on creating/converting to HTML:

> > http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/SweaveConvert

> > For using ESS with mixed HTML/R files, see this:

> > https://stat.ethz.ch/pipermail/ess-help/2006-December/003826.html

>  A couple of other bits of information:
> (1) I started writing a LaTeX-to-Rwiki translator,
> anyone who wants to try it should contact me.
> (It sounds, though, as though the wiki that Doug
> Bates is going to use accepts HTML rather than
> the Rwiki format?

> (2) I have trouble getting through to Ian Hutchinson's
> web site to get tth.  If you do get there, check out
> ttm (TeX-to-MathML) as well.

I found the links that Tom sent to be very helpful and did install
both the tth and  tex4ht Debian packages.

Both approaches are interesting but there are some inherent
limitations to HTML and the course management system that are
difficult to overcome.  I haven't found a way to include an expression
like $\bar{x}$ without resorting to images and the course management
system goes to great lengths to hide the file hierarchy so including
an image in part of the wiki is difficult.

Perhaps it is better if I create stand-alone PDF documents but then
linking becomes difficult.


>   The SGML approach does seem like the wave of
> the future, but in the meanwhile TeX/tth works well
> for me.
>
>   Ben Bolker
> --
> View this message in context: 
> http://www.nabble.com/Sweave-ESS-like-tools-for-HTML-tf4635106.html#a13272187
> Sent from the R devel mailing list archive at Nabble.com.
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] nlminb - Port library codes

[Douglas Bates]

You could always look in the SVN logs but I'll save you the trouble.
I did the initial installation of the nlminb function.

On Fri, May 1, 2009 at 12:52 PM, John C Nash  wrote:
> Having looked at documentation and codes, I'm still looking to find out who
> did the installation of the Port libraries. As I work on optimization code
> improvements, there are often choices of settings (tolerances etc.), and I'd
> prefer to learn if there are particular reasons for choices before diving in
> and making any adjustments.
>
> It is also worth giving credit where it is due. Given the complexity of the
> Port code structure, adapting them to R must have been a lot of work.
>
> JN
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] beginner's guide to C++ programming with R packages?

[Douglas Bates]

On Fri, Jun 26, 2009 at 2:43 PM, Paul Johnson wrote:
> Hello, again.
>
> I'm interested to learn how programmers develop & test C/C++ code with
> R packages in Linux.  I've been reading R source and the manual on
> Writing R Extensions but there are just a couple of details I can't
> understand.  I wish I could watch over a developer's shoulder to see
> how people actually do this.
>
> I've tested a bit.  I am able to take package.tar.gz file, open it up,
> fiddle the source code, and then run
>
> R CMD check package-dir
>
> from the directory above "package-dir" ,
>
> R CMD build package-dir
>
> and
>
> R CMD INSTALL

For your purposes it is probably better to avoid building a new tar.gz
file and combine the last two steps as

R CMD INSTALL package-dir

The install process uses the make utility which will check which of
the object files need to be recompiled.  If you only modify one source
file it will be the only file recompiled.
> on the tarball that is produced. Then in R, I can load the package and use it.
>
> That part is "all good", but somewhat tedious.  I don't want to
> entirely recompile and reinstall the whole package just to test one
> function.  I notice that R CMD check creates a new directory called
> "package.Rcheck" and the shared objects and example code of the
> package are in there.  Can I force R to use those *.so files instead
> of the ones in /usr/lib/R ?
>
>
> I also wonder "what is wrong with gprof?   In the Writing R Extensions
> manual, it describes "oprofile" and "sprof" for Linux. I will try
> them, but they are unfamilar to me.  I've used gprof in the past in C
> projects, and it is a pretty painless thing to add a compiler flag
> -pg, run the program, and then review gmon.out.  The Writing R
> Extensions manual does not mention gprof in its section on Linux, but
> it does mention it under Solaris.  There is a somewhat ambiguous
> statement:
>
> 3.4.2 Solaris
>
> On 64-bit (only) Solaris, the standard profiling tool gprof collects
> information from shared libraries compiled with -pg.
>
> Does "(only)" here mean to differentiate Solaris from other Linux/Unix
> systems?  Or does it differentiate 64bit Solaris from other Solaris?
>
> But this draws me back to the basic question.  I don't want to run R
> CMD INSTALL 20 times per hour.  How do developers "actually" test
> their code?
>
> pj
> --
> Paul E. Johnson
> Professor, Political Science
> 1541 Lilac Lane, Room 504
> University of Kansas
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Problem using model.frame with argument subset in own function

[Douglas Bates]

On Sat, Aug 8, 2009 at 1:31 PM, Gavin Simpson wrote:
> Dear List,

> I am writing a formula method for a function in a package I maintain. I
> want the method to return a data.frame that potentially only contains
> some of the variables in 'data', as specified by the formula.

The usual way to call model.frame (the method that Thomas Lumley has
called "the standard, non-standard evaluation) is to match the call to
foo, replace the name of the function being called with
as.name("model.frame") and force an evaluation in the parent frame.
it looks like

mf <- match.call()
if (missing(data)) data <- environment(formula)
## evaluate and install the model frame
m <- match(c("formula", "data", "subset", "weights", "na.action", "offset"),
   names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
fr <- eval(mf, parent.frame())

The point of all of this manipulation is to achieve the kind of result
you need where the subset argument is evaluated in the correct
environmnent.

> The problem I am having is in writing the function and wrapping it
> around model.frame. Consider the following data frame:
>
> dat <- data.frame(A = runif(10), B = runif(10), C = runif(10))
>
> And the wrapper function:
>
> foo <- function(formula, data = NULL, ..., subset = NULL,
>                na.action = na.pass) {
>    mt <- terms(formula, data = data, simplify = TRUE)
>    mf <- model.frame(formula(mt), data = data, subset = subset,
>                      na.action = na.action)
>    ## real function would do more stuff here and pass mf on to
>    ## other functions
>    mf
> }
>
> This is how I envisage the function being called. The real world use
> would have a data.frame with tens or hundreds of components where only a
> few need to be excluded. Hence wanting formulas of the form below to
> work.
>
> foo(~ . - B, data = dat)
>
> The aim is to return only columns A and C in an object returned by
> model.frame. However, when I run the above, I get the following error:
>
>> foo(~ A + B, data = dat)
> Error in xj[i] : invalid subscript type 'closure'
>
> I've tracked this down to the line in model.frame.default
>
>    subset <- eval(substitute(subset), data, env)
>
> After evaluating this line, subset contains:
>
> Browse[1]> subset
> function (x, ...)
> UseMethod("subset")
> 
>
> Not NULL, and hence the error later on when calling the internal
> model.frame code.
>
> So the question is, what am I doing wrong?
>
> If I leave the subset argument out of the definition of foo and rely
> upon the default in model.frame.default, the function works as
> expected.
>
> Perhaps the question should be, how do I modify foo() to allow it to
> have a formal subset argument, passed to model.frame?
>
> Any other suggestions gratefully accepted.
>
> Thanks in advance,
>
> G
> --
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>  Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
>  ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
>  Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
>  Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
>  UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Calling C functions with value parameters

[Douglas Bates]

On Mon, Aug 17, 2009 at 9:23 AM, Jeffrey J. Hallman wrote:
> One hassle I could do without is the necessity of writing C wrapper functions
> like this:
>
> void fameInit(int *status){
>  cfmini(status);
>  return;
> }
>
> when I want to call a library function (cfmini, in this case) that takes an
> int argument.  The .C interface only lets me pass a pointer to an int, rather
> than the int itself.

> Is there any chanch that .C could be enhanced to allow passing arguments by
> value to the compiled code?  It would make some of my stuff much simpler.

I suppose that if you design a new interface and submit patches to the
code in R that does the interfacing, it could be done.  However, I
imagine you would find it easier to continue to write short interface
functions like that. :-)

(By the way, the "return;" is redundant in a C function that returns void.)

You should remember that internally there are no scalar objects in R -
everything is a vector.  If you use the .Call interface, which
initially is more complicated but eventually helps you simplify your
code, arguments are passed and values returned as SEXPs and you can
use convenient utilities like asInteger which make your R code simpler
because they extract a scalar integer value from an SEXP, coercing if
necessary.

So, if you have a C function that takes two integer arguments and
returns an integer you can write the C interface function as

SEXP R_callable_foo(SEXP a, SEXP b)
{
return ScalarInteger(foo(asInteger(a), asInteger(b)));
}


It is not surprising that you need to write interface code - it is
remarkable that it works at all.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] clearNames and unname

[Douglas Bates]

On Thu, Sep 3, 2009 at 10:02 AM, Martin
Maechler wrote:
>> "MM" == Martin Maechler 
>>     on Thu, 3 Sep 2009 16:14:24 +0200 writes:
>
>> "HW" == Hadley Wickham 
>>     on Thu, 3 Sep 2009 08:16:27 -0500 writes:
>
>    HW> Just noticed these two functions (clearNames is stats and unname in
>    HW> base) that do the same thing.
>
>    MM> clearNames looks like an accident.  unname() has existed for
>    MM> much longer {and as someone else said},
>    MM> has been more general, too.
>
> Hmm, the first part in the above may be wrong:
> clearNames() has been part of the 'nls' package (before that
> was integrated into 'stats'),
> and the code of 'nls' may well have been almost as old as the
> unname() function :
> I see both in  R 0.90.1  (December 15, 1999).
>
> I see only about 5 usages of clearNames() in CRAN packages, and
> none in the 300+ Bioconductor "Rpacks"  packages.
>
> Consequently, I'm proposing to deprecate  clearNames() {for
> unname()}.

Fine with me.  I had forgotten that we even wrote clearNames() for the
nls package.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Non-GPL packages for R

[Douglas Bates]

On Fri, Sep 11, 2009 at 1:48 PM, rudjer  wrote:
>
> Comrades,
>
> When talk turns to the purity of the revolution, and purge of packages then
> the guillotine can't be far behind.  We all remember Lenin berating the
> "renegade Kautsky" for his "pragmatism," and we know where that led...
>
> So let me put in a good word for pragmatism, and incidentally for saving one
> of my
> own packages, SparseM, and perhaps eventually my neck.  Last week Kurt asked
> me to look into a SparseM licensing quirk based on an inquiry from the
> Fedora
> folks.  SparseM is GPL except for one routine cholesky.f written at Oakridge
> Lab by E. Ng and B. Peyton.  Our version of the code was redistributed in
> the
> package PCx which was copywrited by the U. of Chicago, who specified that
> commercial users should contact someone at Argonne National Lab.  Since the
> beginning we have retained this language in the License file of SparseM,
> even
> though the code in question was not actually developed as a part of PCx.
>
> I contacted one of the original PCx developers who responded as follows:
>
>        The routine you mention was distributed with PCx but not part
>        of it as you see from the legalese and not covered by the PCx
>        copyright.  I tried to interest the authors of that code
>        in legal issues in around 1997 but could not get them
>        motivated (frankly I also can't get too interested).
>
> To which I heartily concurred.  If someone who is worried about getting sued
> would like to dig into this can of worms, then fine.  But life is too short
> for the rest of us.  This is quite a murky business, we shouldn't create
> incentives to make it murkier by covering up relevant language on licensing.
> But surely we can also all agree that CRAN has been a fantastic success, and
> adding new constraints on its operation is ill-advised.

It is unfortunately common in the numerical analysis community,
especially those still using Fortran, to have a rather vague approach
to licensing. ("I'll send a copy of my code to anyone who asks for it
but put in some language that if someone is going to get fantastically
rich from it then they owe me money too")  In the Open Source
community licensing is very important - it is what makes Open Source
software, including CRAN, possible. Most non-lawyers don't find the
study and discussion of licenses to be terribly fascinating but they
are the foundation of Open Source software. If the authors of Fortran
subroutines feel that it is too much of a bother to pay attention to
licenses (or to learn post-1950's programming languages) then
evolution will run its course and they will be left behind.  It's
annoying in that so little software from the numerical analysis
community is covered by suitable licenses but that will change.

Tim Davis's C and C++-based sparse matrix code that is incorporated in
the Matrix package is licensed under the GPL or LGPL.  Why mess around
with antiquated software and vague or non-existent licenses when there
are better alternatives?  It is painful to need to recode old Fortran
routines in modern programming languages and under real licenses but
it is the only way we will ever bring numerical analysis into the
post-Beatles era.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R + C + Lapack toy regression example

[Douglas Bates]

On Wed, Sep 23, 2009 at 2:39 AM, Vinh Nguyen  wrote:
> dear list,
>
> since matrix manipulations is often of interest in statistical
> computations, i'd like to get a working example of using Lapack for
> regression.  However, i run into an error.
>
> My matrix-lapack-example.c file:
> #include 
>
> void reg(const char* trans, const int* m, const int* n,
>         const int* nrhs, double* a, const int* lda,
>         double* b, const int* ldb,
>         double* work, const int* lwork, int* info)
> {
>  F77_CALL(dgels)(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info);
> }
>
> My matrix-lapack-example.R file:
> dyn.load( "matrix-lapack-example.so" )
>
> regression <- function( X, Y ){
>  m <- nrow( X )
>  n <- ncol( X )
>  stopifnot( m >= n , is.vector( Y ), n != length( Y ) )
>  mn <- min( m, n )
>  lwork <- mn + max( mn, 1 )
>  ### to be used with dgels
>  ### http://www.netlib.org/lapack/double/dgels.f
>  rslt <- .C( "reg"
>             , trans=as.character( "N" )

In your C code you are treating that object as a C character string
(char*) but that is not what is being passed.  Look at section 5.2 of
"Writing R Extensions" - a character variable in R is passed by .C as
a char**, not a char*.

I find it much easier to use the .Call interface instead of the .C
interface for applications like this.  You do need to learn some of
the mysteries of the functions declared in Rinternals.h but, once you
do, it is much easier to pass a matrix from R with and extract all the
fussy details within the C code.  Several of the C source code files
in the Matrix package are devoted to exactly the kind of operation you
want to carry out.  Look at the function lsq_dense_QR in
Matrix/src/dense.c, for example. (Although now that I look at it
myself I see some questionable programming practices - I should have
PROTECTed the result of coerceVector but it happens that it would not
have needed protection.  Rather than coercing I should just check
isInteger on the dim attribute.)
>             , m=as.integer( m ), n=as.integer( n )
>             , nrhs=as.integer( 1 ), a=as.double( X )
>             , lda=as.integer( m ), b=as.double( Y )
>             , ldb=as.integer( m ) , work=double( lwork )
>             , lwork=as.integer( lwork )
>             , info=integer( 1 ) )
>  ##F77_NAME(dgels)(const char* trans, const int* m, const int* n,
>  ##                const int* nrhs, double* a, const int* lda,
>  ##                double* b, const int* ldb,
>  ##                double* work, const int* lwork, int* info);
>  return( rslt )
> }
>
> n <- 100
> x <- rnorm( 100 )
> y <- 2.5 + 3*x + rnorm( n )
> X <- cbind( 1, x )
>
> temp <- regression( X=X, Y=y )
>
>
> dgels does linear least squares.  the C code compile fines with a
> warning (ld: warning, duplicate dylib
> /usr/local/lib/libgcc_s.1.dylib).  in R, i get the following when i
> run regression():
>> temp <- regression( X=X, Y=y )
> Parameter 1 to routine DGELS  was incorrect
> Mac OS BLAS parameter error in DGELS , parameter #0, (unavailable), is 0
>
> Process R exited abnormally with code 1 at Wed Sep 23 00:21:59 2009
>
> am i doing something wrong?  is my as.character() used wrong here?
>
> i know R uses fortran code for lm.  but suppose i wanted to do
> regression in C/C++.  is this lapack method the recommended / "best
> practice" route?  i just want to know whats the recommended method for
> doing matrix manipulations in C.  thanks!
>
> vinh
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R + C + Lapack toy regression example

[Douglas Bates]

On Thu, Sep 24, 2009 at 2:09 PM, Vinh Nguyen  wrote:
> On Thu, Sep 24, 2009 at 11:49 AM, Simon Urbanek
>  wrote:
>> As Doug pointed out you don't want to be using .C(). As for matrix
>> manipulations - they are usually done directly on the objects which are
>> vectors stored in column-major order.
>
> i meant .Call().   also, sorry for the poor word choice, i meant
> matrix computations, not matrix manipulations.  i guess i just want
> some recommendations for future references on what C library to use if
> i were doing something like re-implement linear regression
> (hypothetically speaking).  i assumed it would be lapack.  is this the
> recommended approach?

Yes, Lapack would be the recommended code for numerical linear
algebra.  If happens to be Fortran code but it is callable from C/C++.
 One big advantage of Lapack (besides its having been written by some
of the top people in the field and being very well tested) is that it
uses the BLAS (Basic Linear Algebra Subroutines) as levels 1, 2 and 3.

However, implementing linear models software is more subtle than just
the numerical part.  In fact, least squares calculations are (I think)
the only part of R where Linpack is used instead of  Lapack, by
default.  Most of the time the symbolic analysis of a linear model
formula done in R produces a model matrix with full column rank but
under some circumstances it doesn't.  (The easiest such case to
explain is a two-way layout with missing cells.)  The way these cases
are handled is to use a special pivoting algorithm in the
decomposition.  Neither Linpack nor Lapack provided a pivoting scheme
so early in the R project Ross modified the Linpack version of the QR
decomposition to use such a scheme.  It is still the one in use.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] unexpected behaviour of isSymmetric() (PR#14000)

[Douglas Bates]

On Mon, Oct 12, 2009 at 6:41 AM, Romain Francois
 wrote:
> On 10/12/2009 02:05 AM, m...@stat.ubc.ca wrote:
>>
>> Full_Name: Mike Danilov
>> Version: 2.9.0
>> OS: Fedora Core 9
>> Submission from: (NULL) (142.103.121.198)
>>
>>
>> When checking for the symmetry of a matrix, function isSymmetric.matrix()
>> gets
>> confused by the discrepancy of colnames/rownames if its argument. See the
>> code
>> snippet below. Perhaps it's a problem of the matrix product which copies
>> colnames of the first argument but not the rownames of the second to its
>> value.
>> Not sure which one should be fixed but the way it is now it seems
>> illogical that
>> X'X is deemed to be non-symmetric.
>>
>> x<- c(1,2,3)
>> names(x)<- c("v1","v2","v3")
>> isSymmetric(x%*%t(x)) ## returns FALSE instead of TRUE
>
> It seems to be concerned with the names
>
>> y <- x %*% t(x)
>> y
>     v1 v2 v3
> [1,]  1  2  3
> [2,]  2  4  6
> [3,]  3  6  9
>
>> isSymmetric( y )
> [1] FALSE
>
> # dropping the dimnames
>> isSymmetric( structure( y, dimnames = NULL ) )
> [1] TRUE
>
> # pass the ... along this path : isSymmetric > all.equal > all.equal.numeric
>> isSymmetric( y, check.attributes = F )
> [1] TRUE
>
> # set the dimnames equal
>> dimnames( y ) <- rep( list( names(x) ), 2 )
>> isSymmetric( y )
> [1] TRUE
>
> Not sure this is expected behaviour
>
> Romain

I think the problem is more with the propagation of the column names
in the construction x %*% t(x).  If you use the tcrossprod function to
create x %*% t(x) more carefully then the results are sensible


> x<- c(1,2,3); names(x)<- c("v1","v2","v3")
> isSymmetric(tcrossprod(x))
[1] TRUE
> tcrossprod(x)
 [,1] [,2] [,3]
[1,]123
[2,]246
[3,]369

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Confusion regarding allocating Matrices.

[Douglas Bates]

On Fri, Oct 23, 2009 at 8:39 AM, Abhijit Bera  wrote:
> Hi
>
> I'm having slight confusion.

Indeed.

> I plan to grow/realloc a matrix depending on the data available in a C
> program.

> Here is what I'm tried to do:

> Data=allocMatrix(REALSXP,3,4);
> SEXP Data;

Those lines should be in the other order, shouldn't they?

Also, you need to PROTECT Data or bad things will happen.

> REAL(Data)[8]=0.001123;
> REAL(Data)[20]=0.001125;
> printf("%f %f\n\n\n\n",REAL(Data)[8],REAL(Data)[20]);

> Here is my confusion:

> Do I always require to allocate the exact number of data elements in a R
> Matrix?

Yes.

> In the above code segment I have clearly exceeded the number of
> elements that have been allocated but my program doesn't crash.

Remember that when programming in C you have a lot of rope with which
to hang yourself.   You have corrupted a memory location beyond that
allocated to the array but nothing bad has happened  - yet.

> I don't find any specific R functions for reallocation incase my data set
> grows. How do I reallocate?

You allocate a new matrix, copy the contents of the current matrix to
the new matrix, then release the old one.  It gets tricky in that you
should unprotect the old one and protect the new one but you need to
watch the order of those operations.

This approach is not a very good one.  If you really need to grow an
array it is better to allocate and reallocate the memory within your C
code using calloc and realloc then, at the end of the calculations,
allocate an R matrix and copy the results over.

Also, you haven't said whether you are growing the matrix by row or by
column or both.  If you are adding rows then you can't just reallocate
storage because R stores matrices in column-major order. The positions
of the elements in a matrix with n+1 rows are different from those in
a matrix with n rows.

> Is it necessary to reallocate or is R handling
> the memory management for the matrix that I have allocated?
>
> Regards
>
> Abhijit Bera
>
>        [[alternative HTML version deleted]]
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Confusion regarding allocating Matrices.

[Douglas Bates]

On Fri, Oct 23, 2009 at 9:23 AM, Douglas Bates  wrote:
> On Fri, Oct 23, 2009 at 8:39 AM, Abhijit Bera  wrote:
>> Hi
>>
>> I'm having slight confusion.
>
> Indeed.
>
>> I plan to grow/realloc a matrix depending on the data available in a C
>> program.
>
>> Here is what I'm tried to do:
>
>> Data=allocMatrix(REALSXP,3,4);
>> SEXP Data;
>
> Those lines should be in the other order, shouldn't they?
>
> Also, you need to PROTECT Data or bad things will happen.
>
>> REAL(Data)[8]=0.001123;
>> REAL(Data)[20]=0.001125;
>> printf("%f %f\n\n\n\n",REAL(Data)[8],REAL(Data)[20]);

And I forgot to mention, it is not a good idea to write REAL(Data)
many times like this.  REAL is a function, not a macro and you are
calling the same function over and over again unnecessarily.  It is
better to write

double *dat = REAL(Data);

and use the dat pointer instead of REAL(Data).

>> Here is my confusion:
>
>> Do I always require to allocate the exact number of data elements in a R
>> Matrix?
>
> Yes.
>
>> In the above code segment I have clearly exceeded the number of
>> elements that have been allocated but my program doesn't crash.
>
> Remember that when programming in C you have a lot of rope with which
> to hang yourself.   You have corrupted a memory location beyond that
> allocated to the array but nothing bad has happened  - yet.
>
>> I don't find any specific R functions for reallocation incase my data set
>> grows. How do I reallocate?
>
> You allocate a new matrix, copy the contents of the current matrix to
> the new matrix, then release the old one.  It gets tricky in that you
> should unprotect the old one and protect the new one but you need to
> watch the order of those operations.
>
> This approach is not a very good one.  If you really need to grow an
> array it is better to allocate and reallocate the memory within your C
> code using calloc and realloc then, at the end of the calculations,
> allocate an R matrix and copy the results over.
>
> Also, you haven't said whether you are growing the matrix by row or by
> column or both.  If you are adding rows then you can't just reallocate
> storage because R stores matrices in column-major order. The positions
> of the elements in a matrix with n+1 rows are different from those in
> a matrix with n rows.
>
>> Is it necessary to reallocate or is R handling
>> the memory management for the matrix that I have allocated?
>>
>> Regards
>>
>> Abhijit Bera
>>
>>        [[alternative HTML version deleted]]
>>
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Confusion regarding allocating Matrices.

[Douglas Bates]

On Fri, Oct 23, 2009 at 2:02 PM, Abhijit Bera  wrote:
> Sorry, I made a mistake while writing the code. The declaration of Data
> should have been first.

> I still have some doubts:

Because you are making some sweeping and incorrect assumptions about
the way that the internals of R operate.  R allows for arrays to be
dynamically resized but this is accomplished internally by allocating
new storage, copying the current contents to this new location and
installing the values of the new elements.  It is an expensive
operation, which is why it is discouraged.

Your design is deeply flawed.  Go back to the drawing board.

> When you say calloc and realloc are you talking about R's C interface Calloc
> and Realloc or the regular calloc and realloc?

Either one.

> I want to feed data directly into a R matrix and grow it as required. So at
> one time I might have 100 rows coming in from a data source. The next time I
> might have 200 rows coming in from a data source. I want to be able to
> expand the R-matrix instead of creating a regular C float matrix and then
> make an R-matrix based on the new size. I just want to have one R object and
> be able to expand it's size dynamically.

R stores floating-point numbers as the C data type double, not float.
It may seem pedantic to point out distinctions like that but not when
you are writing programs.  Compilers are the ultimate pedants - they
are real sticklers for getting the details right.

As I said, it just doesn't work the way that you think it does.  The
fact that there is an R object with a certain name before and after an
operation doesn't mean it is the same R object.

> I was reading the language specs. It says that one could declare an object
> in R like this:
>
> m=matrix(nrows=10,ncols=10)
>
> and then one could assign
>
> m[101]=1.00
>
> to expand the object.
>
> but this has one problem when I do a
>
> dim(m)
>
> I get
>
> NULL instead of 10 10
>
> So what is happening here?
>
>
> I am aware that R matrices are stored in column major order.
>
> Thanks for the tip on using float *dat= REAL(Data);
>
> Regards
>
> Abhijit Bera
>
>
>
> On Fri, Oct 23, 2009 at 7:27 PM, Douglas Bates  wrote:
>>
>> On Fri, Oct 23, 2009 at 9:23 AM, Douglas Bates 
>> wrote:
>> > On Fri, Oct 23, 2009 at 8:39 AM, Abhijit Bera 
>> > wrote:
>> >> Hi
>> >>
>> >> I'm having slight confusion.
>> >
>> > Indeed.
>> >
>> >> I plan to grow/realloc a matrix depending on the data available in a C
>> >> program.
>> >
>> >> Here is what I'm tried to do:
>> >
>> >> Data=allocMatrix(REALSXP,3,4);
>> >> SEXP Data;
>> >
>> > Those lines should be in the other order, shouldn't they?
>> >
>> > Also, you need to PROTECT Data or bad things will happen.
>> >
>> >> REAL(Data)[8]=0.001123;
>> >> REAL(Data)[20]=0.001125;
>> >> printf("%f %f\n\n\n\n",REAL(Data)[8],REAL(Data)[20]);
>>
>> And I forgot to mention, it is not a good idea to write REAL(Data)
>> many times like this.  REAL is a function, not a macro and you are
>> calling the same function over and over again unnecessarily.  It is
>> better to write
>>
>> double *dat = REAL(Data);
>>
>> and use the dat pointer instead of REAL(Data).
>>
>> >> Here is my confusion:
>> >
>> >> Do I always require to allocate the exact number of data elements in a
>> >> R
>> >> Matrix?
>> >
>> > Yes.
>> >
>> >> In the above code segment I have clearly exceeded the number of
>> >> elements that have been allocated but my program doesn't crash.
>> >
>> > Remember that when programming in C you have a lot of rope with which
>> > to hang yourself.   You have corrupted a memory location beyond that
>> > allocated to the array but nothing bad has happened  - yet.
>> >
>> >> I don't find any specific R functions for reallocation incase my data
>> >> set
>> >> grows. How do I reallocate?
>> >
>> > You allocate a new matrix, copy the contents of the current matrix to
>> > the new matrix, then release the old one.  It gets tricky in that you
>> > should unprotect the old one and protect the new one but you need to
>> > watch the order of those operations.
>> >
>> > This approach is not a very good one.  If you really need to grow an
>> > array it is better to allocate and reallocate the memory within your C
>&g

Re: [Rd] Summary methods

[Douglas Bates]

On Sun, Nov 8, 2009 at 2:26 PM, Doran, Harold  wrote:
> I've defined the following for objects of a class called jml
>
> summary.jml <- function(object, ...){
>        tab <- cbind(Estimate = coef(object),
>                        StdError = object$se,
>                        Infit = object$Infit,
>                        Outfit = object$Outfit)
>        res <- list(call = object$call, coefficients = tab,
>                        N = nrow(object$Data), iter = object$Iterations)
>        class(res) <- "summary.jml"
>        res
> }
>
> print.summary.jml <- function(x, ...){
>   cat("Call:\n")
>   print(x$call)
>   cat("\n")
>   cat("Number of iterations to completion:", x$iter, "\n")
>   cat("Number of individuals used in estimation:", x$N, "\n")
>   cat("\n")
>   printCoefmat(x$coefficients)
>   }
>
> Use of the methods on a fitted jml object yields:
>
>> summary(aa)
> Call:
> jml2.formula(formula = ~V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 +
>    V9 + V10, data = itemDat, bias = F)
>
> Number of iterations to completion: 6
> Number of individuals used in estimation: 486
>
>                 StdError     Infit Outfit
>  [1,] -0.380346  0.103002  1.007466 0.9935
>  [2,]  0.025939  0.104052  1.003050 1.0083
>  [3,]  2.563784  0.171174  0.941453 0.9414
>  [4,] -2.930519  0.156923  1.010786 1.0515
>  [5,]  1.139241  0.118932  0.978101 1.1424
>  [6,] -1.461751  0.111563  1.030612 1.2709
>  [7,]  0.486202  0.107986  1.008374 1.0394
>  [8,] -0.497102  0.103117  0.961431 0.9012
>  [9,] -0.486478  0.103099  1.001752 0.9829
> [10,]  1.541029  0.129214  1.010011 0.9150
>
> Two questions. First, why is the name of the first column empty instead of 
> "Estimate" as I have specified in the summary method?

Because you are using cbind to create the table.  Use data.frame
instead.  I think that will also help with the alignment issue.

> Second, I am struggling to get the row names of the coefficients to align 
> with the formula call. For example, instead of
>
>                 StdError     Infit Outfit
>  [1,] -0.380346  0.103002  1.007466 0.9935
>
> I would prefer
>
>                 StdError     Infit Outfit
>  V1 -0.380346  0.103002  1.007466 0.9935
>
> This also occurs in my print method
>
> print.jml <- function(x, digits = 2, ...){
>   cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
>   cat("Coefficients:\n")
>                print.default(format(coef(x), digits = digits), print.gap=2,
>                quote = FALSE)
>   invisible(x)
>   }
>
> Which produces
>
>> win
> Call:
> jml2.default(dat = itemDat[, 1:10])
>
> Coefficients:
>             [,1]
>  [1,] -0.38034638
>  [2,]  0.02593937
>  [3,]  2.56378422
>  [4,] -2.93051899
>  [5,]  1.13924076
>  [6,] -1.46175119
>  [7,]  0.48620247
>  [8,] -0.49710150
>  [9,] -0.48647770
> [10,]  1.54102895
>
> Thank you
> Harold
>
>> sessionInfo()
> R version 2.10.0 (2009-10-26)
> i386-pc-mingw32
>
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
> [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
> [5] LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] MASS_7.3-3         lme4_0.999375-32   Matrix_0.999375-31 lattice_0.17-26
> [5] MiscPsycho_1.4     statmod_1.4.1
>
> loaded via a namespace (and not attached):
> [1] grid_2.10.0  plink_1.2-2  tools_2.10.0
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] snark of the day: new "features" (??) in Mathematica 7

[Douglas Bates]

On Mon, Jan 11, 2010 at 10:34 AM, Ben Bolker  wrote:

>  who wants to write rgl code to do these?

> http://tinyurl.com/yzojfn2
> http://tinyurl.com/ylrz2p8

>  :-)

I think that on this one we should follow Nancy Reagan's advice and
"Just say no".

I wonder if they have read Tufte's descriptions of "chart junk".  If
so they missed the point that chart junk isn't a good thing.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] unable to compile Recommended packages

[Douglas Bates]

Did you run

$RHOME/src/tools/rsync-recommended

to obtain the most recent versions of the recommended packages before
trying to compile them?

On Fri, Mar 5, 2010 at 3:01 PM, Alex Bryant  wrote:
> Hi folks,  I'm having a problem with installing  R on Solaris.  Has anyone 
> seen a similar issue?  I don't find any hits on the search engines.  Thanks.
>
>
> Environment:
> SunOS  5.10 Generic_118822-25 sun4u sparc SUNW,Sun-Fire-280R
>
> During the installation/compilation of R-2.10.1 I  get the following error:
>
> begin installing recommended package boot
> ERROR: cannot extract package from 'boot.tgz'
> *** Error code 1
> The following command caused the error:
> MAKE="make" R_LIBS= ../../../bin/R CMD INSTALL --no-lock -l 
> "../../../library" boot.tgz > boot.ts.out 2>&1 || (cat boot.ts.out && exit 1)
> make: Fatal error: Command failed for target `boot.ts'
> Current working directory /home/ireview/R-2.10.1/src/library/Recommended
> *** Error code 1
> The following command caused the error:
> make stamp-recommended
> make: Fatal error: Command failed for target `recommended-packages'
> Current working directory /home/ireview/R-2.10.1/src/library/Recommended
> *** Error code 1
> The following command caused the error:
> (cd src/library/Recommended && make)
> make: Fatal error: Command failed for target `stamp-recommended'
>
>
> -Alex
>
>
> //
> // Alex Bryant
> // Software Developer
> // Integrated Clinical Systems, Inc.
> // 908-996-7208
>
>
> 
> Confidentiality Note: This e-mail, and any attachment to...{{dropped:13}}
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Defining a method in two packages

[Douglas Bates]

I think a simpler solution, which I have implemented in lme4 but not
yet released is to have

importFrom(nlme, ranef, fixef)

in the NAMESPACE file of packages that implement methods for those
generics (and, of course, add nlme to the Imports: specification in
the DESCRIPTION file).  As nlme is a required package I don't think
this is too much of a burden.

On Tue, Mar 9, 2010 at 8:52 PM, Kevin Coombes  wrote:
> Wouldn't it make sense to simply create a "ranef" package whose only role in
> the universe is to create the generic function that lme4, coxme, and anyone
> else who needs it could just import, without getting tons of additional and
> (depending on the application) irrelevant code?
>
> Best,
>   Kevin
>
> Uwe Ligges wrote:
>>
>>
>> On 08.03.2010 17:16, Terry Therneau wrote:
>>>
>>> Brian&  Uwe,
>>>   Thanks for responding.  Let me see if I can refine the query and move
>>> towards a solution.
>>>
 From Uwe:
 Of course, after loading lme4, you can still use the ranef from coxme:
 coxme::ranef(fit)
>>>
>>> Of course, but I'm interested in other users (as well as myself) and
>>> prefer to avoid the 'secret handshake' form of a call.
>>>
 In your own package, you could simply import the generic from coxme.
>>>
>>> I don't understand this.
>>
>> You could import the generic from the other package and define your won
>> methods for it in order to make dispatching work correctly.
>>
>>
>>
 From Brian:
 My solution would though be NOT to reuse a name that is already
 established in another package (nlme has used it for many years).
 The design problem is that generic foo() in package B might have
 nothing to do with foo() in package A.  When it does, we expect B ...
>>>
>>> I disagree completely.  It is precisely because of nlme and lmer
>>> prominence that I want to reprise their methods: my users have a much
>>> better chance of remembering how to do things.  If I followed this logic
>>> to its conclusion one should never define a print() method because it
>>> might conflict with the base definition.
>>>   The consequence is that I am under obligation to NOT make my method
>>> something different than Doug's, if I want to satisfy the goal of user
>>> level consistency.  Several aspects of coxme purposefully mimic lmer,
>>> even in cases (such as print.coxme) where his layout is not precisely
>>> what I would have chosen.
>>
>> Then please folow my suggestion and import the generic from the packages
>> mentioned above in your namespace. Then you could extend it by your own
>> methods wihtout having to define another generic of the same name and avoid
>> the conflicts.
>>
>>
>>>   I really do not want to require lme4 just to pick up the methods
>>> definition.  It's a huge package, and there is no code in common.  Both
>>> packages work very hard to be efficient via sparse matrix methods, but
>>> the actual details are completely different due to the mathematical
>>> structure of our underlying likelihoods.  Use of both in the same
>>> analysis would be rare, so my issue won't be common.
>>
>> Well, then things become complicated if not impossible.
>>
>>
 The situation can be alleviated by making S3 methods visible.  Thus if
 coxme exported coxme:::ranef.coxme and lme4 had a default method
>>>
 ranef<- function (object, ...) UseMethod("ranef")
>>>
>>>  I have no objection to exporting my method.  If a joint change to lme4
>>> and coxme is the best solution, I will take the discussion off line with
>>> Doug.  Is this the best way forward?
>>
>> I think so.
>>
>> Best wishes,
>> uwe
>>
>>
>>>
>>> Terry
>>>
>>>
>>>
>>>
>>
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Lapack, determinant, multivariate normal density, solution to linear system, C language

[Douglas Bates]

On Mon, Apr 12, 2010 at 10:27 PM, shotwelm  wrote:
> r-devel list,
>
> I have recently written an R package that solves a linear least squares
> problem, and computes the multivariate normal density function.

For both of those applications you can use a Cholesky decomposition of
the symmetric matrix.  If the Cholesky decomposition fails then you
have a singular least squares problem or a singular
variance-covariance matrix for your multivariate normal density
function.

Have you tried comparing the speed of your code to

prod(diag(chol(mm))^2

or, probably better, is to use the logarithm of the determinant

2 * sum(log(diag(chol(mm)))

If you use the Matrix package class dpoMatrix to solve the linear
system it will cache the results of the Cholesky decomposition when
solving the system so later evaluation of the determinant will be very
fast - although I suspect you would need to be working with matrices
of sizes in the hundreds or doing the same operation thousands of
times before you would notice a difference.

If you really insist on doing this in compiled code you just need to
call F77_CALL(dpotrf) then accumulate the product of the diagonal
elements of the resulting factor.

You could use packed storage but the slight advantage in memory usage
(at best, 1/2 of the full storage usage) is not worth the pain of
writing code to navigate the packed storage locations.

> The bulk
> of the code is written in C, with interfacing code to the BLAS and
> Lapack libraries. The motivation here is speed. I ran into a problem
> computing the determinant of a symmetric matrix in packed storage.
> Apparently, there are no explicit routines for this as part of Lapack.
> While there IS an explicit routine for this in Linpack, I did not want
> to use the older library. Also, right before I needed the determinant of
> the matrix A, I had used the Lapack routine dspsv to solve the linear
> system Ax=b, which does much of the work of computing a determinant
> also. In fact, the solution I came up with involves the output of this
> routine (which might be obvious to Lapack designers, but not me)
>
> My modest Googleing turned up very little unique material (as is typical
> with BLAS/Lapack/Linpack queries). Hence, I am writing the r-devel list
> partly to document the solution I've come up with, but mainly to elicit
> additional wisdom from seasoned R programmers.
>
> My solution to the problem is illustrated in the appended discussion and
> C code. Thanks for your input.
>
> -Matt Shotwell
>
> --
>
> The Lapack routine dspsv solves the linear system of equations Ax=b,
> where A is a symmetric matrix in packed storage format. The dspsv
> function performs the factorization A=UDU', where U is a unitriangular
> matrix and D is a block diagonal matrix where the blocks are of
> dimension 1x1 or 2x2. In addition to the solution for x, the dspsv
> function also returns the matrices U and D. The matrix D may then be
> used to compute the determinant of A. Recall from linear algebra that
> det(A) = det(UDU') = det(U)det(D)det(U'). Since U is unitriangular,
> det(U) = 1. The determinant of D is the product of the determinants of
> the diagonal blocks. If a diagonal block is of dimension 1x1, then the
> determinant of the block is simply the value of the single element in
> the block. If the diagonal block is of dimension 2x2 then the
> determinant of the block may be computed according to the well-known
> formula b11*b22-b12*b21, where bij is the value in the i'th row and j'th
> column of the block.
>
>  int i, q, info, *ipiv, one = 1;
>  double *b, *A, *D, det;
>
>  /*
>  ** A and D are upper triangular matrices in packed storage
>  ** A[] = a00, a01, a11, a02, a12, a22, a03, a13, a23, a33, ...
>  ** use the following macro to address the element in the
>  ** i'th row and j'th column for i <= j
>  */
>  #define UMAT(i, j) (i + j * ( j + 1 ) / 2)
>
>  /*
>  ** additional code should be here
>  ** - set q
>  ** - allocate ipiv...
>  ** - allocate and fill A and b...
>  */
>
>  /*
>  ** solve Ax=b using A=UDU' factorization where
>  ** A represents a qxq matrix, b a 1xq vector.
>  ** dspsv outputs the elements of the matrix D
>  ** is upper triangular packed storage
>  ** in the memory addressed by A. That is, A is
>  ** replaced by D when dspsv returns.
>  */
>  F77_CALL(dspsv)("U", &q, &one, A, ipiv, b, &q, &info);
>  if( info > 0 ) { /*issue warning, system is singular*/ }
>  if( info < 0 ) { /*issue error, invalid argument*/ }
>
>  /*
>  ** compute the determinant det = det(A)
>  ** if ipiv[i] > 0, then D(i,i) is a 1x1 block diagonal
>  ** if ipiv[i] = ipiv[i-1] < 0, then D(i-1,i-1),
>  ** D(i-1,i), and D(i,i) form the upper triangle
>  ** of a 2x2 block diagonal
>  */
>  D = A;
>  det = 1.0;
>  for( i = 0; i < q; i++ ) {
>    if( ipiv[ i ] > 0 ) {
>      det *= D[ UMAT(i,i) ];
>    } else if( i > 0 && ipiv[ i ] < 0 && ipiv[ i-1 ] == ipiv[ i ] ) {
>      det *= D[ UMAT(i,i) ] * D[ UMAT(i-1,i-1) ] -\
>             D[ UM

[Rd] The XO laptop from the One Laptop Per Child (OLPC) program

[Douglas Bates]

There was recently a question on the R-help list about the eee pc.  I
had a related question about the XO laptop from OLPC (laptop.org).
Has anyone looked at the development environment sufficiently to
determine if it would be possible to create an executable image for R?
 The laptop itself only supports Python, Javascript, etc. but it is
running a real Linux operating system.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] S3 vs S4 for a simple package

[Douglas Bates]

On Jan 7, 2008 1:34 PM, John Chambers <[EMAIL PROTECTED]> wrote:
> Prof Brian Ripley wrote:
> > On Mon, 7 Jan 2008, Robin Hankin wrote:
> >
> >
> >> I am writing a package and need to decide whether to use S3 or S4.
> >>
> >> I have a single class, "multipol"; this needs methods for "[" and "[<-"
> >> and I also need a print (or show) method and methods for arithmetic +-
> >> */^.
> >>
> >> In S4, an object of class "multipol" has one slot that holds an array.
> >>
> >> Objects of class "multipol" require specific arithmetic operations;
> >> a,b being
> >> multipols means that a+b and a*b are defined in peculiar ways
> >> that make sense in the context of the package. I can also add and
> >> multiply
> >> by scalars (vectors of length one).
> >>
> One thing you cannot do in S3 is to have methods that depend on anything
> but the first argument.  Do you want something sensible for 1 + a when
> a is a "multipol"? The default call to the primitive version may or may
> not give you what you want.

I agree with John.  If method dispatch on multiple arguments is
potentially useful to you, and it sounds as if it is, then S4 is
clearly the winner.  In the Matrix package, for example, determining a
good method for evaluating A %*% B or solve(A, B), where A and B are
Matrix objects, or matrix objects would be extremely difficult if we
were not using S4 classes and methods.

> >> My impression is that S3 is perfectly adequate for this task, although
> >> I've not yet finalized the coding.
> >>
> >> S4 seems to be "overkill" for such a simple system.
> >>
> >> Can anyone give me some motivation for persisting with S4?
> >>
> >> Or indeed reassure me that S3 is a good design decision?
> >>
> >
> > Does performance matter?: S4 dispatch is many times slower than S3
> > dispatch for such functions. (It is several times slower in general, but
> > the difference is particularly marked for primitives.)
> >
> Well, the question is whether performance of _method dispatch_ matters,
> which it tends not to in many cases. And it would be good to have some
> data to clarify "many times slower".  Generally, looking up inherited
> methods is likely to take a while, but only the first time does R code
> need to work out the inheritance. On repeated calls with the same
> signature, dispatch should be basically a lookup in an environment.

I can imagine circumstances in which the time spent on method dispatch
would be important but I haven't encountered them and we have been
working with S4 classes and methods in the Matrix and lme4 packages
for a long time now.  Much more important to me is the amount of time
that I spend designing and debugging code and that is where S4 is, for
me, the clear winner.  Because of the formal class definitions and the
validation of objects  I can write C code assuming that the contents
of particular slots in an object have the desired class.  When I end
up dealing with S3 objects not S4 objects I find it very clunky
because there is no formal class definition and, if I am to play it
safe, I should always check for the existence of certain components in
the object then check or coerce the type then check the length, etc.
With S4 classes I only need to do that once in a validate method.

With S4 I feel that I can use a single class system in my R code and
in the underlying C code.  I don't need to create exotic structs or
C++ classes in the compiled code because the S4 class system in R is
enforcing the class definition for me.

I just read Tim Hesterberg's message from later in this thread and I
agree with him that the formal definition of classes makes it awkward
to revise the slots.  Fortunately, I have a very understanding user
base for the lme4 package and they tolerate my changing class
definitions between releases because I manage to convince them somehow
that the package is improved by the changes.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] 0.45<0.45 = TRUE (PR#10744)

[Douglas Bates]

In the "R Programming Style" thread on R-help Ronald Rau gave a list
of aphorisms from "Elements of Programming Style" by Kernighan and
Plauger.  These include

10.0 times 0.1 is hardly ever 1.0

I think that should be included in FAQ 7.31

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R CMD check should check date in description

[Douglas Bates]

On Fri, Apr 4, 2008 at 2:26 PM, Kurt Hornik <[EMAIL PROTECTED]> wrote:
> > hadley wickham writes:
>
>  > I'm always forgetting to update the date in DESCRIPTION.  Would it be
>  > possible to add a warning to R CMD check if it's old?
>
>  I recently thought about this.  I see several issues.
>
>  * How can we determine if it is "old"?  Relative to the time when the
>   package was uploaded to a repository?
>
>  * Some developers might actually want a different date for a variety of
>   reasons ...
>
>  * What we currently say in R-exts is
>
>  The optional `Date' field gives the release date of the current
>  version of the package.  It is strongly recommended to use the
>  -mm-dd format conforming to the ISO standard.
>
>   Many packages do not comply with the latter (but I have some code to
>   sanitize most of these), and "release date" may be a moving target.
>
>  The best that I could think of is to teach R CMD build to *add* a Date
>  field if there was none.

I agree.  I think that the Date field in the DESCRIPTION file should
reflect the date and time zone at which the package was built or some
other date of the package maintainer's choice.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Reaching R-forge - tunnel ssh through ssh?

[Douglas Bates]

This question is only related to R-devel in that I need to work out
how to reach R-forge.R-project.org to update and commit files of R
packages for which I am a developer.

My desktop computer at work runs Ubuntu 8.04.  Because this is not a
supported operating system this computer is behind a firewall.  I can
use a proxy server for http, https and rsync connections.  However,
the proxy server does not allow outgoing (or incoming) ssh
connections.  I can use ssh to log in to a local computer that is not
behind the firewall and then use svn+ssh to contact R-forge.  I have
managed to get myself into a bind with that scheme, however.  I added
and deleted files on the checked out copy on the shared file system
using my desktop and then could not check in the changes from the
"visible" machine because of svn version skew.  The "visible" machine
has an out-of-date version of svn.

Does anyone know how to tunnel ssh over ssh and, further, how I would
specify that in the svn checkin?  If my desktop computer is A, the
"visible" computer is B and the R-forge server is C then I want to use
ssh to set up a tunnel between A and B so that A can reach C for an
svn+ssh checkin.

Off-list replies would be welcome.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] "for" loop wiht S4 Objects

[Douglas Bates]

On Fri, Jan 23, 2009 at 1:25 AM, Biczok Rudolf
 wrote:
> Hi all,

> I'm working with the S4-Class system and I have a little problem with

> Implementing iteration functionality in my S4 class but it don't work:

>> setClass("foo",representation(bar="list"))
>> x <- new("foo",bar=list(1,2,3))
>>for(e in x) cat(e)
> invalid type/length (S4/1) in vector allocation

I'm not sure why you would expect a for loop to work there.  Wouldn't
a function in the "apply" family be more natural?  At least in my code
it is much more common to use lapply to loop over the elements of a
list than it is to use 'for'.

I once suggested to John Chambers that it would be a good idea to have
a way of applying a function to the slots of an S4 object and he
persuaded me otherwise.  However, you can generate something like that
with

lapply(slotNames(x), function(nm) cat(slot(x, nm)))

or create a list out of the S4 object

S4toL <- function(x) {nms <- slotNames(x); names(nms) <- nms;
lapply(nms, slot, object = x)}

then use lapply on the result of S4toL

>
>
>
>
>
> But when I extend from a primitive vector it works:
>
>
>
>> setClass("foo",contains="list")
>
>
>
>> x <- new("foo",.Data=list(1,2,3))
>
>
>
>>for(e in x) cat(e)
>
> 123
>
>
>
> This is ok, but is there any other way to do this (for e.g. with a
> generic function)?
>
>
>
> Thanks,
>
> Rudolf
>
>
>[[alternative HTML version deleted]]
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] [SoC09-Info] An IPopt interface for R

[Douglas Bates]

There have been several messages on R-devel mentioning the interior
point optimization software Ipopt, https://projects.coin-op/Ipopt/.
This C++ library is released under a license called the Common Public
License.

I have two questions that readers of R-devel may be able to answer.

1) Would creating an Ipopt interface for R be duplicating existing
efforts?  That is, has someone already done so or is already working
on this project.

2) Could a package incorporating Ipopt be released under some version
of the GPL or LGPL?

Assuming that the answers to those questions are no and yes, I would
be willing to mentor a student in a Google Summer of Code project to
create such a package.  However, I don't have a whole lot of time to
do so and would not object at all if another potential mentor, or
perhaps a co-mentor, were to volunteer.  Failing that I will create an
application, including a programming exercise.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R tutorial

[Douglas Bates]

On Sun, Feb 22, 2009 at 7:44 AM, Christophe Dutang  wrote:
> Dear all,

> I have just found a 'good' tutorial R for datamining. I think it should be
> on the contributed docs.
> http://cran.r-project.org/other-docs.html

> Here is the link
> http://www.liaad.up.pt/~ltorgo/DataMiningWithR/

> What do you think?

A quick glance indicates that the last updates were in 2003.   Most
introductory material from then would still apply but details may have
changed and often there are now better ways of doing things than there
were then.

Do you know if the author has plans for updating the book at all?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [SoC09-Info] An IPopt interface for R

[Douglas Bates]

On Sun, Feb 22, 2009 at 3:05 PM, Douglas Bates  wrote:
> There have been several messages on R-devel mentioning the interior
> point optimization software Ipopt, https://projects.coin-op/Ipopt/.
> This C++ library is released under a license called the Common Public
> License.
>
> I have two questions that readers of R-devel may be able to answer.
>
> 1) Would creating an Ipopt interface for R be duplicating existing
> efforts?  That is, has someone already done so or is already working
> on this project.
>
> 2) Could a package incorporating Ipopt be released under some version
> of the GPL or LGPL?
>
> Assuming that the answers to those questions are no and yes, I would
> be willing to mentor a student in a Google Summer of Code project to
> create such a package.  However, I don't have a whole lot of time to
> do so and would not object at all if another potential mentor, or
> perhaps a co-mentor, were to volunteer.  Failing that I will create an
> application, including a programming exercise.

Having looked in more detail at the Ipopt sources and installation
instructions I think there will be a problem.  While Ipopt itself is
covered by the Common Public License it requires a sparse indefinite
solver such as MUMPS, Pardiso or one of the Harwell Subroutine
Libraries, http://www.coin-or.org/Ipopt/documentation/node6.html .  As
far as I can see, none of these are covered by a license that allows
redistribution so they can't be included in a package on CRAN or, I
think, in an SoC project.

It's too bad.  I am encouraged by projects like COIN, www.coin-or.org,
and Ipopt that do recognize the importance of open source software
covered by a suitable license but so much of the numerical analysis
community still misses the boat on that one.

By the way, the URL in the quoted message was garbled.  It should have been

https://projects.coin-or.org/Ipopt

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] 1-based arrays and copying vectors

[Douglas Bates]

On 6/19/05, Rob Lopaka Lee <[EMAIL PROTECTED]> wrote:
> I'm interfacing to C code that uses 1-based indexing on arrays -- it
> ignores the zeroth element.  Thus, input vectors from R must be moved up
> one, and output arrays must be moved down one.
> 
> What is the best way to deal with this using R internal code?

If the C code can be relied upon to ignore the zero'th element of the
array then write a wrapper that gets the array, say int v[], from R
and passes &v[-1] to your C code.  That has the effect of shifting all
the addresses back by one position.
 
> My current approach is:
> 
> For an input R vector of length n, allocate a new vector(v) of length n+1
> and copy input into v[1] to v[1+n].  Call with new vector.
> 
> For an output array(a) of length n, allocate a new vector of length n-1
> and copy a[1] to a[n] into v[0] to v[n-1].
> 
> If this is the best approach, is there an idiom for copying vectors?

In C you could use the memcpy function or the Memcpy macro defined in
.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Lapack calls from fortran

[Douglas Bates]

On 6/22/05, Paul Gilbert <[EMAIL PROTECTED]> wrote:
> I am trying to call the lapack routine dgesv in the R distribution from
> the fortran code for a package, but when I dyn.load("dse1.so") I get an
> error about undefined symbol: dgesv_
> 
> I thought the proper incantation was
> 
> MAKEFLAG="FLIBS=/path to/R-2.1.1/libRlapack.so" R CMD SHLIB dse1.f
> 
> but apparantly not. How does one do this?

See section 1.2.1 of "Writing R Extensions".  You should add a file
called Makevars to the src directory and include the line

PKG_LIBS = ${LAPACK_LIBS} ${BLAS_LIBS}

in that file

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Running ./tools/rsync-recommended through a proxy

[Douglas Bates]

My computers at my office will no longer be able to connect directly
to web sites etc.  I will be going through a proxy server.  The
particular server is running squid on port 3128.

I have managed to configure web browsers, ssh, apt, svn and a whole
lot of other tools to use the proxy server but I haven't been able to
configure rsync.  My usual method of updating my copy of the R-devel
sources is via

cd my_R_devel_sources
svn up
./tools/rsync-recommended

cd my_R_build_directory
...

Can anyone offer suggestions on how to get rsync-recommended to work
through a proxy?  I have set

export RSYNC_PROXY="machine.name:3128"

which I understand from the documentation is the magic environment
variable (I don't give the name of the server explicitly because it is
an open proxy).  However, I still get

[EMAIL PROTECTED]:/usr/src/r-devel$ ./tools/rsync-recommended
bad response from proxy - HTTP/1.0 403 Forbidden
rsync: failed to connect to machine.name: Success (0)
rsync error: error in socket IO (code 10) at clientserver.c(94)
*** rsync failed to update Recommended files ***
Creating links

Is it likely that the proxy server is not passing connections to port 873?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] How difficult is it to wrap a large C++ library with R?

[Douglas Bates]

On 7/4/05, Bo Peng <[EMAIL PROTECTED]> wrote:
> > > * Wrap C++ class hierarchy. Virtual functions need to be supported.
> > > (SWIG can generate Python shadow classes that behave almost exactly
> > > like the underlying C++ classes)
> >
> > This is hard to do in R, because the R object model is quite different
> > from that of C++ (whereas Python's is more similar).
> 
> I guess I will have to learn R OOP in detail first. It would be
> terrible if I will have to re-write every C++ class in R (and without
> inheritance)!

What makes you think that there is no inheritance of classes in R?
 
> > > * Direct access to C++ data structure. For example, an object may keep
> > > a C array internally. I will need a method (maybe wrap this data into
> > > a R object) to read/write this array directly.
> >
> > That's not too hard provided you use C++ code to do the actual access.
> > That is, you write an R function that calls C++ code to do the work.
> > It's a lot harder if you want to keep it all in R, because it doesn't
> > understand C++ type definitions, alignment conventions, etc.
> 
> So this has to be done through functions. In Python, I can create an
> object to wrap C array that can be accessed just like regular list. I
> guess it would be difficult to create a new data type in R.
> 
> Anyway, where can I find the most complete information about
> implementation of R/OOP at C/C++ level? (Or should I write classes in
> R and call C++ functions for every member function?) I have not read
> http://cran.r-project.org/doc/manuals/R-exts.pdf in detail but my
> impression is that almost everything is done at a function (not
> object) level using .C().

Well the R object model is based on generic functions and methods are
functions.  The S4 class structure differs from the S3 structure in
that there is a formal definition of an S4 class but not for an S3
class.  However, even in S4, methods do not belong to a class. 
Methods are intermediate to classes and generic functions and are
associated with a signature for the arguments of the generic function.

If you want to be able to use S4 classes effectively you should use
the .Call interface to compiled code and not the .C interface.  If you
pass a classed object through the .Call interface to compiled code you
can extract or set the values of slots in the object through the
GET_SLOT and SET_SLOT macros.  Admittedly this is not as convenient
for the C++ programmer as being able to pass a pointer to an object
from a C++ class but it does make it easier to think of objects in R
and the corresponding objects in compiled code.  The Matrix package on
CRAN contains many examples of this use of this technique.

Another technique that you may want to consider is the use of
ExternalPointers.  If you examine the C source file lmer.c from the
aforementioned Matrix package you will see where I create a pointer to
a C struct in the compiled code and pass it back to R so it can later
be passed to the C code and I can extract the struct.  The same could
be done for an instance of a C++ class.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] calling fortran functions CHOL and DPOTRF form Fortran

[Douglas Bates]

On 7/5/05, Gilles GUILLOT <[EMAIL PROTECTED]> wrote:
> Hi all,
> 
> I'm working out some Fortran code for which
> I want to compute the Choleski decomposition of a covariance matrix
> in Fortran.
> 
> I tried to do it by two methods :
> 
> 1) Calling the lapack function DPOTRF.
> I can see the source code and check that my call is correct,
> but it does not compile with:
> system("R CMD SHLIB ~/main.f")
> dyn.load("~/main.so")
> 
> I get:
> Error in dyn.load(x, as.logical(local), as.logical(now)) :
> unable to load shared library
> [...]
>  undefined symbol: dpotrf_
> 
> Could anyone point out how to compile with such a call to lapack?

Make a small package from your source code then follow the instructions in

http://tolstoy.newcastle.edu.au/R/devel/05/06/1382.html

> 
> 2) Calling the Fortran function CHOL
> The previous combination of R CMD SHLIB  and dyn.load
> compiles fine but  I would like to see the source code of this function
> to check if my call is correct. I've not been able to find it in  the R source
> distribution.
> Could anyone point the place of that file?

${RSRC}/src/appl/chol.f

> 
> Thanks in advance,
> 
> Gilles
> 
> 
> --
> _
> Gilles GUILLOT
> 
> INRA -Département Mathématiques et Informatique Appliquées
> Unité Mixte de Recherche INRA - INAPG - ENGREF
> 
> Institut National Agronomique de Paris-Grignon
> 16, rue Claude Bernard
> 75231 Paris cedex 5 France
> 
> phone 33 1 44 08 18 42
> fax   33 1 44 08 16 66
> http://www.inapg.fr/ens_rech/mathinfo/personnel/guillot/welcome.html
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Spurious output like in R-devel

[Douglas Bates]

I have noticed spurious message of the form

when developing code for the Matrix package and testing under R-devel.
 These messages are not present when testing under R-2.1.1

I have not reported this because I didn't know if it was caused by my
code in the Matrix package or some other code.  Today I checked the
newly uploaded vcd package under R-devel and got the same symptoms.

Try

library('vcd')
example(mosaic')

to reproduce the problem.  If you want to see the same problem in the
Matrix package run

R-devel CMD check Matrix

and look at the output in Matrix.Rcheck/tests/lmer.Rout

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Spurious output like in R-devel

[Douglas Bates]

On 7/12/05, Duncan Temple Lang <[EMAIL PROTECTED]> wrote:
> Hi Doug.
> 
>   I noticed this also after the recent change to the
> code to handle missing PACKAGE arguments for
> .C/.Call/.Fortran routines.
>   There is a line in the dotcode.c (1510) that
> has a call to Rf_PrintValue(env)  and that is
> what is causing it.
> 
> I haven't had a chance to get to this for the last week or so but
> am doing so now.

Thanks Duncan.  In my code it turns out that the two .Call calls that
were producing the spurious output did not have PACKAGE = "Matrix".
Once I fixed that I no longer got the spurious output.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] extracting and manipulating components from nlme and gnls objects

[Douglas Bates]

On 7/19/05, Nicholas Lewin-Koh <[EMAIL PROTECTED]> wrote:
> Hello,
> I am writing a set of functions to do prediction and calibration
> intervals
> for at least the subset of selfstarting models if not some more general
> ones.
> 
> I need to be able to extract the varFunction from a fit object
> and evaluate it at a predicted point. Are there any examples around?
> 
> Also in a self start model, say SSlogis, how would I get the gradient
> at a point?

I think that the output of 

example(SSlogis)

should answer that question for you.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R_AllocatePtr

[Douglas Bates]

On 7/19/05, Paul Roebuck <[EMAIL PROTECTED]> wrote:
> Had been looking into Luke Tierney's R_AllocatePtr() and
> was left with a question about exactly when does R reclaim
> heap memory. Implication of 'simpleref.nw' is that one can
> allocate C data on the R heap, and as long as pointer object
> is alive, the data and pointer will remain valid. But it
> calls allocString() which is implemented using R_alloc().
> 
> 'Writing R Extensions" says R will reclaim memory allocated
> using R_alloc() at the end of the .Call().
> 
> Was the intent that its invocation be restricted to within
> a .Call()? Doesn't sound as though it could be used to return
> an SEXP from a .Call() with its data intact.

Sounds like you may want an R external pointer object, which is also
documented in simpleref.nw

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Reading a repeated fixed format

[Douglas Bates]

The Harwell-Boeing format for exchanging matrices is one of those
lovely legacy formats that is based on fixed-format  Fortran
specifications and 80 character records.  (Those of you who don't know
why they would be 80 characters instead of, say, 60 or 100 can ask one
of us old-timers some day and we'll tell you long, boring stories
about working with punched cards.)

Reading this format would take about 10 lines of R code if it were not
for the fact that it allows things like 40 two-digit integers to be
written as one 80 character record with no separators.  This actually
made sense to some people once upon a time.

I could use read.fwf or, better, use some of the code in the read.fwf
function to extract the strings that should have been separated and
convert them to numeric values but I have been trying to think if
there is a more clever way of doing this.  I know the number of
records and the number of elements to read and, if it would help, I
can assemble the records into one long text string.

Can anyone think of a vectorized way to extract successive substrings
of length k or, perhaps, a way to use regular expressions to insert a
blank after every k characters?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] include C functions from nmath in my own C functions

[Douglas Bates]

On 8/11/05, yyan liu <[EMAIL PROTECTED]> wrote:
>  Hi:
>   I followed the README in src/nmath/standalone/
>  to make the use the command "make shared" to make the
>  libRmath.so file. I also add the directories containg
>  libRmath.so to  LD_LIBRARY_PATH by using command
>  "export
> D_LIBRARY_PATH=$LD_LIBRARY_PATH:$/home/zhliu/Backup/R-2.0.1/src/nmath/standalon
> e
>  "
>  However, when I try to run the following codes by the
> command "gcc test.c -lRmath" on Linux Fedora Core 2,
>  /***/
> /* file name test.c */
>  #define MATHLIB_STANDALONE 1
>  #include 
> 
>  int
>  main()
>  {
>  /* something to force the library to be included */
> qnorm(0.7, 0.0, 1.0, 0, 0);
> return 0;
>  }
>  /**/
> 
>  the compiler gives me the following error message. It
> seems definitions of some R functions can not be found
> in the libRmath.so file. Anyone has any idea about
> this
> problem? Thank you very much!
> 
> 
>  /usr/local/lib/libRmath.so: undefined reference to
>  `expm1'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `log'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `sqrt'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `rint'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `cos'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `sin'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `pow'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `sinh'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `log10'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `exp'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `tan'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `log1p'
>  /usr/local/lib/libRmath.so: undefined reference to
>  `hypot'
>  collect2: ld returned 1 exit status
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

Add -lm to the end of the call to gcc.  You are missing functions from
the math library.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] An lgrind definition for the S language

[Douglas Bates]

I seem to recall discussion of an language definition file for S for
use with the lgrind utility but I can't find any trace of it in an R
Site Search.  The lgrind utility takes a file of code in a particular
programming language and prepares it for "pretty printing" in LaTeX. 
In my version the available language definitions are

$ lgrind -s
When specifying a language case is insignificant. You can use the
name of the language, or, where available, one of the synonyms in
parantheses. Thus the following are legal and mark Tcl/Tk, Pascal
and Fortran input, respectively:
   lgrind -ltcl/tk ...
   lgrind -lpaSCAL ...
   lgrind -lf ...
The list of languages currently available in your lgrindef file:
Ada   MLisp  (Emacs Mock Lisp)
Asm   SML/NJ  (ML)
Asm68 Scheme  (scm)
BASIC model
Batch  (bat)  Modula2  (mod2, m2)
C Pascal  (pas, p, bp)
C++  (CC) PERL  (pl)
csh   PostScript  (ps)
FORTRAN  (f77, f) PROLOG
Gnuplot   Python  (py)
Icon  RATFOR
IDL   RLaB
ISP   Russell
Java  SAS
Kimwitu++  (kimw) SDL
LaTeX sh
LDL   SICStus
Lex   src
Linda SQL
make  Tcl/Tk  (tcl, tk)
MASM  VisualBasic  (vbasic)
MATLABVMSasm
Mercury   yacc  (y)

Does anyone know of a similar facility for S code?

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Questions about calls and formulas

[Douglas Bates]

On 8/22/05, Erich Neuwirth <[EMAIL PROTECTED]> wrote:
> I am trying to adapt boxplot.formula (in graphics) to accept an
> additional parameter, weights.
> I already managed to adapt boxplot.default to do this.
> 
> boxplot.formula prepares the data for a call to boxplot.default and to
> achieve that does the following: It takes a formula like
> 
> x~g*h
> 
> as the first argument, and then by using
> 
> m <- match.call(expand.dots = FALSE)
> 
> saves the call. It transforms the call
> 
> m$na.action <- na.action # force use of default for this method
> m[[1]] <- as.name("model.frame")
> 
> and then  evaluates the modified call
> mf <- eval(m, parent.frame())
> 
> print(m)
> gives
> model.frame(formula = x ~ g * h)
> 
> Then it uses components of mf for the call to boxplot.default.
> 
> m has a component m$formula containing the parsed model formula.
> mode(m$formula) is "call".
> In our case, deparse(m$formula) gives a string representation of the
> formula: "x~g*h".
> I want to replace the response variable (in our case x) by the weights
> variable, which in the string expression can be done easily with
> strsplit and paste. Then I need to reconvert the modified string to a call.
> 
> So I create newmodelstring<-"weights~g*h" and try
> 
> m$formula<-as.call(parse(newmodelstring))
> 
> print(m)
> gives
> model.frame(formula = weights ~ g * h())
> 
> 
> When I try to evaluate the modified m this does not work. When I try to
> evaluate m with this modification I get
> 
> Error in model.frame(formula = weights ~ g * h()) :
> attempt to apply non-function
> 
> Is there a way to get rid of the empty parentheses at the
> end of the formula? I think then my code could work.
> 
> --
> Erich Neuwirth, Didactic Center for Computer Science
> University of Vienna
> Visit our SunSITE at http://sunsite.univie.ac.at
> Phone: +43-1-4277-39902 Fax: +43-1-4277-9399
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 

I think the preferred way to do this is using substitute although
formulas are a bit tricky in that you need to eval them after the
substitution to make sure that the object has class "formula".


> (foo <- eval(substitute(x ~ g * h, list(x = as.name("weights")
weights ~ g * h
> class(foo)
[1] "formula"

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Master's project to coerce linux nvidia drivers to run generalised linear models

[Douglas Bates]

On 1/23/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> I wonder if it would make more sense to get a relatively
> low level package to run on it so that all packages that
> used that low level package would benefit.  The Matrix
> package and the functions runmean and sum.exact in
> package caTools are some things that come to mind.
> Others may have other ideas along these lines.

Martin and I are delighted to learn that the Matrix package is a
"relatively low-level" package :-)

We were of the opinion that the amount of code and design work that
went into it made it a little more sophisticated than that.

More seriously, the approach to speeding up model fitting that has
been most successful to date is to speed up the BLAS (Basic Linear
Algebra Subroutines), especially the Level-3 BLAS.  The bulk of the
computation in the Matrix package takes place in either Lapack (for
dense matrices) or CHOLMOD (for sparse matrices) code and those are
based on calls to the Levels 1, 2 and 3 BLAS.  The Atlas package and
K. Goto's BLAS are designed to obtain the highest level of performance
possible from the CPU on these routines.  I think the easiest way of
incorporating the power of the GPU into the model fitting process
would be to port the BLAS to the GPU.  I also imagine that someone
somewhere has already started on that.

>
> On 1/23/06, Oliver LYTTELTON <[EMAIL PROTECTED]> wrote:
> >
> >
> > Hi,
> >
> > I am working with a friend on a master's project. Our laboratory does a
> > lot of statistical analysis using the R stats package and we also have a
> > lot of under-utilised nvidia cards sitting in the back of our networked
> > linux machines. Our idea is to coerce the linux nvidia driver to run
> > some of our statistical analysis for us. Our first thought was to
> > specifically code up a version of glm() to run on the nvidia cards...
> >
> > Thinking that this might be of use to the broader community we thought
> > we might ask for feedback before starting?
> >
> > Any ideas...
> >
> > Thanks,
> >
> > Olly
> >
> > __
> > R-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Peculiar timing result

[Douglas Bates]

I have been timing a particular model fit using lmer on several
different computers and came up with a peculiar result - the model fit
is considerably slower on a dual-core Athlon 64 using Goto's
multithreaded BLAS than on a single-core processor.

Here is the timing on a single-core Athlon 64 3000+ running under
today's R-devel with version 0.995-5 of the Matrix package.

> library(Matrix)
> data(star, package = 'mlmRev')
> system.time(fm1 <- lmer(math~gr+sx+eth+cltype+(yrs|id)+(1|tch)+(yrs|sch), 
> star, control = list(nit=0,grad=0,msV=1)))
  0  241720.:  1.16440 0.335239  0.0  1.78732 0.867209 0.382318  0.0
  1  239722.:  1.94952 5.0e-10 0.00933767  1.65999 0.858003
0.341520 0.00908757
  2  239580.:  1.95924 0.0884059 0.00933767  1.65308 0.857487
0.339296 0.00954718
  3  239215.:  2.60877 0.0765848 0.0177699  1.45739 0.843562
0.275100 0.0236849
  4  239204.:  2.62582 0.106670 0.0239698  1.41976 0.841086
0.261033 0.0267073
  5  239176.:  2.63149 0.0787924 0.0367185  1.37952 0.838527
0.245076 0.0301134
  6  239141.:  2.64949 0.108534 0.0594586  1.28846 0.832543
0.208404 0.0375456
  7  239049.:  2.64794 0.0789214 0.121782  1.10436 0.819711
0.126101 0.0524965
  8  239004.:  2.66044 0.117957 0.181505 0.932202 0.798982
0.0718116 0.0628958
  9  238944.:  2.66310 0.0819653 0.334477 0.631735 0.740855
0.258359 0.0806590
 10  238893.:  2.72626 0.0975205 0.653432 0.703912 0.666067
0.109922 0.201809
 11  238892.:  2.74381 0.46 0.666155 0.693544 0.662000 0.104060 0.207591
 12  23.:  2.75052 0.0990238 0.689199 0.694588 0.655781
0.106516 0.216460
 13  238861.:  2.80325 0.126935  1.05912 0.733914 0.556162 0.159296 0.360938
 14  238832.:  2.82656 0.117617  1.59471 0.607916 0.441371
0.0749944 0.976142
 15  238811.:  2.87350 0.136332  1.59046 0.653141 0.353763 0.226061  1.79285
 16  238810.:  2.87663 0.125135  1.58992 0.656808 0.352605 0.220488  1.79282
 17  238806.:  2.89342 0.141551  1.58607 0.676523 0.344212 0.181833  1.79268
 18  238804.:  2.90080 0.125137  1.56624 0.682921 0.261295 0.180598  1.74325
 19  238802.:  2.91950 0.128569  1.56836 0.680436 0.336051 0.159940  1.80400
 20  238801.:  2.92795 0.134762  1.56597 0.685121 0.331695 0.145547  1.80414
 21  238801.:  2.93741 0.125667  1.56139 0.687827 0.332700 0.138854  1.81495
 22  238800.:  2.94588 0.131757  1.55294 0.687909 0.330414 0.137834  1.82826
 23  238799.:  2.96867 0.129716  1.52943 0.688678 0.323171 0.139912  1.84615
 24  238799.:  2.98994 0.133378  1.52188 0.700038 0.337387 0.131403  1.77623
 25  238799.:  3.00312 0.135308  1.51475 0.697550 0.311750 0.145683  1.78037
 26  238799.:  3.00461 0.129920  1.51083 0.697666 0.306722 0.138745  1.81188
 27  238799.:  3.00504 0.134882  1.50539 0.696745 0.302949 0.135897  1.84405
 28  238799.:  3.00422 0.134013  1.47947 0.698115 0.303243 0.133806  1.86486
 29  238799.:  3.00819 0.134378  1.48185 0.701871 0.307097 0.134637  1.84996
 30  238799.:  3.01313 0.134279  1.49098 0.702883 0.304788 0.133682  1.86254
 31  238799.:  3.01291 0.134253  1.49060 0.701818 0.303155 0.133771  1.84613
 32  238799.:  3.01142 0.134314  1.48921 0.701782 0.303589 0.134439  1.84653
 33  238799.:  3.01174 0.134315  1.48926 0.701641 0.304120 0.134145  1.84635
 34  238799.:  3.01175 0.134304  1.48942 0.701740 0.303762 0.134185  1.84649
 35  238799.:  3.01173 0.134307  1.48937 0.701724 0.303809 0.134206  1.84647
[1] 43.10  3.78 48.41  0.00  0.00


(If you run the timing yourself and don't want to see the iteration
output, take the msV=1 out of the control list.  I keep it in there so
I can monitor the progress.)

If I time the same model fit on a dual-core Athlon 64 X2 3800+ with
the same version of R, BLAS and Matrix package, the timing ends up
with something like

90 140 235 0 0

I do see that the multithreading is working on a calculation that is
essentially BLAS-bound such as

> mm <- Matrix(rnorm(1e6), nc = 1e3)
> system.time(crossprod(mm))
[1] 0.57 0.02 0.60 0.00 0.00

On the X2 processor it still takes about 0.6 seconds user time but
only 0.3 seconds elapsed time when the machine is otherwise idle and
both cores are available for the calculation.

Any suggestions why the dual-core processor is so much slower than the
single core processor?

By the way, I would be interested in accumulating timings of this
model fit on other systems.  If you do time it please send me
(off-list) a summary of the version of R, version of the accelerated
BLAS if you use them, processor speed and configuration (i.e.
multiprocessor, multicore, etc.) and, if you know it, memory speed.

This is an example of a complex multilevel model with crossed grouping
factors fit to a relatively large (3 observations on 1
students, 1400 teachers and 80 schools) data set.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Peculiar timing result

[Douglas Bates]

I don't think this calculation is memory-bound at all and I would be
surprised if changing to a 32-bit environment would change things.  I
do have a 32-bit chroot environment on these machines (needed for
things like wine and acroread) so I'll try that out but I think I will
need to use Atlas as the accelerated BLAS rather than Goto's BLAS.  I
imagine the Opteron/Athlon 64 version of Goto's BLAS assumes a 64-bit
environment.

On 3/3/06, Paul Gilbert <[EMAIL PROTECTED]> wrote:
> Doug
>
> This is probably not your reason, but I am finding my dual core Athlon
> 64 is much slower running 64 bit Linux and R than it was running 32 bit
> Linux and R. All the programs are bigger.  (Some, like the clock applet,
> are a lot bigger for no obvious reason.)  The difference is enough to
> put my meager 1GB machine into swapping much more, with the result that
> it is a lot slower.
>
> Paul
>
> Douglas Bates wrote:
>
> >I have been timing a particular model fit using lmer on several
> >different computers and came up with a peculiar result - the model fit
> >is considerably slower on a dual-core Athlon 64 using Goto's
> >multithreaded BLAS than on a single-core processor.
> >
> >Here is the timing on a single-core Athlon 64 3000+ running under
> >today's R-devel with version 0.995-5 of the Matrix package.
> >
> >
> >
> >>library(Matrix)
> >>data(star, package = 'mlmRev')
> >>system.time(fm1 <- lmer(math~gr+sx+eth+cltype+(yrs|id)+(1|tch)+(yrs|sch), 
> >>star, control = list(nit=0,grad=0,msV=1)))
> >>
> >>
> >  0  241720.:  1.16440 0.335239  0.0  1.78732 0.867209 0.382318  
> > 0.0
> >  1  239722.:  1.94952 5.0e-10 0.00933767  1.65999 0.858003
> >0.341520 0.00908757
> >  2  239580.:  1.95924 0.0884059 0.00933767  1.65308 0.857487
> >0.339296 0.00954718
> >  3  239215.:  2.60877 0.0765848 0.0177699  1.45739 0.843562
> >0.275100 0.0236849
> >  4  239204.:  2.62582 0.106670 0.0239698  1.41976 0.841086
> >0.261033 0.0267073
> >  5  239176.:  2.63149 0.0787924 0.0367185  1.37952 0.838527
> >0.245076 0.0301134
> >  6  239141.:  2.64949 0.108534 0.0594586  1.28846 0.832543
> >0.208404 0.0375456
> >  7  239049.:  2.64794 0.0789214 0.121782  1.10436 0.819711
> >0.126101 0.0524965
> >  8  239004.:  2.66044 0.117957 0.181505 0.932202 0.798982
> >0.0718116 0.0628958
> >  9  238944.:  2.66310 0.0819653 0.334477 0.631735 0.740855
> >0.258359 0.0806590
> > 10  238893.:  2.72626 0.0975205 0.653432 0.703912 0.666067
> >0.109922 0.201809
> > 11  238892.:  2.74381 0.46 0.666155 0.693544 0.662000 0.104060 
> > 0.207591
> > 12  23.:  2.75052 0.0990238 0.689199 0.694588 0.655781
> >0.106516 0.216460
> > 13  238861.:  2.80325 0.126935  1.05912 0.733914 0.556162 0.159296 
> > 0.360938
> > 14  238832.:  2.82656 0.117617  1.59471 0.607916 0.441371
> >0.0749944 0.976142
> > 15  238811.:  2.87350 0.136332  1.59046 0.653141 0.353763 0.226061  
> > 1.79285
> > 16  238810.:  2.87663 0.125135  1.58992 0.656808 0.352605 0.220488  
> > 1.79282
> > 17  238806.:  2.89342 0.141551  1.58607 0.676523 0.344212 0.181833  
> > 1.79268
> > 18  238804.:  2.90080 0.125137  1.56624 0.682921 0.261295 0.180598  
> > 1.74325
> > 19  238802.:  2.91950 0.128569  1.56836 0.680436 0.336051 0.159940  
> > 1.80400
> > 20  238801.:  2.92795 0.134762  1.56597 0.685121 0.331695 0.145547  
> > 1.80414
> > 21  238801.:  2.93741 0.125667  1.56139 0.687827 0.332700 0.138854  
> > 1.81495
> > 22  238800.:  2.94588 0.131757  1.55294 0.687909 0.330414 0.137834  
> > 1.82826
> > 23  238799.:  2.96867 0.129716  1.52943 0.688678 0.323171 0.139912  
> > 1.84615
> > 24  238799.:  2.98994 0.133378  1.52188 0.700038 0.337387 0.131403  
> > 1.77623
> > 25  238799.:  3.00312 0.135308  1.51475 0.697550 0.311750 0.145683  
> > 1.78037
> > 26  238799.:  3.00461 0.129920  1.51083 0.697666 0.306722 0.138745  
> > 1.81188
> > 27  238799.:  3.00504 0.134882  1.50539 0.696745 0.302949 0.135897  
> > 1.84405
> > 28  238799.:  3.00422 0.134013  1.47947 0.698115 0.303243 0.133806  
> > 1.86486
> > 29  238799.:  3.00819 0.134378  1.48185 0.701871 0.307097 0.134637  
> > 1.84996
> > 30  238799.:  3.01313 0.134279  1.49098 0.702883 0.304788 0.133682  
> > 1.86254
> > 31  238799.:  3.01291 0.134253  1.49060 0.701818 0.303155 0.133771  
> > 1.84613
> > 32  238799.:  3.01142 0.1343

Re: [Rd] Peculiar timing result

[Douglas Bates]

On 3/11/06, Prof Brian Ripley <[EMAIL PROTECTED]> wrote:
> Here is a summary of some results on a dual Opteron 252 running FC3
>
> 64-bit gcc 3.4.5
> R's blas34.83  3.45 38.56
> ATLAS   36.70  3.28 40.14
> ATLAS multithread   76.85  5.39 82.29
> Goto 1 thread   36.17  3.44 39.76
> Goto multithread   178.06 345.97 467.99
> ACML49.69  3.36 53.23
>
> 64-bit gcc 4.1.0
> R's blas34.98  3.49 38.55
> 32-bit gcc 3.4.5
> R's blas33.72  3.27 36.99
> 32-bit gcc 4.1.0
> R's blas34.62  3.25 37.93
>
> The timings are not that repeatable, but the message seems clear that
> this problem does not benefit from a tuned BLAS and the overhead from
> multiple threads is harmful.  (The gcc 4.1.0 results took fewer
> iterations, which skews the results in its favour.)
>
> And my 2GHz Pentium M laptop under Windows gave 39.96  3.68 44.06.
>
> Clearly the Goto BLAS has a problem here: the results are slower on a dual
> 252 than a dual 248 (see below).

Thanks for the information on the timings.  It happens that this
message ended up in my R-help folder and I only got around to reading
that folder today.

Is it ok with you if I forward this message to Simon Urbanek?  I am
having similar difficulties in the timing with R on a dual-core Intel
MacBook.
>
>
> On Fri, 3 Mar 2006, Prof Brian Ripley wrote:
>
> > On Fri, 3 Mar 2006, Douglas Bates wrote:
> >
> >> I have been timing a particular model fit using lmer on several
> >> different computers and came up with a peculiar result - the model fit
> >> is considerably slower on a dual-core Athlon 64 using Goto's
> >> multithreaded BLAS than on a single-core processor.
> >
> > Is there a Goto BLAS tuned for that chip?  I can only see one tuned for an
> > (unspecified) Opteron.  L1 and L2 cache sizes do sometimes matter a lot
> > for tuned BLAS, and (according to the AMD site I just looked up) the X2
> > 3800+ only has a 512Kb per core L2 cache.  Opterons have a 1Mb L2 cache.
> >
> > Also, the very large system time seen in the dual-core run is typical of
> > what I see when pthreads is not working right, and I suggest you try a
> > limit of one thread (see the R-admin manual).  On our dual-processor
> > Opteron 248 that ran in 44 secs instead of 328.
> >
> >> Here is the timing on a single-core Athlon 64 3000+ running under
> >> today's R-devel with version 0.995-5 of the Matrix package.
> >>
> >>> library(Matrix)
> >>> data(star, package = 'mlmRev')
> >>> system.time(fm1 <- lmer(math~gr+sx+eth+cltype+(yrs|id)+(1|tch)+(yrs|sch), 
> >>> star,
> > control = list(nit=0,grad=0,msV=1)))
> >> [1] 43.10  3.78 48.41  0.00  0.00
> >>
> >>
> >> (If you run the timing yourself and don't want to see the iteration
> >> output, take the msV=1 out of the control list.  I keep it in there so
> >> I can monitor the progress.)
> >>
> >> If I time the same model fit on a dual-core Athlon 64 X2 3800+ with
> >> the same version of R, BLAS and Matrix package, the timing ends up
> >> with something like
> >>
> >> 90 140 235 0 0
> > 
> >
> >
>
> --
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] gsummary function (nlme library) (PR#8782)

[Douglas Bates]

The documentation for gsummary describes the argument FUN as

 FUN: an optional summary function or a list of summary functions
  to be applied to each variable in the frame.  The function or
  functions are applied only to variables in 'object' that vary
  within the groups defined by 'groups'.  Invariant variables
  are always summarized by group using the unique value that
  they assume within that group.  If 'FUN' is a single function
  it will be applied to each non-invariant variable by group to
  produce the summary for that variable.  If 'FUN' is a list of
  functions, the names in the list should designate classes of
  variables in the frame such as 'ordered', 'factor', or
  'numeric'.  The indicated function will be applied to any
  non-invariant variables of that class.  The default functions
  to be used are 'mean' for numeric factors, and 'Mode' for
  both 'factor' and 'ordered'.  The 'Mode' function, defined
  internally in 'gsummary', returns the modal or most popular
  value of the variable.  It is different from the 'mode'
  function that returns the S-language mode of the variable.

so the behavior you noticed is documented.

The "summary" in "gsummary" is used in the sense of a representative
value, not in the more general sense of a numerical summary of any
sort.  If the values do not vary within a group then the common value
within the group is, according to our definition, the representative
value.


On 4/19/06, [EMAIL PROTECTED] <[EMAIL PROTECTED]> wrote:
> Full_Name: Ben Saville
> Version: 2.1
> OS: Windows XP
> Submission from: (NULL) (152.2.94.145)
>
>
> I'm using the gsummary function to calculate a sum of V1 (column one) from my
> data 'mytest' by group (V2,or column 2).  If V1 (the variable of interest) is
> all the same value (in this case all 2's), I do not get back the correct
> summation.  If there is at least one difference in V1 (all 2's except for one
> 1), it gives me correct values.  So either I am doing something wrong or there
> is a bug in the gsummary function.
>
># Incorrect sums
> mytest <- as.data.frame(matrix(c(2,rep(2,8),1,1,2,2,2,3,3,3,3),ncol=2))
> mytest
> gsummary(mytest,form=V1~1|V2, FUN=sum)[,1]
>
># Correct sums
> mytest <- as.data.frame(matrix(c(1,rep(2,8),1,1,2,2,2,3,3,3,3),ncol=2))
> mytest
> gsummary(mytest,form=V1~1|V2, FUN=sum)[,1]
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] R Project accepted for the Google Summer of Code

[Douglas Bates]

We are very pleased that an R project has been selected as one of the GNU
projects for the Google Summer of Code 2006 (http://code.google.com/soc/).

Miguel Angel R. Manese, an M.S. Statistics student at the University of
the Philippines, will be working with Douglas Bates and Brian Ripley on a
project entitled `SQLite Data Frames for R' to let R store very large
datasets in an SQLite database transparently, implementing primitive
operations for data stored in SQLite so that they behave exactly like
ordinary data frames to the users.  It is likely that the project will
result in the first instance in an R package (which may in due course
become part of the tarball).

Congratulations, Miguel!

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel