Re: [Rd] R-devel Digest, Vol 112, Issue 8
I'll not be able to comment on the use of C like this, but will warn that I wrote the routines that became Nelder-Mead, CG, and BFGS in optim() in the mid 1970s. CG never did as well as I would like, but the other two routines turned out pretty well. However, in nearly 40 years, there are a few improvements, particularly in handling bounds and masks (fixed parameters). For all-R routines see Rcgmin and Rvmmin. Rather than directly use the optim() routines in C, you may want to use some more modern ones, but the choice may be dependent on your problem. JN On 06/09/2012 06:00 AM, r-devel-requ...@r-project.org wrote: > Message: 2 > Date: Fri, 8 Jun 2012 09:40:17 -0400 > From: Edward Worbis > To: r-devel@R-project.org > Subject: [Rd] Working with optim in C > Message-ID: > > Content-Type: text/plain > > I've searched to find examples of how to work with the C versions of > optim. > > I've separated out the function just to test on it alone, and currently I'm > attempting to use fmmin as follows: > > > !~~CODE ~~! __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] How to set up an object to share data across diverse functions
In the past 6 months I've been struggling with an issue that has been raised periodically on the lists. This is the need to share information across a group of functions, possibly from different packages. So far I've found solutions that are either quite clumsy or else don't work as have (likely incorrectly) read the manuals. I'll be writing up my experiences as a way to help others, but am wondering if there is an accepted way to establish an object that can be read and written by a number of functions. I'm currently using list2env to create a list under the "umbrella" function, and if necessary pass this to other functions by name through the dot arguments, though preferably I'd rather use a pre-arranged and fixed name (this does not seem to work unless the object is in the global environment, and then I've occasionally had locked bindings errors, which I'll freely admit I don't understand at all). I also have to be very careful with dot arguments, especially if they are matrices. I've a couple of examples, but they are pretty long, so I'll save for off-list transmission. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Converting expression to a function
Previously, I've posted queries about this, and thanks to postings and messages
in
response have recently had some success, to the extent that there is now a
package called
nlmrt on the R-forge project https://r-forge.r-project.org/R/?group_id=395 for
solving
nonlinear least squares problems that include small or zero residual problems
via a
Marquardt method using a call that mirrors the nls() function. nls()
specifically warns
against zero residual problems.
However, I would still like to be able to convert expressions with example
vectors of
parameters to functions that optim() and related functions can use. The code
below gets
"almost" there, but
1) Can the code be improved / cleaned up?
2) Can the eval() of the output of the Form2resfun be avoided?
3) Can the extraction of the parameter names be embedded in the function rather
than put
separately?
Off-list responses are likely best at this stage, while the tedious details are
sorted
out. I will post a summary in a couple of weeks of the results. Collaborations
re: this
and the larger package welcome, as there is considerable testing and tuning to
do, but
preliminary experience is encouraging.
John Nash
# - code block ---
rm(list=ls()) # clear workspace
Form2resfun <- function(f, p ) {
cat("In Form2resfun\n")
xx <- all.vars(f)
fp <- match(names(p), xx) # Problem in matching the names of params
xx2 <- c(xx[fp], xx[-fp])
ff <- vector("list", length(xx2))
names(ff) <- xx2
sf<-as.character(f)
if ((length(sf)!=3) && (sf[1]!="~")) stop("Bad model formula
expression")
lhs<-sf[2] # NOTE ORDER formula with ~ puts ~, lhs, rhs
rhs<-sf[3]
# And build the residual at the parameters
resexp<-paste(rhs,"-",lhs, collapse=" ")
fnexp<-paste("crossprod(",resexp,")", sep="")
ff[[length(ff) + 1]] <- parse(text=fnexp)
# want crossprod(resexp)
myfn<-as.function(ff, parent.frame())
}
# a test
y<-c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443,
38.558, 50.156, 62.948, 75.995, 91.972) # for testing
t<-1:length(y) # for testing
f<- y ~ b1/(1+b2*exp(-1*b3*t))
p<-c(b1=1, b2=1, b3=1)
b<-p
npar<-length(b)
for (i in 1:npar){
bbit<-paste(names(b)[[i]],"<-",b[[i]])
eval(parse(text=bbit))
}
tfn<-Form2resfun(f, b)
ans<-eval(tfn(t=t,y=y, b))
print(ans)
# - end code block ---
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] bug in sum() on integer vector
I agree that where the overflow occurs is not critical (one can go back to
cumsum and find
out). I am assuming that Uwe still wants to know there has been an overflow at
some point
i.e., a warning. This could become more "interesting" as parallel computation
causes
different summation orderings on sums of large numbers of items.
JN
On 12/14/2011 03:58 PM, Uwe Ligges wrote:
>
>
> On 14.12.2011 17:19, peter dalgaard wrote:
>>
>> On Dec 14, 2011, at 16:19 , John C Nash wrote:
>>
>>>
>>> Following this thread, I wondered why nobody tried cumsum to see where the
>>> integer
>>> overflow occurs. On the shorter xx vector in the little script below I get
>>> a message:
>>>
>>> Warning message:
>>> Integer overflow in 'cumsum'; use 'cumsum(as.numeric(.))'
>>>>
>>>
>>> But sum() does not give such a warning, which I believe is the point of
>>> contention. Since
>>> cumsum() does manage to give such a warning, and show where the overflow
>>> occurs, should
>>> sum() not be able to do so? For the record, I don't class the non-zero
>>> answer as an error
>>> in itself. I regard the failure to warn as the issue.
>>
>> It (sum) does warn if you take the two "halves" separately. The issue is
>> that the
>> overflow is detected at the end of the summation, when the result is to be
>> saved to an
>> integer (which of course happens for all intermediate sums in cumsum)
>>
>>> x<- c(rep(180003L, 1000), -rep(120002L, 1500))
>>> sum(x[1:1000])
>> [1] NA
>> Warning message:
>> In sum(x[1:1e+07]) : Integer overflow - use sum(as.numeric(.))
>>> sum(x[1001:2500])
>> [1] NA
>> Warning message:
>> In sum(x[1001:1.5e+07]) : Integer overflow - use sum(as.numeric(.))
>>> sum(x)
>> [1] 4996000
>>
>> There's a pretty easy fix, essentially to move
>>
>> if(s> INT_MAX || s< R_INT_MIN){
>> warningcall(call, _("Integer overflow - use sum(as.numeric(.))"));
>> *value = NA_INTEGER;
>> }
>>
>> inside the summation loop. Obviously, there's a speed penalty from two FP
>> comparisons
>> per element, but I wouldn't know whether it matters in practice for anyone.
>>
>
>
> I don't think I am interested in where the overflow happens if I call sum()...
>
> Uwe
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] bug in sum() on integer vector
Following this thread, I wondered why nobody tried cumsum to see where the
integer
overflow occurs. On the shorter xx vector in the little script below I get a
message:
Warning message:
Integer overflow in 'cumsum'; use 'cumsum(as.numeric(.))'
>
But sum() does not give such a warning, which I believe is the point of
contention. Since
cumsum() does manage to give such a warning, and show where the overflow
occurs, should
sum() not be able to do so? For the record, I don't class the non-zero answer
as an error
in itself. I regard the failure to warn as the issue.
For info, on my Ubnuntu Lucid 10.04 system that has 4 GB of RAM but no swap,
the last line
of the script to do the int64 sum chugs for about 2 minutes then gives "Killed"
and
returns to the terminal prompt. It also seems to render some other applications
unstable
(I had Thunderbird running to read R-devel, and this started to behave
strangely after the
crash, and I had to reboot.) I'm copying Romain as package maintainer, and I'll
be happy
to try to work off-list to figure out how to avoid the "Killed" result. (On a
16GB
machine, I got the 0 answer.)
Best,
John Nash
Here's the system info and small script.
>> sessionInfo()
> R version 2.14.0 (2011-10-31)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C
> [3] LC_TIME=en_US.utf8LC_COLLATE=en_US.utf8
> [5] LC_MONETARY=en_US.utf8LC_MESSAGES=en_US.utf8
> [7] LC_PAPER=CLC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] int64_1.1.2
>>
## sumerr.R 20111214
library(int64)
x <- c(rep(180003L, 1000), -rep(120002L, 1500))
xx <- c(rep(180003L, 1000), -rep(120002L, 1500))
sum(x)
sum(as.double(x))
sum(xx)
sum(as.double(xx))
cumsum(xx)
cumsum(as.int64(xx))
tmp<-readline("Now try the VERY SLOW int64")
sum(as.int64(x))
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Detecting typo in function argument
With some chagrin after spending a couple of hours trying to debug a script, I
realized I
had typed in something like
ans<-optimx(start, myfn, mygr, lower<-lo, upper=up)
that is, the "<-" rather than "=". The outcome on my machine was a non-obvious
error
several layers deep in the call stack. For info, optim() seems to stop much
more quickly.
The error is "obvious", but I'm wondering if there is a simple way to trap or
warn of it.
For completeness, I include the commands I used to force the error. Note that
it will only
work fully with the latest (R-forge) version of optimx/optfntools because of
the form of
gr="gfwd" that allows a choice of different numerical gradient routines.
This is a curiosity rather than a necessity, but if there is a simple way to
check, I'll
put it in my codes.
Cheers,
JN
rm(list=ls())
start<-rep(3,6)
lo<-rep(2,6)
up<-rep(4,6)
flb <- function(x)
{ p <- length(x); sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) }
ans<-optim(start, flb, lower=lo, upper=up)
ans
ans<-optim(start, flb, lower<-lo, upper=up)
ans
ans1<-optim(start, flb, lower<-lo, upper=up)
ans1
require(optimx)
ans1x<-optimx(start, flb, lower<-lo, upper=up)
ans1x<-optimx(start, flb, gr="gfwd",lower<-lo, upper=up)
ans1<-optim(start, flb, gr=NULL,lower<-lo, upper=up)
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] testing on platforms, was Fortran error
I think the poster is interested in being able to try the build/check on a Mac in the fashion that Winbuilder does. That is, rather than have CRAN do all the platform checks, is there a way to submit a package to be tested for just one platform? It may be useful to have such a facility so package builders can test before submission, and from my own experience I believe would help to render packages more platform neutral. I suspect that all the pieces are in place to do this, but perhaps a nice front end is needed. Is that something for a Google Summer of Code task? Best, JN On 11/21/2011 06:00 AM, r-devel-requ...@r-project.org wrote: > What does 'submit a test mean'? Any test you write can check the > value of Sys.info()["sysname"]: it is "Darwin" on Mac OS X. > Or R.version$platform, which will be something like > i386-apple-darwin9.8.0. __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] overly long lines in pdf output of manual from package
Thanks for quick reply. To help others I'll put in short comments in edited thread. On 11/02/2011 03:54 PM, Prof Brian Ripley wrote: > On Wed, 2 Nov 2011, John C Nash wrote: > >> In re-factoring my optimx package, I'm finding that the pdf output has some >> lines that are >> outside the margins (and there are warnings in R CMD check optimx). >> 1) How do I adjust the DESCRIPTION file to avoid too long a line. I can add >> \cr elements to the Rd... > > Wrap lines? And is this R 2.14.0, as the issues are easier with the fonts > used there? R 2.14.0 -- installed from M Rutter repository today. I don't know why it wraps sometimes and not others, but have occasionally seen similar over-long lines with other TEX files with the TeXlive I'm using, even with the texlive-fonts-extra installed. > >> 2) Is there an option control that does the job? I've not found one in >> searching, but >> could easily have chosen poor search terms. > > R has ways to reformat description files. See ?write.dcf. Thanks. I would not have found that! Hence example lines would be Depends: pkg1, pkg2, pkg3, pkg4, pkg5, pkg6. i.e., spaces on second line. >> 3) To avoid running the whole of R CMD check ... > > R CMD R2pdf on a package source directory is what R CMD check uses. > When in directory just above that with optimx package in it: R CMD Rd2pdf optimx works quickly and well. Note Rd2pdf not R2pdf (I only noticed just now). I had misunderstood the Rd2pdf manual about directories. My bad. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] overly long lines in pdf output of manual from package
In re-factoring my optimx package, I'm finding that the pdf output has some lines that are outside the margins (and there are warnings in R CMD check optimx). Clearly I can fix this by editing the tex file that is generated, but the build infrastructure would still get things wrong in automatic processing. So that gives rise to 3 questions: 1) How do I adjust the DESCRIPTION file to avoid too long a line (optimx uses a lot of other packages). I can add \cr elements to the Rd files to control those lines, but DESCRIPTION doesn't allow this. 2) Is there an option control that does the job? I've not found one in searching, but could easily have chosen poor search terms. 3) To avoid running the whole of R CMD check to rebuild the manual, do I need to dig into R CMD check to find the bits that just build the manual, or is there an option to Rd2pdf to do this? I've found how to process a list of files with Rd2pdf, but need to add the DESCRIPTION stuff. Best John Nash __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] list of operations that preserve attributes
I've been trying to prepare some wrappers for optimization objective functions and gradients that use try() to check for computable functions. I'm also trying to do scaling as per optim() so that a number of methods that are available on CRAN and R-Forge can use parameter and function scaling. This got me into minor trouble when I scaled a gradient or Hessian and the "inadmissible" attribute I had created suddenly disappeared. When I discovered why -- that some operations don't preserve attributes -- I started to look for a list of which ops preserve attributes and which don't. The code snippet below shows that matrix multiplication (%*%) does not, while multiplication by a vector (simple *) does. I'm not really surprised that some ops don't preserve attributes, particularly those that are binary like multiplication, but it might be nice to know which operations and special functions do so. rm(list=ls()) m<-matrix(1:4,nrow=2, ncol=2) print(m) attributes(m) attr(m,"check")<-"***" attributes(m) bigm<-10*m str(bigm) bigm1<-diag(c(1,1))%*%m str(bigm1) bigm1<-c(1,2)*m str(bigm1) print(bigm1) arraym<-as.array(m) str(arraym) tanm<-tan(m) str(tanm) Best, John Nash __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] control list gotcha
This is mainly a reminder to others developing R packages to be careful not to supply control list items that are not used by the called package. Optimx is a wrapper package that aims to provide a common syntax to a number of existing optimization packages. Recently in extending optimx package I inadvertently introduced a new control for optimx which is NOT in any of the wrapped optimization packages. There are probably other methods of keeping things tidy, but I copy the control list and null out unwanted elements for each of the called packages. I missed this in a couple of places in the R-forge development version of optimx (I'm working on fixing these, but they are still there at the moment). The "nasty" here was that the package mostly works, with plausible but not very good results for some of the optimizers. If it crashed and burned, it would have been noticed sooner. There is also a potential interaction with a use of the dot-dot-dot variable to pass scaling information. If there are ideas on how to quickly reveal errors related to calling sequences involving control lists and "...", I'd welcome them (off-list?), and be prepared to summarize findings in a vignette. Best, JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] An example of very slow computation
This message is about a curious difference in timing between two ways of
computing the
same function. One uses expm, so is expected to be a bit slower, but "a bit"
turned out to
be a factor of >1000. The code is below. We would be grateful if anyone can
point out any
egregious bad practice in our code, or enlighten us on why one approach is so
much slower
than the other. The problem arose in an activity to develop guidelines for
nonlinear
modeling in ecology (at NCEAS, Santa Barbara, with worldwide participants), and
we will be
trying to include suggestions of how to prepare problems like this for
efficient and
effective solution. The code for nlogL was the "original" from the worker who
supplied the
problem.
Best,
John Nash
--
cat("mineral-timing.R == benchmark MIN functions in R\n")
# J C Nash July 31, 2011
require("microbenchmark")
require("numDeriv")
library(Matrix)
library(optimx)
# dat<-read.table('min.dat', skip=3, header=FALSE)
# t<-dat[,1]
t <- c(0.77, 1.69, 2.69, 3.67, 4.69, 5.71, 7.94, 9.67, 11.77, 17.77,
23.77, 32.77, 40.73, 47.75, 54.90, 62.81, 72.88, 98.77, 125.92, 160.19,
191.15, 223.78, 287.70, 340.01, 340.95, 342.01)
# y<-dat[,2] # ?? tidy up
y<- c(1.396, 3.784, 5.948, 7.717, 9.077, 10.100, 11.263, 11.856, 12.251,
12.699,
12.869, 13.048, 13.222, 13.347, 13.507, 13.628, 13.804, 14.087, 14.185, 14.351,
14.458, 14.756, 15.262, 15.703, 15.703, 15.703)
ones<-rep(1,length(t))
theta<-c(-2,-2,-2,-2)
nlogL<-function(theta){
k<-exp(theta[1:3])
sigma<-exp(theta[4])
A<-rbind(
c(-k[1], k[2]),
c( k[1], -(k[2]+k[3]))
)
x0<-c(0,100)
sol<-function(t)100-sum(expm(A*t)%*%x0)
pred<-sapply(dat[,1],sol)
-sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE))
}
getpred<-function(theta, t){
k<-exp(theta[1:3])
sigma<-exp(theta[4])
A<-rbind(
c(-k[1], k[2]),
c( k[1], -(k[2]+k[3]))
)
x0<-c(0,100)
sol<-function(tt)100-sum(expm(A*tt)%*%x0)
pred<-sapply(t,sol)
}
Mpred <- function(theta) {
# WARNING: assumes t global
kvec<-exp(theta[1:3])
k1<-kvec[1]
k2<-kvec[2]
k3<-kvec[3]
# MIN problem terbuthylazene disappearance
z<-k1+k2+k3
y<-z*z-4*k1*k3
l1<-0.5*(-z+sqrt(y))
l2<-0.5*(-z-sqrt(y))
val<-100*(1-((k1+k2+l2)*exp(l2*t)-(k1+k2+l1)*exp(l1*t))/(l2-l1))
} # val should be a vector if t is a vector
negll <- function(theta){
# non expm version JN 110731
pred<-Mpred(theta)
sigma<-exp(theta[4])
-sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE))
}
theta<-rep(-2,4)
fand<-nlogL(theta)
fsim<-negll(theta)
cat("Check fn vals: expm =",fand," simple=",fsim," diff=",fand-fsim,"\n")
cat("time the function in expm form\n")
tnlogL<-microbenchmark(nlogL(theta), times=100L)
tnlogL
cat("time the function in simpler form\n")
tnegll<-microbenchmark(negll(theta), times=100L)
tnegll
ftimes<-data.frame(texpm=tnlogL$time, tsimp=tnegll$time)
# ftimes
boxplot(log(ftimes))
title("Log times in nanoseconds for matrix exponential and simple MIN fn")
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] R-devel Digest, Vol 101, Issue 2
I've been finding that the "loose ends" in many of these older codes cause more trouble than it is worth in their use. When I encounter them, I've attempted to re-program the algorithm in R. A lot of the Fortran code is because of the software structure the author used and nothing to do with the job to be done. If you can prepare an R function for this, you'd be doing the R community a favour. You may also find that a judicious combination of optimize() and grid search gets the task done satisfactorily. Best, JN > Date: Sat, 2 Jul 2011 01:43:03 +0530 > From: Mohit Dayal > To: R-devel@r-project.org > Subject: [Rd] AS Algorithms > > I would like to use one of the AS Algorithms that used to be published in > the journal Applied Statistics of the Royal Statistical Society (Series C). > FORTRAN code based on these are available on the Statlib website at > > http://lib.stat.cmu.edu/modules.php?op=modload&name=PostWrap&file=index&page=apstat/ [snip] > > *BTW I am looking at AS 133 : Finding the global maximum or minimum of a > function of 1 variable. > * > __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Default values in control list of function
In building a function for a package, I'd like to set the defaults in a control
list,
e.g.,
makeg<-function(parameters, eps = 1.0e-7, control=list(showwork=TRUE,
rubbish=1.0e+101)){
etc.
}
This does not provide showwork or rubbish within the function if control() is
not fully
specified. Copying others, I've previously set the control defaults within the
function
then matched names, and that certainly works. I'm happy to keep using that
approach, but
would like to have a pointer to good practice for doing this (I'm prepared to
write or
help write a short tutorial if none exists).
I can understand that trying to do the default replacement at a second or
greater level
could give problems. A search through the language definition and some googling
didn't
show up any obvious warning about this, but I'm sure there are some comments
somewhere.
Can anyone suggest some URLs?
One gets used to the nicety of being able to specify defaults and forgets that
even R has
some limits.
JN
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] return(); was Suggestions for "good teaching" packages
I tend to code with return(), at least in development, because I've once stepped in the cowpad of ans<- list() then attr(ans ) and forgot to do another ans so got only part of what I wanted. Perhaps its just my thinking style, but I agree with some others who suggest that it's not such a bad idea to be explicit about what one is doing. I prefer pedestrian code that I can understand easily and quickly fix/modify rather than highly optimized and uncommented brilliance that I cannot reuse. Given the overhead of return(), I'll likely switch to ans # return(ans) to make my programs clear, especially to non-R folk migrating in. I have also been writing optimization functions. Modularizing might be a nice student exercise, as well as avoiding early return()s, but Canada isn't wide enough for all the indents of the else clauses when methods crash at different stages and we want to return a very simple structure with partial data etc. Reminds me of the great "GOTO" debate some 30+ years ago. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Google Summer of Code 2011 - credit where it is due
In my reminder that GSoC project proposals are requested (to R wiki developers' projects page for GSoC 2011), I mentioned that Dirk Eddelbuettel had acted as leader for the R Foundation activity on GSoC prior to handing the torch to Claudia Beleites and I for this year. I should have mentioned that for 2009 Manuel Eugster wore the hat, and in 2008 Friedrich Leisch. GSoC has been running since 2005, but we were not involved for the first three years. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Google Summer of Code 2011
The 2011 Google Summer of Code will soon be open for organizations to submit potential projects for which students may apply (with detailed plans) for funding. We have some proposals in process at http://rwiki.sciviews.org/doku.php?id=developers:projects:gsoc2011 Note that projects do need to have mentors. Google has so far had a 1 mentor per project policy, but informally there have been support teams. Claudia Beleites and I are also administering the gsoc-r discussion on google groups mentioned at the bottom of the page linked above. Thanks to Dirk E. who has run this up to now and passed the torch to us. John Nash __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Competing with one's own work
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. Best, John Nash __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] [R] DBLEPR?
Thanks for comments from several folk, fix from Doug Bates and start to finding new email for ucminf maintainer. Summary of responses: DBLEPR and its relations are briefly mentioned but with no or minimal examples in Writing R Extensions as well as Venables and Ripley book. I should have emphasized I was especially seeking examples, and Brian Ripley has pointed out that cluster, mda and tweedie among others use these routines. Fortran WRITE or PRINT routines, even if unused by code, may interfere with C output, especially in Windows. (See R-devel post from Brian Ripley). Doug Bates provided a different approach and code that appears to fix ucminf. We have a lead on contacting the maintainer, who has moved, so hopefully the package can be updated. My conclusion from this: xxxPR routines are still usable, but likely not the best approach for output from Fortran routines. Once again, thanks to all. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] DBLEPR?
My reaction is leaning heavily towards "Virtuoso!" as opposed to "Show Off!". Thanks very much. JN On 11/16/2010 05:39 PM, Douglas Bates wrote: > Try this. > > On Tue, Nov 16, 2010 at 4:06 PM, Prof. John C Nash wrote: >> We've tried to contact Stig since July. Possibly he changed emails. >> >> My thought was to use Rprintf as suggested and was looking into doing that >> to see if our >> optimx problems would go away. Will keep it as open issue while we give a >> bit more time >> for response, and welcome input from others. >> >> JN >> >> >> On 11/16/2010 04:52 PM, Douglas Bates wrote: >>> 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] DBLEPR?
We've tried to contact Stig since July. Possibly he changed emails. My thought was to use Rprintf as suggested and was looking into doing that to see if our optimx problems would go away. Will keep it as open issue while we give a bit more time for response, and welcome input from others. JN On 11/16/2010 04:52 PM, Douglas Bates wrote: > 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] DBLEPR?
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. 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. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] BLAS benchmarks on R 2.12.0 and related performance issues
This thread pointed out that the "plain vanilla" library for linear algebra outperformed the fancy ones for the original poster -- and he had mentioned this, but still got "you ought to " advice that was inappropriate and ignored his stated experience. I've been surprised sometimes myself with performance results, and I'm now getting wary of making pronouncements. I always thought statistics was about using data sensibly to uncover truths. This is a case where we can actually get the data pretty easily. Perhaps the R posting guide needs an addition to say "When discussing performance, you really should present some data and not just speculation." Ultimately we need good performance benchmarks. They are difficult to set up properly and tedious to run. Maybe a good subject for a Google Summer of Code project for next year or some undergraduate projects. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] iteration count in optim()
Date: Fri, 6 Aug 2010 11:14:37 +0200 From: Christophe Dutang To: r-devel@r-project.org Subject: [Rd] on the optim function Message-ID: <7e004a07-03e1-4ded-a506-6c564edb6...@gmail.com> Content-Type: text/plain; charset=us-ascii Dear useRs, I have just discovered that the R optim function does not return the number of iterations. I still wonder why line 632-634 of optim C, the iter variable is not returned (for the BFGS method for example) ? Is there any trick to compute the iteration number with function call number? Kind regards Christophe -- Christophe Dutang Ph.D. student at ISFA, Lyon, France website: http://dutangc.free.fr For BFGS the number of iterations is the number of gradient evaluations i.e., counts[2] For most of the optimization tools available for R (not just in optim), these counts are a nightmare for developers who want some consistent naming and meaning, as Ravi Varadhan and I can attest in our efforts to build the optimx() wrapper. It can also be debated for some methods what constitutes an "iteration". Best, JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] transpose of complex matrices in R
When one is working with complex matrices, "transpose" very nearly always means *Hermitian* transpose, that is, A[i,j] <- Conj(A[j,i]). One often writes A^* for the Hermitian transpose. I believe that I have actually (many years ago) used a true complex transpose, but I agree that one more often needs the conjugate transpose. I would not be in favour of changing t() because I feel transpose means transpose -- certainly where there are non-square matrices. But I'd be in favour of adding a conjt() (or some similar function) that does the conjugate transpose efficiently. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Location of source code for readline()
A few days ago on R-help I asked about a cross-platform timeout version of readline(). Some suggestions, but only partial joy so far. I can get the Gnu bash 'read -t ...' to work in Windows by using the 'bash -c ' construct, but then R's system() function does not seem to allow this to pass through. Similarly a Perl and Free Pascal routine that I tried, the latter being a single executable that did the prompt and the timeout. (I can send code offline if anyone interested -- not fully protected against bad inputs, however.) Now I'm wondering where the code for readline is located in the R source. I've tracked as far as the 'do_readln' in names.c, but now want to find the actual code to see if I can patch it, though I am a real novice in C. Suggestions welcome. My application, FYI, is to have a script that will display something, and wait for a keypress (for readline it seems to need the Enter key) but timeout after a preset number of seconds. The setTimeLimit "almost" works -- but not for readline. I'm thinking of a modified readline like readline(prompt='Do you want to continue?', timeout=6). Note that the issue seems to be Windows. I haven't a Mac to try, but Linux can be made to function by various methods at the top. Sigh. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] R CMD check issue with soft-linked directory
Good catch Simon. Changed fstab to force exec on mount of the drive in question and things worked. Thanks, JN Simon Urbanek wrote: On Mar 10, 2010, at 12:43 , Prof. John C Nash wrote: I've been having some strange problems with R CMD check in the last couple of days, but now believe I have localized the issue. I had been running Ubuntu Hardy on one drive and then upgraded to Jaunty, but put Jaunty on a different drive. I continue to be able to boot Hardy when I wish. I soft-linked my R working area i.e., /home/john/Rstuff >/media/lnx804/home/john/Rstuff I can still build packages fine, but get j...@nsrv-jaunty:~/jtest$ R CMD check minqa * checking for working pdflatex ... OK * using log directory '/media/store2/jn/test/minqa.Rcheck' * using R version 2.10.1 (2009-12-14) * using session charset: UTF-8 * checking for file 'minqa/DESCRIPTION' ... OK * checking extension type ... Package * this is package 'minqa' version '1.02' * checking package name space information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking for executable files ... OK * checking whether package 'minqa' can be installed ... OK * checking package directory ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... ERROR Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared library '/media/store2/jn/test/minqa.Rcheck/minqa/libs/minqa.so': /media/store2/jn/test/minqa.Rcheck/minqa/libs/minqa.so: failed to map segment from shared object: Operation not permitted Error : .onLoad failed in 'loadNamespace' for 'minqa' Error: package/namespace load failed for 'minqa' Execution halted It looks like this package has a loading problem: see the messages for details. The above test was run with a newly created jtest/ directory on yet another drive (partition) to see if there was possibly corruption. However, the issue seems to be one of using a softlink. I would be very surprised if this was softlink's fault per se. It looks more like a permission issue in your system to me -- the first thing I would check is that your /media/... is not mounted with noexec... Cheers, Simon I would not be upset if R CMD check simply told me that this isn't the right way to do things i.e., don't use the softlinked directory. Does anyone know if this is a recognized issue and is it possible to put some sort of reasonable error message into R CMD check? 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
[Rd] R CMD check issue with soft-linked directory
I've been having some strange problems with R CMD check in the last couple of days, but now believe I have localized the issue. I had been running Ubuntu Hardy on one drive and then upgraded to Jaunty, but put Jaunty on a different drive. I continue to be able to boot Hardy when I wish. I soft-linked my R working area i.e., /home/john/Rstuff >/media/lnx804/home/john/Rstuff I can still build packages fine, but get j...@nsrv-jaunty:~/jtest$ R CMD check minqa * checking for working pdflatex ... OK * using log directory '/media/store2/jn/test/minqa.Rcheck' * using R version 2.10.1 (2009-12-14) * using session charset: UTF-8 * checking for file 'minqa/DESCRIPTION' ... OK * checking extension type ... Package * this is package 'minqa' version '1.02' * checking package name space information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking for executable files ... OK * checking whether package 'minqa' can be installed ... OK * checking package directory ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... ERROR Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared library '/media/store2/jn/test/minqa.Rcheck/minqa/libs/minqa.so': /media/store2/jn/test/minqa.Rcheck/minqa/libs/minqa.so: failed to map segment from shared object: Operation not permitted Error : .onLoad failed in 'loadNamespace' for 'minqa' Error: package/namespace load failed for 'minqa' Execution halted It looks like this package has a loading problem: see the messages for details. The above test was run with a newly created jtest/ directory on yet another drive (partition) to see if there was possibly corruption. However, the issue seems to be one of using a softlink. I would not be upset if R CMD check simply told me that this isn't the right way to do things i.e., don't use the softlinked directory. Does anyone know if this is a recognized issue and is it possible to put some sort of reasonable error message into R CMD check? JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Idea for install.packages on (certain) linux distributions
Some time ago, I had some email discussion with Dirk E. about putting a front-end on
install.packages to first look at the debian repositories for R and use them before trying
to install from source. The code for this would not be very large. As Uwe points out in
another posting, the issue then becomes one of repository maintenance. And as more types
of installers get included, the code and the chance of package mismatches get more risky.
However, where we have repositories, it may be useful to have example code to try such an
approach. In pseudo-code this could be implemented without damaging install.packages as:
[start of my.install.packages]
if (exist(local.install) ) {
local.install(package)
} else {
install.packages(package)
}
If anyone gets enthused about this, I'd suggest posting on R-wiki. Note that local.install
will have to be pretty aggressive at checking the OS version etc.
JN
Date: Mon, 01 Mar 2010 16:33:08 +0100
From: "Carlos J. Gil Bellosta"
To: r-devel@r-project.org
Subject: [Rd] Idea for install.packages on (certain) linux
distributions
Message-ID: <4b8bde34.6070...@datanalytics.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hello,
I do not know whether this idea would be considered useful or not. Or
easy to implement.
But it is the case that on certain Linux distributions there are OS
packages with precompiled R packages. However, install.packages (and
related functions) download the source ones.
Would it be possible to add an extra "repository" (or option) on
install.packages that would direct R to use the OS level package manager
(apt-get, yum or the like) so as to install the precompiled packages
from the distribution mirrors instead of the CRAN ones?
Best regards,
Carlos J. Gil Bellosta
http://www.datanalytics.com
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] [R] R CMD check when package directory is symlinked
I've done a couple of searches and not found mention of the issue below, though some older scripts mention getting absolute paths for R CMD check. If the issue is "new" or unfamiliar I'll be happy to follow up and document it, but suspect it is in some sense already known and I've missed the right search strategy. The workaround is pretty simple -- copy the files to a new absolute directory. However, if I can save some grief for others, I will. I run Ubuntu Jaunty 9.04, but until recently was running Hardy 8.04. I symlinked my R-optimtest directory from /home/john/ to my "old" home directory's version. When I changed a package and tried R CMD check I got the following. * checking whether the package can be loaded ... ERROR Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared library '/media/lnx804/home/john/R-optimtest/work/minqa.Rcheck/minqa/libs/minqa.so': /media/lnx804/home/john/R-optimtest/work/minqa.Rcheck/minqa/libs/minqa.so: failed to map segment from shared object: Operation not permitted Error : .onLoad failed in 'loadNamespace' for 'minqa' Error: package/namespace load failed for 'minqa' Execution halted I copied the entire minqa tree to a new "localR" directory and things work fine. So it looks like R is unhappy with the expanded directory name i.e., from /home/john/R-optimtest/work/... to /media/lnx804/home/john/R-optimtest/work/... Cheers, JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Buglet in optim() SANN
Indeed Brian is correct about the functioning of SANN and the R
documentation. I'd misread the "maxit" warning. Things can stay as they
are for now.
The rest of this msg is for information and an invitation to off-list
discussion.
I realize my posting opens up the can of worms about what "convergence"
means. As someone who has occasionally published discussions on
convergence versus termination, I'd certainly prefer to set the
'convergence' flag to 1 for SANN, since we have only a termination at
the maximum number of function evaluations and not necessarily a result
that can be presumed to be "optim"al. Or perhaps put a note in the
description of the 'convergence' flag to indicate the potential
misinterpretation with SANN where results need the user to externally
check if they are likely to be usable as an optimum.
It may be better to call the non-zero results for "convergence" a
"termination indicator" rather than an "error code". Some related
packages like ucminf give more than one non-zero indicators for results
that are generally usable as optima. They are informational rather than
errors. Writing our optimx wrapper for a number of methods has forced us
to think about how such information is returned and reported through a
flag like "convergence". There are several choices and plenty of room
for confusion.
Right now a few of us are working on improvements for optimization, but
the first goal is to get things working OK for smooth, precisely defined
functions. However, we have been looking at methods like SANN for
multimodal and noisy (i.e., imprecisely defined) functions. For those
problems, knowing when you have an acceptable or usable result is never
easy.
Comments and exchanges welcome -- off-list of course.
Cheers, JN
Prof Brian Ripley wrote:
> As the posting guide says, please read the help carefully before
> posting. It does say:
>
> ‘maxit’ The maximum number of iterations. Defaults to ‘100’ for
> the derivative-based methods, and ‘500’ for ‘"Nelder-Mead"’.
> For ‘"SANN"’ ‘maxit’ gives the total number of function
> evaluations. There is no other stopping criterion.
>^^^^
> Defaults to ‘1’.
>
> so this is indicating 'successful convergence' as documented.
>
> On Tue, 20 Oct 2009, Prof. John C Nash wrote:
>
>> I think SANN method in optim() is failing to report that it has not
>> converged. Here is an example
>>
>> genrose.f<- function(x, gs=NULL){ # objective function
>> ## One generalization of the Rosenbrock banana valley function (n
>> parameters)
>> n <- length(x)
>>if(is.null(gs)) { gs=100.0 }
>> fval<-1.0 + sum (gs*(x[1:(n-1)]^2 - x[2:n])^2 + (x[2:n] - 1)^2)
>>return(fval)
>> }
>>
>> xx<-rep(pi,10)
>> test<-optim(xx,genrose.f,method="SANN",control=list(maxit=1000,trace=1))
>> print(test)
>>
>>
>> Output is:
>>
>>> source("tsann.R")
>> sann objective function values
>> initial value 40781.805639
>> iter 999 value 29.969529
>> final value 29.969529
>> sann stopped after 999 iterations
>> $par
>> [1] 1.0135254 0.9886862 1.1348609 1.0798927 1.0327997 1.1087146 1.1642130
>> [8] 1.3038754 1.8628391 3.7569285
>>
>> $value
>> [1] 29.96953
>>
>> $counts
>> function gradient
>>1000 NA
>>
>> $convergence
>> [1] 0 <-- THIS SHOULD BE 1 ACCORDING TO THE DOCS
>
> It _should_ be 0 according to the help page.
>
>>
>> $message
>> NULL
>>
>> Note terribly important, but maybe fixable.
>>
>> Cheers,
>>
>> John Nash
>>
>> __
>> 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] Advice on how to arrange fix of buglet
Recently I reported a small bug in optim's SANN method failing to report that it had exceeded the maximum function evaluation limit in the convergence code. This is a small enough matter that I was reluctant to create a full-blown bug report. Indeed in the optimx package Ravi Varadhan and I have been developing on r-forge (under the OptimizeR project) it was a minimal work around to fix the matter in our wrapper that incorporates optim() and a number of other tools. While I don't normally do C code, I could likely figure out a fix for optim too. My query is about how to best get this done without causing a lot of work for others i.e., where to I send patches etc. I expect there are a number of similar issues for different areas of R and its documentation, and a clarification from someone in the core team could streamline things. Maybe the bug system is still the right place? Cheers, JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Buglet in optim() SANN
I think SANN method in optim() is failing to report that it has not
converged. Here is an example
genrose.f<- function(x, gs=NULL){ # objective function
## One generalization of the Rosenbrock banana valley function (n
parameters)
n <- length(x)
if(is.null(gs)) { gs=100.0 }
fval<-1.0 + sum (gs*(x[1:(n-1)]^2 - x[2:n])^2 + (x[2:n] - 1)^2)
return(fval)
}
xx<-rep(pi,10)
test<-optim(xx,genrose.f,method="SANN",control=list(maxit=1000,trace=1))
print(test)
Output is:
> source("tsann.R")
sann objective function values
initial value 40781.805639
iter 999 value 29.969529
final value 29.969529
sann stopped after 999 iterations
$par
[1] 1.0135254 0.9886862 1.1348609 1.0798927 1.0327997 1.1087146 1.1642130
[8] 1.3038754 1.8628391 3.7569285
$value
[1] 29.96953
$counts
function gradient
1000 NA
$convergence
[1] 0 <-- THIS SHOULD BE 1 ACCORDING TO THE DOCS
$message
NULL
Note terribly important, but maybe fixable.
Cheers,
John Nash
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Why license filtering is important (was non-GPL ...)
Using the acronym "GPL" in comments on the management of package repositories led the discussion well away from the issues I wanted to shed light upon, so I have changed the subject tag. Examples of my concerns: 1) I use a package with a non-free component and learn how to work with it efficiently. I'm a retired academic, but still try to recover some costs of the work I do for R and other open projects by doing some contract and training work. A client has a problem and wants a solution. If there is a "non-commercial" restriction, I'm put in the awkward position of having to get a permission -- possibly from a defunct organization or one that has long forgotten how to grant such a permission. 2) There may be new tools available that ARE free software. Dirk's blacklist (perhaps we could have a less loaded name?) may suggest opportunities for gradually moving to packages that can be used without need to check license details. I have used such tasks in the past as student projects where they are relatively straightforward. 3) Many workers are not aware of the consequences of license status of their codes (I was not for some years). The development of CRAN and similar repositories has been and can be a positive force allowing for improvement and understanding of methods. 4) We definitely should retain the ability to access non-free codes -- somehow folk have misread my comments otherwise. I'll use them for learning and when there is no alternative, but if at all possible, I'll choose the free ones for production work so I don't get caught out as above. A comment: In looking at SparseM, I first used the pdf -- it simply says GPL(>=2) as the license. (I'm sure I've got similar bugs in my own docs.) I dug deeper and found the LICENSE file, and also looked at cholesky.f, which is unfortunately one of the bigger files. I was hoping it was somewhat related to work I've done over the years on Cholesky decompositions in the hope I could offer a substitute as per concern (2) above, but it is not of the same nature as the algorithms I worked on as far as I can determine. Maybe someone else is able to step in on that. And a big tip of the hat to Dirk E. and Charles B. for the cran2deb work -- I already admired the work. Now I've started looking at some of the package files for license info, I'm amazed they've got as far as they have. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Non-GPL packages for R
License filters will work for me. My offer stands to help on documentation,or to act as a "stooge" to test tools in this area. Thanks to those who responded. And for myself, "GPL compatible" was my intended expression. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Non-GPL packages for R
The responses to my posting yesterday seem to indicate more consensus than I expected: 1) CRAN should be restricted to GPL-equivalent licensed packages 2) r-forge could be left "buyer beware" using DESCRIPTION information 3) We may want a specific repository for restricted packages (RANC?) How to proceed? A short search on Rseek did not turn up a chain of command for CRAN. I'm prepared to help out with documentation etc. to move changes forward. They are not, in my opinion, likely to cause a lot of trouble for most users, and should simplify things over time. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Non-GPL packages for R
Subject: Non-GPL packages for R Packages that are not licensed in a way that permits re-distribution on CRAN are frequently a source of comment and concern on R-help and other lists. A good example of this problem is the Rdonlp2 package that has caused a lot of annoyance for a number of optimization users in R. They are also an issue for efforts like Dirk Eddelbuettel's cran2deb. There are, however, a number of circumstances where non-GPL equivalent packages may be important to users. This can imply that users need to both install an R package and one or more dependencies that must be separately obtained and licensed. One such situation is where a new program is still under development and the license is not clear, as in the recent work we pursued with respect to Mike Powell's BOBYQA. We wanted to verify if this were useful before we considered distribution, and Powell had been offering copies of his code on request. Thus we could experiment, but not redistribute. Recently Powell's approval to redistribute has been obtained. We believe that it is important that non-redistributable codes be excluded from CRAN, but that they could be available on a repository such as r-forge. However, we would like to see a clearer indication of the license status on r-forge. One possibility is an inclusion of a statement and/or icon indicating such status e.g., green for GPL or equivalent, amber for uncertain, red for restricted. Another may be a division of directories, so that GPL-equivalent packages are kept separate from uncertain or restricted licensed ones. We welcome comments and suggestions on both the concept and the technicalities. John Nash & Ravi Varadhan __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] identical(0, -0)
I'll save space and not include previous messages. My 2 cents: At the very least the documentation needs a fix. If it is easy to do, then Ted Harding's suggestion of a switch (default OFF) to check for sign difference would be sensible. I would urge inclusion in the documentation of the +0, -0 example(s) if there is NOT a way in R to distinguish these. There are occasions where it is useful to be able to detect things like this (and NaN and Inf and -Inf etc.). They are usually not of interest to users, but sometimes are needed for developers to check edge effects. For those cases it may be time to consider a package FPIEEE754 or some similar name to allow testing and possibly setting of flags for some of the fancier features. Likely used by just a few of us in extreme situations. Unfortunately, some of the nice exception handling that was suggested in the standard is apparently rarely implemented in compilers. For info, I was actually a voting member of IEEE 754 because I found a nice "feature" in the IBM Fortran G arithmetic, though I never sat in on any of the meetings. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Automatic Differentiation for R
Martin (see below) gives a good explanation of the difference between AD and
symbolic
differentiations. I'm of the opinion we can use both. However, the real issue
as far
as I'm concerned (from an optimizer's point of view, which may also be that of
ODE and
PDE folk) is that right now none of the offerings that we have are easy to use.
Indeed,
usability is one of the key issues in my current efforts to improve
optimization tools.
There are several people in UK from both AD and R sides who are communicating,
an initiative that Shaun F. helped get going. I'm copying the folk I think are
involved -- forgive omissions and pass things along please. The main purpose of
this
message is to ensure Martin's offer is noted, as in my opinion his knowledge of
the
R internals is very valuable.
Cheers, JN
[MM stumbling over on old thread ... he'd be interested]
> "GaGr" == Gabor Grothendieck
> on Wed, 15 Apr 2009 09:53:18 -0400 writes:
GaGr> Not sure if this is sufficient for your needs but R does include
symbolic
GaGr> differentiation, see ?D, and the Ryacas and rSymPy
GaGr> packages interface R to the yacas and sympy computer algebra
GaGr> systems (CAS) and those system include symbolic differentiation.
No, symbolic differentiation is not enough.
Automatic Differentiation (AD) is something much more general (in one
way) and much less mathematical from a classical view point:
But then, AD is much more generally useful for minimization as, basically,
the input is an R function
f(x) {with x multidimensional}
or f(x1,x2, ..., xp) {with scalar x1, x2, ..}
and the output is again an R function
which computes f() and all {or just selected} partial
derivatives d f / d{xi}.
Now consider that the function f() can contain if() and while()
clauses and conceptually ever language feature of R.
In practice, I'm pretty sure the list of features would have to
be restricted, similarly as they'd have to for an R compiler to
be feasible.
I agree that AD for R would be very nice and could be very
useful.
I'd also be interested to help AD people learn the S4 classes
and methods (hoping that it's close enough to what they call
"operator overloading" something I'd presume to be less general
than the powerful S4 class/methods system).
Martin Maechler, ETH Zurich
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] nlminb - Port library codes
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
[Rd] (Not quite) Automatic Differentiation
Many thanks to Gabor Grothendieck for responding to my posting about Automatic Differentiation (invite from Shaun Forth for interaction with R developers) showing how one might use rSymPy and symbolic (rather than automatic) differentiation to get a function that computes gradients. See http://code.google.com/p/rsympy/#Automatic_Differentiation_(well,_sort_of) for a worked example on the Broyden test function. This is a big step forward. There's still a way to go before we can produce a vectorized gradient code automatically when the size of the problem is variable, but the example may serve to incite some imaginative coders to action. Thanks again Gabor. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Automatic Differentiation for R
In efforts to improve optimization tools for R, one of my interests has been getting automatic differentiation capabilities so that analytic rather than numerical derivatives can be used. They would be helpful in several other areas besides optimization, My timings show factors of the order of 1000s in time improvements by avoiding numerical derivatives in some cases. There has been some work in this e.g., http://code.google.com/p/pbs-software/ is an R interface to ADMB (Automatic Differentiation Model Builder). However, as far as I can see, this is directed essentially to nonlinear least squares modelling, an important but not general problem. Tom Coleman of Waterloo responded favourably with some advice, but the most enthusiastic answer came from Shaun Forth, which I have included below. I read this as an opportunity to develop what could be a profitable collaboration with the AD community. Unfortunately, I cannot take up the invitation to join the AD folk in Oxford due to a pre-existing obligation. Nor am I more than a complete novice with S3 and S4 classes etc. I am, nevertheless, willing to help organize the effort e.g., do some of the communications, chasing grant money, getting Google Summer of Code applications filled in etc. Can the R community come up with a few people who can provide the AD workers with appropriate information? If so, is there a reasonable chance to generate sufficient funding for a student? I suspect the answer in both cases is yes, but that we need some form of "booster cables" to get things going. (In Canada, booster cables are used to get cars started in winter by connecting a running vehicle's battery to that of a dead one.) I suggest communications off-list until there is progress to report. Possibly there is a better forum for this -- suggestions welcome. John Nash included msg from Shaun Forth --- Hi John, My computational statistics colleague Trevor Ringrose has asked me to consider AD in R in the past. As you may or may not be aware AD is implemented in one of two ways: overloading using OO features of the target language, or source transformation using compiler tools (after several man years of development) to read in the target code and spit out the differentiated code. Last time I looked I didn't think the object oriented features of R were up to overloading but on checking today I can see that this might now be possible (I can see overloading of arithmetic operators and functions for example now which I didn't see last time). I'd certainly be interested in following this up particularly on the overloading side but would need to get funding for a PhD student to do the graft. It would be particularly interesting doing this in an open source language because we could then perhaps tweak some of the core language features if they didn't seamlessly support the AD (we can't do this in Matlab and that is a pain!). My immediate suggestion is that you, or some other more local (to UK) R expert talks at the next European AD workshop in Oxford http://www.autodiff.org/?module=Workshops&submenu=EuroAD/8/main We're a very friendly group and I'm sure there are others who might like to tackle R or perhaps we could put together a multigroup project. If someone could give a talk on R, its language features including the OO aspects, and some optimisation examples with associated code, the group there would be able to give you the best feedback on the planet on the possibilities. Please do treat this as a positive response and let's keep in touch on this. Regards Shaun Dr Shaun Forth Applied Mathematics & Scientific Computation Cranfield Defence and Security Cranfield University, Shrivenham Campus Swindon SN6 8LA, England tel: +44 (0)1793 785311 fax: +44 (0)1793 784196 email: s.a.fo...@cranfield.ac.uk http://www.amorg.co.uk # __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
