[Rd] Minor glitch in optim()

[J C Nash]

Having worked with optim() and related programs for years, it surprised me
that I haven't noticed this before, but optim() is inconsistent in how it
deals with bounds constraints specified at infinity. Here's an example:

# optim-glitch-Ex.R
x0<-c(1,2,3,4)
fnt <- function(x, fscale=10){
  yy <- length(x):1
  val <- sum((yy*x)^2)*fscale
}
grt <- function(x, fscale=10){
  nn <- length(x)
  yy <- nn:1
  #gg <- rep(NA,nn)
  gg <- 2*(yy^2)*x*fscale
  gg
}

npar <- 4
lower <- -Inf
l2 <- rep(-Inf,npar)

a1 <- optim(x0, fnt, grt, lower=lower, method="BFGS") # works
a1
a1a<- optim(x0, fnt, grt, lower=l2, method="BFGS") # changes method!
a1a

The first call uses BFGS method without warning. The second gives
a warning that L-BFGS-B should be used, and from the output uses
this.

This is a bit of an edge case. My own preference would be for optim()
to simply fail if bounds of any type are specified without L-BFGS-B
as the method. I believe that gives clarity, even though infinite
bounds imply an unconstrained problem.

The behaviour where a scalar infinite bound is treated as unconstrained
but a vector is not is inconsistent, however, and I think that at some
point should be fixed. Possibly the easiest way is to treat infinite
bounds specified as a vector the same as those specified as a scalar.
That is to adjust the code in File src/library/stats/R/optim.R
in the block

if((length(lower) > 1L || length(upper) > 1L ||
   lower[1L] != -Inf || upper[1L] != Inf)
   && !any(method == c("L-BFGS-B","Brent"))) {
warning("bounds can only be used with method L-BFGS-B (or Brent)")
method <- "L-BFGS-B"
}

Possibly

if((any(is.finite(lower) || any(is.finite(upper))
   && !any(method == c("L-BFGS-B","Brent"))) {
warning("bounds can only be used with method L-BFGS-B (or Brent)")
method <- "L-BFGS-B"
}

Best, JN

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

[Rd] histoRicalg -- project to document older methods used by R and transfer knowledge

[J C Nash]

After some thought, I decided r-devel was probably the best of the R lists
for this item. Do feel free to share, as the purpose is to improve documentation
and identify potential issues.

John Nash



The R Consortium has awarded some modest funding for "histoRicalg",
a project to document and transfer knowledge of some older algorithms
used by R and by other computational systems. These older codes
are mainly in Fortran, but some are in C, with the original implementations
possibly in other programming languages. My efforts
were prompted by finding some apparent bugs in codes, which could be either
from the original programs or else the implementations. Two examples
in particular -- in nlm() and in optim::L-BFGS-B -- gave impetus
to the project.

As a first task, I am hoping to establish a "Working Group on
Algorithms Used in R" to identify and prioritize issues and to
develop procedures for linking older and younger workers to enable
the transfer of knowledge. Expressions of interest are welcome,
either to me (nashjc _at_ uottawa.ca) or to the mailing list
(https://lists.r-consortium.org/g/rconsortium-project-histoRicalg).
A preliminary web-site is at https://gitlab.com/nashjc/histoRicalg.

While active membership of the Working Group is desirable, given
the nature of this project, I anticipate that most members will
contribute mainly by providing timely and pertinent ideas. Some
may not even be R users, since the underlying algorithms are used
by other computing systems and the documentation effort has many
common features. We will also need participation of younger
workers willing to learn about the methods that underly the
computations in R.

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

[Rd] trace in uniroot() ?

[J C Nash]

In looking at rootfinding for the histoRicalg project (see 
gitlab.com/nashjc/histoRicalg),
I thought I would check how uniroot() solves some problems. The following short 
example

ff <- function(x){ exp(0.5*x) - 2 }
ff(2)
ff(1)
uniroot(ff, 0, 10)
uniroot(ff, c(0, 10), trace=1)
uniroot(ff, c(0, 10), trace=TRUE)


shows that the trace parameter, as described in the Rd file, does not seem to
be functional except in limited situations (and it suggests an
integer, then uses a logical for the example, e.g.,
 ## numerically, f(-|M|) becomes zero :
 u3 <- uniroot(exp, c(0,2), extendInt="yes", trace=TRUE)
)

When extendInt is set, then there is some information output, but trace alone
produces nothing.

I looked at the source code -- it is in R-3.5.1/src/library/stats/R/nlm.R and
calls zeroin2 code from R-3.5.1/src/library/stats/src/optimize.c as far as I
can determing. My code inspection suggests trace does not show the iterations
of the rootfinding, and only has effect when the search interval is allowed
to be extended. It does not appear that there is any mechanism to ask
the zeroin2 C code to display intermediate work.

This isn't desperately important for me as I wrote an R version of the code in
package rootoned on R-forge (which Martin Maechler adapted as unirootR.R in
Rmpfr so multi-precision roots can be found). My zeroin.R has 'trace' to get
the pattern of different steps. In fact it is a bit excessive. Note
unirootR.R uses 'verbose' rather than 'trace'. However, it would be nice to be
able to see what is going on with uniroot() to verify equivalent operation at
the same precision level. It is very easy for codes to be very slightly
different and give quite widely different output.

Indeed, even without the trace, we see (zeroin from rootoned here)

> zeroin(ff, c(0, 10), trace=FALSE)
$root
[1] 1.386294

$froot
[1] -5.658169e-10

$rtol
[1] 7.450581e-09

$maxit
[1] 9

> uniroot(ff, c(0, 10), trace=FALSE)
$root
[1] 1.38629

$f.root
[1] -4.66072e-06

$iter
[1] 10

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

>

Is the lack of trace a bug, or at least an oversight? Being able to follow 
iterations is a
classic approach to checking that computations are proceeding as they should.

Best, JN

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

[Rd] typo in Ubuntu download page

[J C Nash]

In https://cran.r-project.org/bin/linux/ubuntu/

Administration and Maintances of R Packages
   ^^

Minor stuff, but if someone who can edit is on the page,
perhaps it can be changed to "Maintenance"

Best, JN

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

Re: [Rd] trace in uniroot() ?

[J C Nash]

Despite my years with R, I didn't know about trace(). Thanks.

However, my decades in the minimization and root finding game make me like 
having
a trace that gives some info on the operation, the argument and the current 
function value.
I've usually found glitches are a result of things like >= rather than > in 
tests etc., and
knowing what was done is the quickest way to get there.

This is, of course, the numerical software developer view. I know "users" (a 
far too vague
term) don't like such output. I've sometimes been tempted with my svd or 
optimization codes to
have a return message in bold-caps "YOUR ANSWER IS WRONG AND THERE'S A LAWYER 
WAITING TO
MAKE YOU PAY", but I usually just satisfy myself with "Not at a minimum/root".

Best, JN

On 2018-08-13 06:00 PM, William Dunlap wrote:
> I tend to avoid the the trace/verbose arguments for the various root finders 
> and optimizers and instead use the trace
> function or otherwise modify the function handed to the operator.  You can 
> print or plot the arguments or save them.  E.g.,
> 
>> trace(ff, print=FALSE, quote(cat("x=", deparse(x), "\n", sep="")))
> [1] "ff"
>> ff0 <- uniroot(ff, c(0, 10))
> x=0
> x=10
> x=0.0678365490630423
> x=5.03391827453152
> x=0.490045026724842
> x=2.76198165062818
> x=1.09760394309444
> x=1.92979279686131
> x=1.34802524899502
> x=1.38677998493585
> x=1.3862897003949
> x=1.38635073555115
> x=1.3862897003949
> 
> or
> 
>> X <- numeric()
>> trace(ff, print=FALSE, quote(X[[length(X)+1]] <<- x))
> [1] "ff"
>> ff0 <- uniroot(ff, c(0, 10))
>> X
>  [1]  0. 10.  0.06783655
>  [4]  5.03391827  0.49004503  2.76198165
>  [7]  1.09760394  1.92979280  1.34802525
> [10]  1.38677998  1.38628970  1.38635074
> [13]  1.38628970
> 
> This will not tell you why the objective function is being called (e.g. in a 
> line search
> or in derivative estimation), but some plotting or other postprocessing can 
> ususally figure that out.
> 
> 
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com <http://tibco.com>
> 
> On Mon, Jul 30, 2018 at 11:35 AM, J C Nash  <mailto:profjcn...@gmail.com>> wrote:
> 
> In looking at rootfinding for the histoRicalg project (see 
> gitlab.com/nashjc/histoRicalg
> <http://gitlab.com/nashjc/histoRicalg>),
> I thought I would check how uniroot() solves some problems. The following 
> short example
> 
> ff <- function(x){ exp(0.5*x) - 2 }
> ff(2)
> ff(1)
> uniroot(ff, 0, 10)
> uniroot(ff, c(0, 10), trace=1)
> uniroot(ff, c(0, 10), trace=TRUE)
> 
> 
> shows that the trace parameter, as described in the Rd file, does not 
> seem to
> be functional except in limited situations (and it suggests an
> integer, then uses a logical for the example, e.g.,
>  ## numerically, f(-|M|) becomes zero :
>      u3 <- uniroot(exp, c(0,2), extendInt="yes", trace=TRUE)
> )
> 
> When extendInt is set, then there is some information output, but trace 
> alone
> produces nothing.
> 
> I looked at the source code -- it is in R-3.5.1/src/library/stats/R/nlm.R 
> and
> calls zeroin2 code from R-3.5.1/src/library/stats/src/optimize.c as far 
> as I
> can determing. My code inspection suggests trace does not show the 
> iterations
> of the rootfinding, and only has effect when the search interval is 
> allowed
> to be extended. It does not appear that there is any mechanism to ask
> the zeroin2 C code to display intermediate work.
> 
> This isn't desperately important for me as I wrote an R version of the 
> code in
> package rootoned on R-forge (which Martin Maechler adapted as unirootR.R 
> in
> Rmpfr so multi-precision roots can be found). My zeroin.R has 'trace' to 
> get
> the pattern of different steps. In fact it is a bit excessive. Note
> unirootR.R uses 'verbose' rather than 'trace'. However, it would be nice 
> to be
> able to see what is going on with uniroot() to verify equivalent 
> operation at
> the same precision level. It is very easy for codes to be very slightly
> different and give quite widely different output.
> 
> Indeed, even without the trace, we see (zeroin from rootoned here)
> 
> > zeroin(ff, c(0, 10), trace=FALSE)
> $root
> [1] 1.386294
> 
> $froot
> [1] -5.658169e-10
> 
> $rtol
> [1] 7.450581e-09
> 
> $maxit
> [1] 9
> 
> > uniroot(ff, c(0, 10), trace=FALSE)
> $root
> 

Re: [Rd] Package inclusion in R core implementation

[J C Nash]

As the original coder (in mid 1970s) of BFGS, CG and Nelder-Mead in optim(), 
I've
been pushing for some time for their deprecation. They aren't "bad", but we have
better tools, and they are in CRAN packages. Similarly, I believe other 
optimization
tools in the core (optim::L-BFGS-B, nlm, nlminb) can and should be moved to
packages (there are already 2 versions at least of LBFGS that I and Matt Fidler
are merging). And optim::SANN does not match any usual expectations of users.

I'm sure there are other tools for other tasks that can and should move to 
packages
to streamline the work of our core team. However, I can understand that there 
is this
awkward issue of actually doing this. I know I'm willing to help with preparing
"Transition Guide" documentation and scripts, and would be surprised if there 
are
not others. R already has a policy of full support only for current version, so
hanging on to antique tools (the three codes at the top are based on papers all
of which now qualify for >50 years old) seems inconsistent with other messages.

For information: I'm coordinating a project to build understanding of what
older algorithms are in R as the histoRicalg project. See
https://gitlab.com/nashjc/histoRicalg. We welcome participation.

Best, JN

On 2019-03-04 7:59 a.m., Jim Hester wrote:
> Conversely, what is the process to remove a package from core R? It seems
> to me some (many?) of the packages included are there more out of
> historical accident rather than any technical need to be in the core
> distribution. Having them as a core (or recommended) package makes them
> harder update independently to R and makes testing, development and
> contribution more cumbersome.
> 
> On Fri, Mar 1, 2019 at 4:35 AM Morgan Morgan 
> wrote:
> 
>> Hi,
>>
>> It sometimes happens that some packages get included to R like for example
>> the parallel package.
>>
>> I was wondering if there is a process to decide whether or not to include a
>> package in the core implementation of R?
>>
>> For example, why not include the Rcpp package, which became for a lot of
>> user the main tool to extend R?
>>
>> What is our view on the (not so well known) dotCall64 package which is an
>> interesting alternative for extending R?
>>
>> Thank you
>> Best regards,
>> Morgan
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> 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
>

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

Re: [Rd] Package inclusion in R core implementation

[J C Nash]

I concur with Avraham that capabilities need to be ensured e.g., in recommended
packages. I should have mentioned that. My concern is that the core should be
focused on the programming language aspects. The computational math and some of 
the more
intricate data management could better be handled by folk outside the core.

JN

On 2019-03-04 9:12 a.m., Avraham Adler wrote:
> On Mon, Mar 4, 2019 at 5:01 PM J C Nash  <mailto:profjcn...@gmail.com>> wrote:
> 
> As the original coder (in mid 1970s) of BFGS, CG and Nelder-Mead in 
> optim(), I've
> been pushing for some time for their deprecation. They aren't "bad", but 
> we have
> better tools, and they are in CRAN packages. Similarly, I believe other 
> optimization
> tools in the core (optim::L-BFGS-B, nlm, nlminb) can and should be moved 
> to
> packages (there are already 2 versions at least of LBFGS that I and Matt 
> Fidler
> are merging). And optim::SANN does not match any usual expectations of 
> users.
> 
> I'm sure there are other tools for other tasks that can and should move 
> to packages
> to streamline the work of our core team. However, I can understand that 
> there is this
> awkward issue of actually doing this. I know I'm willing to help with 
> preparing
> "Transition Guide" documentation and scripts, and would be surprised if 
> there are
> not others. R already has a policy of full support only for current 
> version, so
> hanging on to antique tools (the three codes at the top are based on 
> papers all
> of which now qualify for >50 years old) seems inconsistent with other 
> messages.
> 
> For information: I'm coordinating a project to build understanding of what
> older algorithms are in R as the histoRicalg project. See
> https://gitlab.com/nashjc/histoRicalg. We welcome participation.
> 
> Best, JN
> 
> On 2019-03-04 7:59 a.m., Jim Hester wrote:
> > Conversely, what is the process to remove a package from core R? It 
> seems
> > to me some (many?) of the packages included are there more out of
> > historical accident rather than any technical need to be in the core
> > distribution. Having them as a core (or recommended) package makes them
> > harder update independently to R and makes testing, development and
> > contribution more cumbersome.
> >
> > On Fri, Mar 1, 2019 at 4:35 AM Morgan Morgan  <mailto:morgan.email...@gmail.com>>
> > wrote:
> >
> >> Hi,
> >>
> >> It sometimes happens that some packages get included to R like for 
> example
> >> the parallel package.
> >>
> >> I was wondering if there is a process to decide whether or not to 
> include a
> >> package in the core implementation of R?
> >>
> >> For example, why not include the Rcpp package, which became for a lot 
> of
> >> user the main tool to extend R?
> >>
> >> What is our view on the (not so well known) dotCall64 package which is 
> an
> >> interesting alternative for extending R?
> >>
> >> Thank you
> >> Best regards,
> >> Morgan
> >>
> 
> 
> I have No arguments with updating code to more correct or modern versions, 
> but I think that as a design decision, base R
> should have optimization routines as opposed to it being an external package 
> which conceptually could be orphaned. Or at
> least some package gets made recommended and adopted by R core.
> 
> Thank you,
> 
> Avi
> -- 
> Sent from Gmail Mobile

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

[Rd] Fwd: Package inclusion in R core implementation

[J C Nash]

Rereading my post below, I realize scope for misinterpretation. As I have said 
earlier,
I recognize the workload in doing any streamlining, and also the immense 
service to us
all by r-core. The issue is how to manage the workload efficiently while 
maintaining
and modernizing the capability. That is likely as challenging as doing the work 
itself.

JN


> I concur with Avraham that capabilities need to be ensured e.g., in 
> recommended
> packages. I should have mentioned that. My concern is that the core should be
> focused on the programming language aspects. The computational math and some 
> of the more
> intricate data management could better be handled by folk outside the core.
> 
> JN
>

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

Re: [Rd] Bug in the "reformulate" function in stats package

[J C Nash]

The main thing is to post the "small reproducible example".

My (rather long term experience) can be written

  if (exists("reproducible example") ) {
 DeveloperFixHappens()
  } else {
 NULL
  }

JN

On 2019-03-29 11:38 a.m., Saren Tasciyan wrote:
> Well, first I can't sign in bugzilla myself, that is why I wrote here first. 
> Also, I don't know if I have the time at
> the moment to provide tests, multiple examples or more. If that is not ok or 
> welcomed, that is fine, I can come back,
> whenever I have more time to properly report the bug.
> 
> I didn't find the existing bug report, sorry for that.
> 
> Yes, it is related. My problem was that I have column names with spaces and 
> current solution doesn't solve it. I have a
> solution, which works for me and maybe also for others.
> 
> Either, someone can register me to bugzilla or I can post it here, which 
> could give some direction to developers. I
> don't mind whichever is preferred here.
> 
> Best,
> 
> Saren
> 
> 
> On 29.03.19 09:29, Martin Maechler wrote:
>>> Saren Tasciyan
>>>  on Thu, 28 Mar 2019 17:02:10 +0100 writes:
>>  > Hi,
>>  > I have found a bug in reformulate function and have a solution for 
>> it. I
>>  > was wondering, where I can submit it?
>>
>>  > Best,
>>  > Saren
>>
>>
>> Well, you could have given a small reproducible example
>> depicting the bug, notably when posting here:
>> Just a prose text with no R code or other technical content is
>> almost always not really appropriate fo the R-devel mailing list.
>>
>> Further, in such a case you should google a bit and hopefully
>> have found
>>     https://www.r-project.org/bugs.html
>>
>> which also mention reproducibility (and many more useful things).
>>
>> Then it also tells you about R's bug repository, also called
>> "R's bugzilla" at https://bugs.r-project.org/
>>
>> and if you are diligent (but here, I'd say bugzilla is
>> (configured?) far from ideal), you'd also find bug PR#17359
>>
>>     https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17359
>>
>> which was reported already on Nov 2017 .. and only fixed
>> yesterday (in the "cleanup old bugs" process that happens
>> often before the big new spring release of R).
>>
>> So is your bug the same as that one?
>>
>> Martin
>>
>>  > --
>>  > Saren Tasciyan
>>  > /PhD Student / Sixt Group/
>>  > Institute of Science and Technology Austria
>>  > Am Campus 1
>>  > 3400 Klosterneuburg, Austria
>>
>>  > __
>>  > 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

Re: [Rd] unstable corner of parameter space for qbeta?

[J C Nash]

Given that a number of us are housebound, it might be a good time to try to
improve the approximation. It's not an area where I have much expertise, but in
looking at the qbeta.c code I see a lot of root-finding, where I do have some
background. However, I'm very reluctant to work alone on this, and will ask
interested others to email off-list. If there are others, I'll report back.

Ben: Do you have an idea of parameter region where approximation is poor?
I think that it would be smart to focus on that to start with.

Martin: On a separate precision matter, did you get my query early in year 
about double
length accumulation of inner products of vectors in Rmpfr? R-help more or
less implied that Rmpfr does NOT use extra length. I've been using David
Smith's FM Fortran where the DOT_PRODUCT does use double length, but it
would be nice to have that in R. My attempts to find "easy" workarounds have
not been successful, but I'll admit that other things took precedence.

Best,

John Nash



On 2020-03-26 4:02 a.m., Martin Maechler wrote:
>> Ben Bolker 
>> on Wed, 25 Mar 2020 21:09:16 -0400 writes:
> 
> > I've discovered an infelicity (I guess) in qbeta(): it's not a bug,
> > since there's a clear warning about lack of convergence of the numerical
> > algorithm ("full precision may not have been achieved").  I can work
> > around this, but I'm curious why it happens and whether there's a better
> > workaround -- it doesn't seem to be in a particularly extreme corner of
> > parameter space. It happens, e.g., for  these parameters:
> 
> > phi <- 1.1
> > i <- 0.01
> > t <- 0.001
> > shape1 = i/phi  ##  0.009090909
> > shape2 = (1-i)/phi  ## 0.9
> > qbeta(t,shape1,shape2)  ##  5.562685e-309
> > ##  brute-force uniroot() version, see below
> > Qbeta0(t,shape1,shape2)  ## 0.9262824
> 
> > The qbeta code is pretty scary to read: the warning "full precision
> > may not have been achieved" is triggered here:
> 
> > 
> https://github.com/wch/r-source/blob/f8d4d7d48051860cc695b99db9be9cf439aee743/src/nmath/qbeta.c#L530
> 
> > Any thoughts?
> 
> Well,  qbeta() is mostly based on inverting pbeta()  and pbeta()
> has *several* "dangerous" corners in its parameter spaces
> {in some cases, it makes sense to look at the 4 different cases
>  log.p = TRUE/FALSE  //  lower.tail = TRUE/FALSE  separately ..}
> 
> pbeta() itself is based on the most complex numerical code in
> all of base R, i.e., src/nmath/toms708.c  and that algorithm
> (TOMS 708) had been sophisticated already when it was published,
> and it has been improved and tweaked several times since being
> part of R, notably for the log.p=TRUE case which had not been in
> the focus of the publication and its algorithm.
> [[ NB: part of this you can read when reading  help(pbeta)  to the end ! ]]
> 
> I've spent many "man weeks", or even "man months" on pbeta() and
> qbeta(), already and have dreamed to get a good student do a
> master's thesis about the problem and potential solutions I've
> looked into in the mean time.
> 
> My current gut feeling is that in some cases, new approximations
> are necessary (i.e. tweaking of current approximations is not
> going to help sufficiently).
> 
> Also not (in the R sources)  tests/p-qbeta-strict-tst.R
> a whole file of "regression tests" about  pbeta() and qbeta()
> {where part of the true values have been computed with my CRAN
> package Rmpfr (for high precision computation) with the
> Rmpfr::pbetaI() function which gives arbitrarily precise pbeta()
> values but only when  (a,b) are integers -- that's the "I" in pbetaI().
> 
> Yes, it's intriguing ... and I'll look into your special
> findings a bit later today.
> 
> 
>   > Should I report this on the bug list?
> 
> Yes, please.  Not all problem of pbeta() / qbeta() are part yet,
> of R's bugzilla data base,  and maybe this will help to draw
> more good applied mathematicians look into it.
> 
> 
> 
> Martin Maechler
> ETH Zurich and R Core team
> (I'd call myself the "dpq-hacker" within R core -- related to
>  my CRAN package 'DPQ')
> 
> 
> > A more general illustration:
> > http://www.math.mcmaster.ca/bolker/misc/qbeta.png
> 
> > ===
> > fun <- function(phi,i=0.01,t=0.001, f=qbeta) {
> > f(t,shape1=i/phi,shape2=(1-i)/phi, lower.tail=FALSE)
> > }
> > ## brute-force beta quantile function
> > Qbeta0 <- function(t,shape1,shape2,lower.tail=FALSE) {
> > fn <- function(x) {pbeta(x,shape1,shape2,lower.tail=lower.tail)-t}
> > uniroot(fn,interval=c(0,1))$root
> > }
> > Qbeta <- Vectorize(Qbeta0,c("t","shape1","shape2"))
> > curve(fun,from=1,to=4)
> > curve(fun(x,f=Qbeta),add=TRUE,col=2)
> 
> > __
> > 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/m

Re: [Rd] unstable corner of parameter space for qbeta?

[J C Nash]

Despite the need to focus on pbeta, I'm still willing to put in some effort.
But I find it really helps to have 2-3 others involved, since the questions back
and forth keep matters moving forward. Volunteers?

Thanks to Martin for detailed comments.

JN


On 2020-03-26 10:34 a.m., Martin Maechler wrote:
>>>>>> J C Nash 
>>>>>> on Thu, 26 Mar 2020 09:29:53 -0400 writes:
> 
> > Given that a number of us are housebound, it might be a good time to 
> try to
> > improve the approximation. It's not an area where I have much 
> expertise, but in
> > looking at the qbeta.c code I see a lot of root-finding, where I do 
> have some
> > background. However, I'm very reluctant to work alone on this, and will 
> ask
> > interested others to email off-list. If there are others, I'll report 
> back.
> 
> Hi John.
> Yes, qbeta() {in its "main branches"}  does zero finding, but
> zero finding of   pbeta(...) - p*   and I tried to explain in my
> last e-mail that the real problem is that already pbeta() is not
> accurate enough in some unstable corners ...
> The order fixing should typically be
> 1) fix pbeta()
> 2) look at qbeta() which now may not even need a fix because its
>problems may have been entirely a consequence of pbeta()'s inaccuracies.
>And if there are cases where the qbeta() problems are not
>only pbeta's "fault", it is still true that the fixes that
>would still be needed crucially depend on the detailed
>working of the function whose zero(s) are sought, i.e.,  pbeta()
> 
> > Ben: Do you have an idea of parameter region where approximation is 
> poor?
> > I think that it would be smart to focus on that to start with.
> 
> 
> 
> Rmpfr  matrix-/vector - products:
> 
> > Martin: On a separate precision matter, did you get my query early in 
> year about double
> > length accumulation of inner products of vectors in Rmpfr? R-help more 
> or
> > less implied that Rmpfr does NOT use extra length. I've been using David
> > Smith's FM Fortran where the DOT_PRODUCT does use double length, but it
> > would be nice to have that in R. My attempts to find "easy" workarounds 
> have
> > not been successful, but I'll admit that other things took precedence.
> 
> Well, the current development version of 'Rmpfr' on R-forge now
> contains facilities to enlarge the precision of the computations
> by a factor 'fPrec' with default 'fPrec = 1';
> notably, instead of  x %*% y   (where the `%*%` cannot have more
> than two arguments) does have a counterpart  matmult(x,y, )
> which allows more arguments, namely 'fPrec', or directly 'precBits';
> and of course there are  crossprod() and tcrossprod() one should
> use when applicable and they also got the  'fPrec' and
> 'precBits' arguments.
> 
> {The %*% etc precision increase still does not work optimally
>  efficiency wise, as it simply increases the precision of all
>  computations by just increasing the precision of x and y (the inputs)}.
> 
> The whole  Matrix and Matrix-vector arithmetic is still
> comparibly slow in Rmpfr .. mostly because I valued human time
> (mine!) much higher than computer time in its implementation.
> That's one reason I would never want to double the precision
> everywhere as it decreases speed even more, and often times
> unnecessarily: doubling the accuracy is basically "worst-case
> scenario" precaution
> 
> Martin
>

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

Re: [Rd] Help useRs to use R's own Time/Date objects more efficiently

[J C Nash]

As with many areas of R usage, my view is that the concern is one
of making it easier to find appropriate information quickly. The
difficulty is that different users have different needs. So if
one wants to know (most of) what is available, the Time Series
Task View is helpful. If one is a novice, it may now be rather
daunting, while I've found, as a long time user of different software,
that I have to dig to find what I need.

In optimization I have tried -- and have had several false starts --
to unify several packages. That could be helpful for time and date
functions.

Another possibility could be to put the "see" and "see also" information
at the TOP of the documentation rather than lower down, and also to
refer to Task Views and possibly other -- eventually R-project --
documentation objects. I happen to favour wiki-like approaches, but
there has not been much movement towards that yet. We R people are
quite strong individualists, but perhaps more team minded thinking
would help. Some of us are getting beyond our best-before date.

However, I support Martin's intent, and hope there will be attempts
in these directions.

Best,

John Nash


On 2020-04-04 5:49 a.m., Martin Maechler wrote:
> This is mostly a RFC  [but *not* about the many extra packages, please..]:
> 
> Noticing to my chagrin  how my students work in a project,
> googling for R code and cut'n'pasting stuff together, accumulating
> this and that package on the way  all just for simple daily time series
> (though with partly missing parts),
> using chron, zoo, lubridate, ...  all for things that are very
> easy in base R *IF* you read help pages and start thinking on
> your own (...), I've noted once more that the above "if" is a
> very strong one, and seems to happen rarely nowadays by typical R users...
> (yes, I stop whining for now).
> 
> In this case, I propose to slightly improve the situation ...
> by adding a few more lines to one help page [[how could that
> help in the age where "google"+"cut'n'paste" has replaced thinking ? .. ]] :
> 
> On R's own ?Dates  help page (and also on ?DateTimeClasses )
> we have pointers, notably
> 
> See Also:
> 
>  ...
>  ...
>  
>  'weekdays' for convenience extraction functions.
> 
> So people must find that and follow the pointer
> (instead of installing one of the dozen helper packages).
> 
> Then on that page, one sees  weekdays(), months() .. julian()
> in the usage ... which don't seem directly helpful for a person
> who needs more.  If that person is diligent and patient (as good useRs are 
> ;-),
> she finds
> 
>Note:
> 
>   Other components such as the day of the month or the year are very
>   easy to compute: just use 'as.POSIXlt' and extract the relevant
>   component.  Alternatively (especially if the components are
>   desired as character strings), use 'strftime'.
> 
> 
> But then, nowadays, the POSIXlt class is not so transparent to the
> non-expert anymore (as it behaves very much like POSIXct, and
> not like a list for good reasons) .. and so 97%  of R users will
> not find this "very easy".
> 
> For this reason, I propose to at add the following to the
> 'Examples:' section of the help file ...
> and I hope that also readers of  R-devel  who have not been
> aware of how to do this nicely,  will now remember (or remember
> where to look?).
> 
> I at least will tell my students in the future to use these or
> write versions of these simple utility functions.
> 
> 
> 
> 
> ## Show how easily you get month, day, year, day (of {month, week, yr}), ... :
> ## (remember to count from 0 (!): mon = 0..11, wday = 0..6,  etc !!)
> 
> ##' Transform (Time-)Date vector  to  convenient data frame :
> dt2df <- function(dt, dName = deparse(substitute(dt)), stringsAsFactors = 
> FALSE) {
> DF <- as.data.frame(unclass(as.POSIXlt( dt )), 
> stringsAsFactors=stringsAsFactors)
> `names<-`(cbind(dt, DF, deparse.level=0L), c(dName, names(DF)))
> }
> dt2df(.leap.seconds)# date+time
> dt2df(Sys.Date() + 0:9) # date
> 
> ##' Even simpler:  Date -> Matrix:
> d2mat <- function(x) simplify2array(unclass(as.POSIXlt(x)))
> d2mat(seq(as.Date("2000-02-02"), by=1, length.out=30)) # has R 1.0.0's 
> release date
> 
> 
> 
> In the distant past / one of the last times I touched on people
> using (base) R's  Date / Time-Date  objects, I had started
> thinking if we should not provide some simple utilities to "base R"
> (not in the 'base' pkg, but rather 'utils') for "extracting" from
> {POSIX(ct), Date} objects ... and we may have discussed that
> within R Core 20 years ago,  and had always thought that this
> shouldn't be hard for useRs themselves to see how to do...
> 
> But then I see that "everybody" uses extension packages instead,
> even in the many situations where there's no gain doing so, 
> but rather increases the dependency-complex

Re: [Rd] New URL redirect checks

[J C Nash]

Does this issue fit in the more general one of centralized vs
partitioned checks? I've suggested before that the CRAN team
seems (and I'll be honest and admit I don't have a good knowledge
of how they work) to favour an all-in-one checking, whereas it
might be helpful to developers and also widen the "CRAN checking"
team to partition checks. Partitioned checks would allow the
particular problems that are raised to be dealt with in a more
focussed action. URLs seem an obvious candidate, since
link checking is used outside of R packages.

I'm sure there are others besides myself who would contribute
to such activities. After all, the partitioned checks could be
contributed packages themselves.

JN



On 2020-09-22 9:50 p.m., Yihui Xie wrote:
> Me too. I have changed some valid URLs in \url{} to \verb{} just to
> avoid these check NOTEs. I do appreciate the check for the validity of
> URLs in packages, especially those dead links (404), but discouraging
> URLs with status code other than 200 (such as 301) feels like
> overdoing the job. After I "hide" links from R CMD check with \verb{}
> , it will be hard to know if these links are still valid in the
> future.
> 
> Regards,
> Yihui
> 
> On Tue, Sep 22, 2020 at 1:17 PM Kevin Wright  wrote:
>>
>> Isn't the whole concept of DOI basically link-shortening/redirecting?
>>
>> For example, this link
>> https://doi.org/10.2134/agronj2016.07.0395
>> redirects to
>> https://acsess.onlinelibrary.wiley.com/doi/abs/10.2134/agronj2016.07.0395
>>
>> As a side note, I got so fed up with CRAN check complaints about (perfectly 
>> valid) re-directs that I refuse to use the \url{} tag anymore.
>>
>> Kevin
>>
>>
>> On Thu, Sep 17, 2020 at 8:32 AM Yihui Xie  wrote:
>>>
>>> I don't have an opinion on the URL shorteners, but how about the
>>> original question? Redirection can be extremely useful in general.
>>> Shortening URLs is only one of its possible applications. FWIW, CRAN
>>> uses (303) redirect itself, e.g.,
>>> https://cran.r-project.org/package=MASS is redirected to
>>> https://cran.r-project.org/web/packages/MASS/index.html Should these
>>> "canonical" CRAN links be disallowed in packages, too? Just as another
>>> example, https://cran.r-project.org/bin/windows/base/release.html is
>>> redirected to the latest Windows installer of R (through the 
>>> tag).
>>>
>>> If the intent of the new URL redirect check is to disallow using URL
>>> shorteners like bit.ly or nyti.ms, that may be fair, but it it is to
>>> disallow using any URLs that are redirected, I think this CRAN policy
>>> may be worth a reconsideration.
>>>
>>> Regards,
>>> Yihui
>>> --
>>> https://yihui.org
>>>
>>>
>>> On Thu, Sep 17, 2020 at 3:26 AM Gábor Csárdi  wrote:

 Right, I am sorry, I did not realize the security aspect here. I guess
 I unconsciously treated CRAN package authors as a trusted source.

 Thanks for the correction and clarification, and to CRAN for
 implementing these checks. :)

 G.

 On Wed, Sep 16, 2020 at 10:50 PM Duncan Murdoch
  wrote:
>
> On 16/09/2020 4:51 p.m., Simon Urbanek wrote:
>> I can't comment for CRAN, but generally, shorteners are considered 
>> security risk so regardless of the 301 handling I think flagging those 
>> is a good idea. Also I think it is particularly bad to use them in 
>> manuals because it hides the target so the user has no idea what hey 
>> will get.
>
> I agree, and we do have \href{}{} in Rd files and similar in other
> formats for giving text of a link different than the URL if the URL is
> inconveniently long.  There's still a bit of a security issue though:
> the built in help browser (at least in MacOS) doesn't show the full URL
> when you hover over the link, as most browsers do.  So one could have
>
> \href{https://disney.org}{https://horrible.web.site}
>
> Duncan Murdoch
>
>
>>
>> Cheers,
>> Simon
>>
>>
>>> On Sep 17, 2020, at 5:35 AM, Gábor Csárdi  
>>> wrote:
>>>
>>> Dear all,
>>>
>>> the new CRAN URL checks flag HTTP 301 redirects. While I understand
>>> the intent, I think this is unfortunate, because several URL shortener
>>> services use 301 redirects, and often a shorter URL is actually better
>>> in a manual page than a longer one that can be several lines long in
>>> the console and also potentially truncated in the PDF manual.
>>>
>>> Some example shorteners that are flagged:
>>>
 db <- tools:::url_db(c("https://nyti.ms";, "https://t.co/mtXLLfYOYE";), 
 "README")
 tools:::check_url_db(db)
>>> URL: https://nyti.ms (moved to https://www.nytimes.com/)
>>> From: README
>>> Status: 200
>>> Message: OK
>>>
>>> URL: https://t.co/mtXLLfYOYE (moved to
>>> https://www.bbc.co.uk/news/blogs-trending-47975564)
>>> From: README
>>> Status: 200
>>> Message: OK
>>>
>>> _

Re: [Rd] URL checks

[J C Nash]

Is this a topic for Google Summer of Code? See
https://github.com/rstats-gsoc/gsoc2021/wiki


On 2021-01-09 12:34 p.m., Dirk Eddelbuettel wrote:
> 
> The idea of 'white lists' to prevent known (and 'tolerated') issues, note,
> warnings, ... from needlessly reappearing is very powerful and general, and
> can go much further than just URL checks.
> 
> I suggested several times in the past that we can look at the format Debian
> uses in its 'lintian' package checker and its override files -- which are
> used across thousands of packages there.  But that went nowhere so I stopped.
> 
> This issue needs a champion or two to implement a prototype as well as a
> potential R Core / CRAN sponsor to adopt it.  But in all those years no smoke
> has come out of any chimneys so ...  ¯\_(ツ)_/¯ is all we get.
> 
> Dirk
>

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

Re: [Rd] URL checks

[J C Nash]

Sorry, Martin, but I've NOT commented on this matter, unless someone has been 
impersonating me.
Someone else?

JN


On 2021-01-11 4:51 a.m., Martin Maechler wrote:
>> Viechtbauer, Wolfgang (SP) 
>> on Fri, 8 Jan 2021 13:50:14 + writes:
> 
> > Instead of a separate file to store such a list, would it be an idea to 
> add versions of the \href{}{} and \url{} markup commands that are skipped by 
> the URL checks?
> > Best,
> > Wolfgang
> 
> I think John Nash and you misunderstood -- or then I
> misunderstood -- the original proposal:
> 
> I've been understanding that there should be a  "central repository" of URL
> exceptions that is maintained by volunteers.
> 
> And rather *not* that package authors should get ways to skip
> URL checking..
> 
> Martin
> 
> 
> >> -Original Message-
> >> From: R-devel [mailto:r-devel-boun...@r-project.org] On Behalf Of 
> Spencer
> >> Graves
> >> Sent: Friday, 08 January, 2021 13:04
> >> To: r-devel@r-project.org
> >> Subject: Re: [Rd] URL checks
> >> 
> >> I also would be pleased to be allowed to provide "a list of known
> >> false-positive/exceptions" to the URL tests.  I've been challenged
> >> multiple times regarding URLs that worked fine when I checked them.  We
> >> should not be required to do a partial lobotomy to pass R CMD check ;-)
> >> 
> >> Spencer Graves
> >> 
> >> On 2021-01-07 09:53, Hugo Gruson wrote:
> >>> 
> >>> I encountered the same issue today with 
> https://astrostatistics.psu.edu/.
> >>> 
> >>> This is a trust chain issue, as explained here:
> >>> https://whatsmychaincert.com/?astrostatistics.psu.edu.
> >>> 
> >>> I've worked for a couple of years on a project to increase HTTPS
> >>> adoption on the web and we noticed that this type of error is very
> >>> common, and that website maintainers are often unresponsive to 
> requests
> >>> to fix this issue.
> >>> 
> >>> Therefore, I totally agree with Kirill that a list of known
> >>> false-positive/exceptions would be a great addition to save time to 
> both
> >>> the CRAN team and package developers.
> >>> 
> >>> Hugo
> >>> 
> >>> On 07/01/2021 15:45, Kirill Müller via R-devel wrote:
>  One other failure mode: SSL certificates trusted by browsers that are
>  not installed on the check machine, e.g. the "GEANT Vereniging"
>  certificate from https://relational.fit.cvut.cz/ .
>  
>  K
>  
>  On 07.01.21 12:14, Kirill Müller via R-devel wrote:
> > Hi
> > 
> > The URL checks in R CMD check test all links in the README and
> > vignettes for broken or redirected links. In many cases this 
> improves
> > documentation, I see problems with this approach which I have
> > detailed below.
> > 
> > I'm writing to this mailing list because I think the change needs to
> > happen in R's check routines. I propose to introduce an "allow-list"
> > for URLs, to reduce the burden on both CRAN and package maintainers.
> > 
> > Comments are greatly appreciated.
> > 
> > Best regards
> > 
> > Kirill
> > 
> > # Problems with the detection of broken/redirected URLs
> > 
> > ## 301 should often be 307, how to change?
> > 
> > Many web sites use a 301 redirection code that probably should be a
> > 307. For example, https://www.oracle.com and https://www.oracle.com/
> > both redirect to https://www.oracle.com/index.html with a 301. I
> > suspect the company still wants oracle.com to be recognized as the
> > primary entry point of their web presence (to reserve the right to
> > move the redirection to a different location later), I haven't
> > checked with their PR department though. If that's true, the 
> redirect
> > probably should be a 307, which should be fixed by their IT
> > department which I haven't contacted yet either.
> > 
> > $ curl -i https://www.oracle.com
> > HTTP/2 301
> > server: AkamaiGHost
> > content-length: 0
> > location: https://www.oracle.com/index.html
> > ...
> > 
> > ## User agent detection
> > 
> > twitter.com responds with a 400 error for requests without a user
> > agent string hinting at an accepted browser.
> > 
> > $ curl -i https://twitter.com/
> > HTTP/2 400
> > ...
> > ...Please switch to a supported browser..
> > 
> > $ curl -s -i https://twitter.com/ -A "Mozilla/5.0 (X11; Ubuntu; 
> Linux
> > x86_64; rv:84.0) Gecko/20100101 Firefox/84.0" | head -n 1
> > HTTP/2 200
> > 
> > # Impact
> > 
> > While the latter problem *could* be fixed by

[Rd] Google Summer of Code for R projects / mentors needed

[J C Nash]

One of the mechanisms by which R has been extended and improved has been through
the efforts of students and mentors in the Google Summer of Code initiatives. 
This
year Toby Hocking (along with others) has continued to lead this effort.

This year, Google has changed the format somewhat so that the projects are 
shorter.
There will likely be more of them.

For success, we have learned that projects need at least 2 mentors -- illness,
life events, world events and holidays can get in the way of monitoring student
work and handling the minor but critical short reports to Google to ensure money
gets to deserving students (and does not get sent otherwise!).

Please consider mentoring and/or proposing a project. See
https://github.com/rstats-gsoc/gsoc2021/wiki/

As an example, my own proposal concerns improving the behaviour and features
of the nls() function. 
https://github.com/rstats-gsoc/gsoc2021/wiki/Improvements-to-nls()

John Nash

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

[Rd] GSoC project "Improvement to nls()"

[J C Nash]

This message is to let R developers know that the project in the Subject is now
a Google Summer of Code project.

Our aim in this project is to find simplifications and corrections to the nls()
code, which has become heavily patched. Moreover, it has some deficiencies in 
that
there is no Marquardt stabilization and it is likely the jacobian (called 
gradient
in R) computations are less than ideal. On the other hand, it has a lot of 
features
and capabilities.

A correction I proposed to avoid the "small residual" issue (when models are 
nearly
perfectly fittable) is now in R-devel. Using a new nls.control parameter one 
can avoid
failure, but the default value of 0 leaves legacy behaviour. We hope to be able 
to use
similar approaches so existing nls() example output is unaltered.

It is likely we will only partially meet our goals:
- to document and possibly simplify the existing code
- to correct some minor issues in documentation or function
- to find a graceful way to incorporate a Marquardt stabilization into the
  Gauss-Newton iteration
- to document, evaluate, and possibly improve the jacobian computation

all within the context that any changes impose minimal nuisance for R workers.

Some of these efforts overlap the nlsr, minpack.lm and likely other packages,
and may suggest improvements there also.

There is a gitlab repository established at 
https://gitlab.com/nashjc/improvenls.
We welcome interest and participation off list except for bugs in the current R
function(s).

John Nash
University of Ottawa

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

[Rd] subset argument in nls() and possibly other functions

[J C Nash]

In mentoring and participating in a Google Summer of Code project "Improvements 
to nls()",
I've not found examples of use of the "subset" argument in the call to nls(). 
Moreover,
in searching through the source code for the various functions related to 
nls(), I can't
seem to find where subset is used, but a simple example, included below, 
indicates it works.
Three approaches all seem to give the same results.

Can someone point to documentation or code so we can make sure we get our 
revised programs
to work properly? The aim is to make them more maintainable and provide 
maintainer documentation,
along with some improved functionality. We seem, for example, to already be 
able to offer
analytic derivatives where they are feasible, and should be able to add 
Marquardt-Levenberg
stabilization as an option.

Note that this "subset" does not seem to be the "subset()" function of R.

John Nash

# CroucherSubset.R -- https://walkingrandomly.com/?p=5254

xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
ydata = 
c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
Cform <- ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)
Cstart<-list(p1=1,p2=0.2)
Cdata<-data.frame(xdata, ydata)
Csubset<-1:8 # just first 8 points

# Original problem - no subset
fit0 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), data=Cdata, 
start=list(p1=1,p2=.2))
summary(fit0)

# via subset argument
fit1 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), data=Cdata, 
start=list(p1=1,p2=.2), subset=Csubset)
summary(fit1)

# via explicit subsetting
Csdata <- Cdata[Csubset, ]
Csdata
fit2 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), data=Csdata, 
start=list(p1=1,p2=.2))
summary(fit2)

# via weights -- seems to give correct observation count if zeros not recognized
wts <- c(rep(1,8), rep(0,2))
fit3 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), data=Cdata, 
weights=wts, start=list(p1=1,p2=.2))
summary(fit3)

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

[Rd] Seeking opinions on possible change to nls() code

[J C Nash]

In our work on a Google Summer of Code project "Improvements to nls()",
the code has proved sufficiently entangled that we have found (so far!)
few straightforward changes that would not break legacy behaviour. One
issue that might be fixable is that nls() returns no result if it
encounters some computational blockage AFTER it has already found a
much better "fit" i.e. set of parameters with smaller sum of squares.
Here is a version of the Tetra example:

time=c( 1,  2,  3,  4,  6 , 8, 10, 12, 16)
conc = c( 0.7, 1.2, 1.4, 1.4, 1.1, 0.8, 0.6, 0.5, 0.3)
NLSdata <- data.frame(time,conc)
NLSstart <-c(lrc1=-2,lrc2=0.25,A1=150,A2=50) # a starting vector (named!)
NLSformula <-conc ~ A1*exp(-exp(lrc1)*time)+A2*exp(-exp(lrc2)*time)
tryit <- try(nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE))
print(tryit)

If you run this, tryit does not give information that the sum of squares
has been reduced from > 6 to < 2, as the trace shows.

Should we propose that this be changed so the returned object gives the
best fit so far, albeit with some form of message or return code to indicate
that this is not necessarily a conventional solution? Our concern is that
some examples might need to be adjusted slightly, or we might simply add
the "try-error" class to the output information in such cases.

Comments are welcome, as this is as much an infrastructure matter as a
computational one.

Best,

John Nash

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

Re: [Rd] Seeking opinions on possible change to nls() code

[J C Nash]

Thanks Martin. I'd missed the intention of that option, but re-reading
it now it is obvious.

FWIW, this problem is quite nasty, and so far I've found no method
that reveals the underlying dangers well. And one of the issues with
nonlinear models is that they reveal how slippery the concept of
inference can be when applied to parameters in such models.

JN


On 2021-08-20 11:35 a.m., Martin Maechler wrote:
>>>>>> J C Nash 
>>>>>> on Fri, 20 Aug 2021 11:06:25 -0400 writes:
> 
> > In our work on a Google Summer of Code project
> > "Improvements to nls()", the code has proved sufficiently
> > entangled that we have found (so far!)  few
> > straightforward changes that would not break legacy
> > behaviour. One issue that might be fixable is that nls()
> > returns no result if it encounters some computational
> > blockage AFTER it has already found a much better "fit"
> > i.e. set of parameters with smaller sum of squares.  Here
> > is a version of the Tetra example:
> 
> time=c( 1,  2,  3,  4,  6 , 8, 10, 12, 16)
> conc = c( 0.7, 1.2, 1.4, 1.4, 1.1, 0.8, 0.6, 0.5, 0.3)
> NLSdata <- data.frame(time,conc)
> NLSstart <-c(lrc1=-2,lrc2=0.25,A1=150,A2=50) # a starting vector (named!)
> NLSformula <-conc ~ A1*exp(-exp(lrc1)*time)+A2*exp(-exp(lrc2)*time)
> tryit <- try(nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE))
> print(tryit)
> 
> > If you run this, tryit does not give information that the
> > sum of squares has been reduced from > 6 to < 2, as
> > the trace shows.
> 
> > Should we propose that this be changed so the returned
> > object gives the best fit so far, albeit with some form of
> > message or return code to indicate that this is not
> > necessarily a conventional solution? Our concern is that
> > some examples might need to be adjusted slightly, or we
> > might simply add the "try-error" class to the output
> > information in such cases.
> 
> > Comments are welcome, as this is as much an infrastructure
> > matter as a computational one.
> 
> Hmm...  many years ago, we had introduced the  'warnOnly=TRUE'
> option to nls()  i.e., nls.control()  exactly for such cases,
> where people would still like to see the solution:
> 
> So,
> 
> --
>> try2 <- nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE, 
>   control = nls.control(warnOnly=TRUE))
> 61215.76(3.56e+03): par = (-2 0.25 150 50)
> 2.175672(2.23e+01): par = (-1.9991 0.3171134 2.618224 -1.366768)
> 1.621050(7.14e+00): par = (-1.960475 -2.620293 2.575261 -0.5559918)
> Warning message:
> In nls(NLSformula, data = NLSdata, start = NLSstart, trace = TRUE,  :
>   singular gradient
> 
>> try2
> Nonlinear regression model
>   model: conc ~ A1 * exp(-exp(lrc1) * time) + A2 * exp(-exp(lrc2) * time)
>data: NLSdata
>lrc1lrc2  A1  A2 
>  -22.89   96.43  156.70 -156.68 
>  residual sum-of-squares: 218483
> 
> Number of iterations till stop: 2 
> Achieved convergence tolerance: 7.138
> Reason stopped: singular gradient
> 
>> coef(try2)
>   lrc1   lrc2 A1 A2 
>  -22.88540   96.42686  156.69547 -156.68461 
> 
> 
>> summary(try2)
> Error in chol2inv(object$m$Rmat()) : 
>   element (3, 3) is zero, so the inverse cannot be computed
>>
> --
> 
> and similar for  vcov(), of course, where the above error
> originates.
> 
> { I think  GSoC (andr other)  students should start by studying and
>   exploring relevant help pages before drawing conclusions
>   ..
>   but yes, I've been born in the last millennium ...
> }
> 
> ;-)
> 
> Have a nice weekend!
> Martin
>

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

Re: [Rd] Bug in optim for specific orders of magnitude

[J C Nash]

Extreme scaling quite often ruins optimization calculations. If you think 
available methods
are capable of doing this, there's a bridge I can sell you in NYC.

I've been trying for some years to develop a good check on scaling so I can 
tell users
who provide functions like this to send (lots of) money and I'll give them the 
best answer
there is (generally no answer at all). Or, more seriously, to inform them that 
they should
not expect results unless they scale. Richard Varga once said some decades ago 
that any
problem was trivially solvable in the right scale, and he was mostly right. 
Scaling is
important.

To see the range of answers from a number of methods, the script below is 
helpful. I had
to remove lbfgsb3c from the mix as it stopped mid-calculation in unrecoverable 
way. Note
that I use my development version of optimx, so some methods might not be 
included in
CRAN offering. Just remove the methods from the ameth and bmeth lists if 
necessary.

Cheers, John Nash

# CErickson221223.R
# optim(c(0,0), function(x) {x[1]*1e-317}, lower=c(-1,-1), upper=c(1,1),
#  method='L-BFGS-B')

tfun <- function(x, xpnt=317){
  if ((length(x)) != 2) {stop("Must have length 2")}
  scl <- 10^(-xpnt)
  val <- x[1]*scl # note that x[2] unused. May be an issue!
  val
}
gtfun <- function(x, xpnt=317){ # gradient
  scl <- 10^(-xpnt)
  gg <- c(scl, 0.0)
  gg
}


xx <- c(0,0)
lo <- c(-1,-1)
up <- c(1,1)
print(tfun(xx))
library(optimx)
ameth <- c("BFGS", "CG", "Nelder-Mead", "L-BFGS-B", "nlm", "nlminb",
 "Rcgmin", "Rtnmin", "Rvmmin", "spg", "ucminf", "newuoa", "bobyqa",
 "nmkb", "hjkb",  "hjn", "lbfgs", "subplex", "ncg", "nvm", "mla",
 "slsqp", "anms")

bmeth <- c("L-BFGS-B", "nlminb", "Rcgmin", "Rtnmin", "nvm",
"bobyqa", "nmkb", "hjkb", "hjn", "ncg", "slsqp")

tstu <- opm(x<-c(0,0), fn=tfun, gr=gtfun, method=ameth, control=list(trace=0))
summary(tstu, order=value)

tstb <- opm(x<-c(0,0), fn=tfun, gr=gtfun, method=bmeth, lower=lo, upper=up,
control=list(trace=0))
summary(tstb, order=value)


On 2022-12-23 13:41, Rui Barradas wrote:

Às 17:30 de 23/12/2022, Collin Erickson escreveu:

Hello,

I've come across what seems to be a bug in optim that has become a nuisance
for me.

To recreate the bug, run:

optim(c(0,0), function(x) {x[1]*1e-317}, lower=c(-1,-1), upper=c(1,1),
method='L-BFGS-B')

The error message says:

Error in optim(c(0, 0), function(x) { :
   non-finite value supplied by optim

What makes this particularly treacherous is that this error only occurs for
specific powers. By running the following code you will find that the error
only occurs when the power is between -309 and -320; above and below that
work fine.

p <- 1:1000
giveserror <- rep(NA, length(p))
for (i in seq_along(p)) {
   tryout <- try({
 optim(c(0,0), function(x) {x[1]*10^-p[i]}, lower=c(-1,-1),
upper=c(1,1), method='L-BFGS-B')
   })
   giveserror[i] <- inherits(tryout, "try-error")
}
p[giveserror]

Obviously my function is much more complex than this and usually doesn't
fail, but this reprex demonstrates that this is a problem. To avoid the
error I may multiply by a factor or take the log, but it seems like a
legitimate bug that should be fixed.

I tried to look inside of optim to track down the error, but the error lies
within the external C code:

.External2(C_optim, par, fn1, gr1, method, con, lower,
 upper)

For reference, I am running R 4.2.2, but was also able to recreate this bug
on another device running R 4.1.2 and another running 4.0.3.

Thanks,
Collin Erickson

[[alternative HTML version deleted]]

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


Hello,

See if this R-Help thread [1] earlier this year is relevant.
In particular, the post by R Core team member Luke Tierney [2], that answers 
much better than what I could.

The very small numbers in your question seem to have hit a limit and this limit 
is not R related.


[1] https://stat.ethz.ch/pipermail/r-help/2022-February/473840.html
[2] https://stat.ethz.ch/pipermail/r-help/2022-February/473844.html


Hope this helps,

Rui Barradas

__
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] Potential bug in fitted.nls

[J C Nash]

FWIW, nlsr::nlxb() gives same answers.

JN

On 2023-01-25 09:59, Dave Armstrong wrote:

Dear Colleagues,

I recently answered [this question]() on StackOverflow that identified
what seems to be unusual behaviour with `stats:::nls.fitted()`. In
particular, a null model returns a single fitted value rather than a
vector of the same fitted value of `length(y)`.  The documentation
doesn’t make it seem like this is the intended behaviour, so I’m not
sure if it’s a bug, a “Wishlist” item or something that is working
as intended even though it seems unusual to me.  I looked through the
bug reporting page on the R project website and it suggested contacting
the R-devel list in cases where the behaviour is not obviously a bug to
see whether others find the behaviour equally unusual and I should
submit a Wishlist item through Bugzilla.

Below is a reprex that shows how the fitted values of a model with just
a single parameter is length 1, but if I multiply that constant by a
vector of ones, then the fitted values are of `length(y)`.  Is this
something that should be reported?

``` r
dat <-
data.frame(y=c(80,251,304,482,401,141,242,221,304,243,544,669,638),
ones = rep(1, 13))
mNull1 <- nls(y ~ a, data=dat, start=c(a=mean(dat$y)))
fitted(mNull1)
#> [1] 347.6923
#> attr(,"label")
#> [1] "Fitted values"

mNull2 <- nls(y ~ a*ones, data=dat, start=c(a=mean(dat$y)))
fitted(mNull2)
#>  [1] 347.6923 347.6923 347.6923 347.6923 347.6923 347.6923 347.6923
347.6923
#>  [9] 347.6923 347.6923 347.6923 347.6923 347.6923
#> attr(,"label")
#> [1] "Fitted values"
```

Created on 2023-01-25 by the [reprex
package](https://reprex.tidyverse.org) (v2.0.1)
[[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] Potential bug in fitted.nls

[J C Nash]

nls() actually uses different modeling formulas depending on the 'algorithm', 
and
there is, in my view as a long time nonlinear modeling person, an unfortunate
structural issue that likely cannot be resolved simply. This is because for 
nonlinear
modeling programs we really should be using explicit model statements
e.g., a linear model should be y ~ a * x + b where x is the (independent)
variable and a and b are parameters. But we put in y ~ x as per lm().
Partially linear approaches and indexed parameter models add complexity and
inconsistency. We're pushing structures beyond their design.

This is one of the topics in a paper currently undergoing
final edit on "Improving nls()" that came out of a Google Summer of Code 
project.
The desired improvements in nls() were mostly frustrated by entanglements in 
the code,
but they have led to a lot of tweaks to package nlsr. Perhaps someone more 
facile
with the intricacies of R internals can succeed. For the paper and nlsr,
all the bits should get sent to CRAN and elsewhere in the next month or so,
(co-author is newly started on a Ph D) but if anyone is anxious to try, they 
can email
me. The nlsr code has been stable for several months, but some documentation 
still
being considered.

Sorting out how to deal with the model expression for nls() and related tools
is a worthwhile goal, but not one that can be settled here. It could make a good
review project for a senior undergrad or master's level, and I'd be happy to
join the discussion.

Cheers, JN


On 2023-01-26 12:55, Bill Dunlap wrote:

Doesn't nls() expect that the lengths of vectors on both sides of the
formula match (if both are supplied)?  Perhaps it should check for that.

-Bill

On Thu, Jan 26, 2023 at 12:17 AM Dave Armstrong  wrote:


Dear Colleagues,

I recently answered [this question]() on StackOverflow that identified
what seems to be unusual behaviour with `stats:::nls.fitted()`. In
particular, a null model returns a single fitted value rather than a
vector of the same fitted value of `length(y)`.  The documentation
doesn’t make it seem like this is the intended behaviour, so I’m not
sure if it’s a bug, a “Wishlist” item or something that is working
as intended even though it seems unusual to me.  I looked through the
bug reporting page on the R project website and it suggested contacting
the R-devel list in cases where the behaviour is not obviously a bug to
see whether others find the behaviour equally unusual and I should
submit a Wishlist item through Bugzilla.

Below is a reprex that shows how the fitted values of a model with just
a single parameter is length 1, but if I multiply that constant by a
vector of ones, then the fitted values are of `length(y)`.  Is this
something that should be reported?

``` r
dat <-
data.frame(y=c(80,251,304,482,401,141,242,221,304,243,544,669,638),
ones = rep(1, 13))
mNull1 <- nls(y ~ a, data=dat, start=c(a=mean(dat$y)))
fitted(mNull1)
#> [1] 347.6923
#> attr(,"label")
#> [1] "Fitted values"

mNull2 <- nls(y ~ a*ones, data=dat, start=c(a=mean(dat$y)))
fitted(mNull2)
#>  [1] 347.6923 347.6923 347.6923 347.6923 347.6923 347.6923 347.6923
347.6923
#>  [9] 347.6923 347.6923 347.6923 347.6923 347.6923
#> attr(,"label")
#> [1] "Fitted values"
```

Created on 2023-01-25 by the [reprex
package](https://reprex.tidyverse.org) (v2.0.1)
 [[alternative HTML version deleted]]

__
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


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

Re: [Rd] uniroot violates bounds?

[J C Nash]

I wrote first cut at unirootR for Martin M and he revised and put in
Rmpfr.

The following extends Ben's example, but adds the unirootR with trace
output.

c1 <- 4469.822
c2 <- 572.3413
f <- function(x) { c1/x - c2/(1-x) }; uniroot(f, c(1e-6, 1))
uniroot(f, c(1e-6, 1))
library(Rmpfr)
unirootR(f, c(1e-6, 1), extendInt="no", trace=1)

This gives more detail on the iterations, and it looks like the Inf is the
issue. But certainly we could do more to avoid "gotchas" like this. If
someone is interested in some back and forth, I'd be happy to give it a
try, but I think progress would be better with more than one contributor.

Best,

John Nash

On 2023-02-18 12:28, Ben Bolker wrote:

c1 <- 4469.822
c2 <- 572.3413
f <- function(x) { c1/x - c2/(1-x) }; uniroot(f, c(1e-6, 1))
uniroot(f, c(1e-6, 1))


    provides a root at -6.00e-05, which is outside of the specified bounds.  The default value of the "extendInt" 
argument to uniroot() is "no", as far as I can see ...


$root
[1] -6.003516e-05

$f.root
[1] -74453981

$iter
[1] 1

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05


   I suspect this fails because f(1) (value at the upper bound) is infinite, although setting interval to c(0.01, 1) 
does work/give a sensible answer ...  (works for a lower bound of 1e-4, fails for 1e-5 ...)


   Setting the upper bound < 1 appears to avoid the problem.

  For what it's worth, the result has an "init.it" component, but the only thing the documentation says about it is " 
component ‘init.it’ was added in R 3.1.0".


   And, I think (?) that the 'trace' argument only produces any output if the 
'extendInt' option is enabled?

   Inspired by 
https://stackoverflow.com/questions/75494696/solving-a-system-of-non-linear-equations-with-only-one-unknown/75494955#75494955


   cheers
    Ben Bolker

__
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] Query: Could documentation include modernized references?

[J C Nash]

A tangential email discussion with Simon U. has highlighted a long-standing
matter that some tools in the base R distribution are outdated, but that
so many examples and other tools may use them that they cannot be deprecated.

The examples that I am most familiar with concern optimization and nonlinear
least squares, but other workers will surely be able to suggest cases elsewhere.
I was the source (in Pascal) of Nelder-Mead, BFGS and CG algorithms in optim().
BFGS is still mostly competitive, and Nelder-Mead is useful for initial 
exploration
of an optimization problem, but CG was never very good, right from the mid-1970s
well before it was interfaced to R. By contrast Rcgmin works rather well
considering how similar it is in nature to CG. Yet I continue to see use and
even recommendations of these tools in inappropriate circumstances.

Given that it would break too many other packages and examples to drop the
existing tools, should we at least add short notes in the man (.Rd) pages?
I'm thinking of something like

   optim() has methods that are dated. Users are urged to consider suggestions
   from ...

and point to references and/or an appropriate Task View, which could, of course,
be in the references.

I have no idea what steps are needed to make such edits to the man pages. Would
R-core need to be directly involved, or could one or two trusted R developers
be given privileges to seek advice on and implement such modest documentation
additions?  FWIW, I'm willing to participate in such an effort, which I believe
would help users to use appropriate and up-to-date tools.

John Nash

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

Re: [Rd] Query: Could documentation include modernized references?

[J C Nash]

Thanks Martin.

Following Duncan's advice as well as some textual input, I have put a proposed 
Rd file for
optim on a fork of the R code at
https://github.com/nashjc/r/blob/master/src/library/stats/man/optim.Rd

This has the diffs given below from the github master. The suggested changes
primarily point to a Task View, which I believe is a sensible approach.

I'll admit to being rather clumsy with git, and will be happy to receive advice
on how to proceed if more work needed on my part.

Cheers,

John Nash




--- optim.Rd2022-03-24 19:02:04.0 -0400
+++ optim.Rd.20230324.txt   2023-03-29 09:23:28.373457291 -0400
@@ -15,6 +15,9 @@
   General-purpose optimization based on Nelder--Mead, quasi-Newton and
   conjugate-gradient algorithms. It includes an option for
   box-constrained optimization and simulated annealing.
+  These methods are quite old and better ones are known for many
+  problems.  See the Optimization and Mathematical Programming task
+  view (Schwendinger and Borchers, 2023) for a survey.
 }
 \usage{
 optim(par, fn, gr = NULL, \dots,
@@ -67,6 +70,8 @@
   Beale--Sorenson updates).  Conjugate gradient methods will generally
   be more fragile than the BFGS method, but as they do not store a
   matrix they may be successful in much larger optimization problems.
+  The \code{"CG"} method has known improvements that are discussed in
+  Schwendinger and Borchers (2023)."

   Method \code{"L-BFGS-B"} is that of Byrd \emph{et. al.} (1995) which
   allows \emph{box constraints}, that is each variable can be given a lower
@@ -230,8 +235,10 @@
 \source{
   The code for methods \code{"Nelder-Mead"}, \code{"BFGS"} and
   \code{"CG"} was based originally on Pascal code in Nash (1990) that was
-  translated by \code{p2c} and then hand-optimized.  Dr Nash has agreed
-  that the code can be made freely available.
+  translated by \code{p2c} and then hand-optimized.  Dr Nash has agreed
+  that the code can be made freely available, but recommends that the more
+  reliable \code{optimx::Rcgmin()} function should be used instead of
+  method \code{"CG"}.

   The code for method \code{"L-BFGS-B"} is based on Fortran code by Zhu,
   Byrd, Lu-Chen and Nocedal obtained from Netlib (file
@@ -269,6 +276,10 @@
   Nocedal, J. and Wright, S. J. (1999).
   \emph{Numerical Optimization}.
   Springer.
+   
+  Florian Schwendinger, Hans W. Borchers (2023). \emph{CRAN Task View:
+  Optimization and Mathematical Programming.} Version 2023-02-16.
+  URL https://CRAN.R-project.org/view=Optimization.
 }

 \seealso{




On 2023-03-31 09:31, Martin Maechler wrote:


Thanks a lot, Duncan, for this (as usual from you) very precise
and helpful information / explanations.

I am "happy"/willing to get involved a bit here, as I do want to
spend some time re-reading about current state of (some, notably
optim-related) optimizers.

(But I will be mostly offline for the next 60 hours or so.)


Martin

--
Martin Maechler
ETH Zurich  and  R Core team


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

Re: [Rd] Concerns with SVD

[J C Nash]

Better check your definitions of SVD -- there are several forms, but all I
am aware of (and I wrote a couple of the codes in the early 1970s for the SVD)
have positive singular values.

JN


On 2023-07-16 02:01, Durga Prasad G me14d059 wrote:

Respected Development Team,

This is Durga Prasad reaching out to you regarding an issue/concern related
to Singular Value Decomposition SVD in R software package. I am attaching a
detailed attachment with this letter which depicts real issues with SVD in
R.

To reach the concern the expressions for the exponential of a matrix using
SVD and
projection tensors are obtained from series expansion. However, numerical
inconsistency is observed between the exponential of matrix obtained using
the function(svd()) used in R software.

However, it is observed that most of the researchers fraternity is engaged
in utilising R software for their research purposes and to the extent of my
understanding such an error in SVD in R software might raise the concern
about authenticity of the simulation results produced and published by
researchers across the globe.

Further, I am very sure that the R software development team is well versed
with the happening and they have any specific and resilient reasons for
doing so. I would request you kindly, to guide me through the concern.

Thank you very much.


__
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] feature request: optim() iteration of functions that return multiple values

[J C Nash]

But why time methods that the author (me!) has been telling the community for
years have updates? Especially as optimx::optimr() uses same syntax as optim()
and gives access to a number of solvers, both production and didactic. This set
of solvers is being improved or added to regularly, with a major renewal almost
complete (for the adventurous, code on https://github.com/nashjc/optimx).

Note also that the default Nelder-Mead is good for exploring function surface 
and
is quite robust at getting quickly into the region of a minimum, but can be 
quite
poor in "finishing" the process. Tools have different strengths and weaknesses.
optim() was more or less state of the art a couple of decades ago, but there are
other choices now.

JN

On 2023-08-08 05:14, Sami Tuomivaara wrote:

Thank you all very much for the suggestions, after testing, each of them would 
be a viable solution in certain contexts.  Code for benchmarking:

# preliminaries
install.packages("microbenchmark")
library(microbenchmark)


data <- new.env()
data$ans2 <- 0
data$ans3 <- 0
data$i <- 0
data$fun.value <- numeric(1000)

# define functions

rosenbrock_env <- function(x, data)
{
   x1 <- x[1]
   x2 <- x[2]
   ans <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
   ans2 <- ans^2
   ans3 <- sqrt(abs(ans))
   data$i <- data$i + 1
   data$fun.value[data$i] <- ans
   ans
}


rosenbrock_env2 <- function(x, data)
{
   x1 <- x[1]
   x2 <- x[2]
   ans <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
   ans2 <- ans^2
   ans3 <- sqrt(abs(ans))
   data$ans2 <- ans2
   data$ans3 <- ans3
   ans
}

rosenbrock_attr <- function(x)
{
   x1 <- x[1]
   x2 <- x[2]
   ans <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
   ans2 <- ans^2
   ans3 <- sqrt(abs(ans))
   attr(ans, "ans2") <- ans2
   attr(ans, "ans3") <- ans3
   ans
}


rosenbrock_extra <- function(x, extraInfo = FALSE)
{
   x1 <- x[1]
   x2 <- x[2]
   ans <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
   ans2 <- ans^2
   ans3 <- sqrt(abs(ans))
   if (extraInfo) list(ans = ans, ans2 = ans2, ans3 = ans3)
   else ans
}


rosenbrock_all <- function(x)
{
   x1 <- x[1]
   x2 <- x[2]
   ans <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
   ans2 <- ans^2
   ans3 <- sqrt(abs(ans))
   list(ans = ans, ans2 = ans2, ans3 = ans3)
}

returnFirst <- function(fun) function(...) do.call(fun,list(...))[[1]]
rosenbrock_all2 <- returnFirst(rosenbrock_all)


# benchmark all functions
set.seed <- 100

microbenchmark(env = optim(c(-1,2), rosenbrock_env, data = data),
env2 = optim(c(-1,2), rosenbrock_env2, data = data),
attr = optim(c(-1,2), rosenbrock_attr),
extra = optim(c(-1,2), rosenbrock_extra, extraInfo = FALSE),
all2 = optim(c(-1,2), rosenbrock_all2),
times = 100)


# correct parameters and return values?
env <- optim(c(-1,2), rosenbrock_env, data = data)
env2 <- optim(c(-1,2), rosenbrock_env2, data = data)
attr <- optim(c(-1,2), rosenbrock_attr)
extra <- optim(c(-1,2), rosenbrock_extra, extraInfo = FALSE)
all2 <- optim(c(-1,2), rosenbrock_all2)

# correct return values with optimized parameters?
env. <- rosenbrock_env(env$par, data)
env2. <- rosenbrock_env(env2$par, data)
attr. <- rosenbrock_attr(attr$par)
extra. <- rosenbrock_extra(extra$par, extraInfo = FALSE)
all2. <- rosenbrock_all2(all2$par)

# functions that return more than one value
all. <- rosenbrock_all(all2$par)
extra2. <- rosenbrock_extra(extra$par, extraInfo = TRUE)

# environment values correct?
data$ans2
data$ans3
data$i
data$fun.value


microbenchmarking results:

Unit: microseconds
   expr minlq  meanmedian uq   max neval
env 644.102 3919.6010 9598.3971 7950.0005 15582.8515 42210.900   100
   env2 337.001  351.5510  479.2900  391.7505   460.3520  6900.800   100
   attr 350.201  367.3010  502.0319  409.7510   483.6505  6772.800   100
  extra 276.800  287.2010  402.4231  302.6510   371.5015  6457.201   100
   all2 630.801  646.9015  785.9880  678.0010   808.9510  6411.102   100

rosenbrock_env and _env2 functions differ in that _env accesses vectors in the 
defined environment by indexing, whereas _env2 doesn't (hope I interpreted this 
right?).  This appears to be expensive operation, but allows saving values 
during the steps of the optim iteration, rather than just at convergence.  
Overall, _extra has consistently lowest median execution time!

My earlier workaround was to write two separate functions, one of which returns 
extra values; all suggested approaches simplify that approach considerably.  I 
am also now more educated about attributes and environments that I did not know 
how to utilize before and that proved to be very useful concepts.  Again, thank 
you everyone for your input!


[[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/listin

Re: [Rd] ADA Compliance

[J C Nash]

Slightly tangential, but about two decades ago I was researching
how multimedia databases might be reasonably structured. To have a
concrete test case, I built a database of English Country (Playford)
dances, which I called Playford's Progeny. (Ben B. will be aware of
this, too.) This proved rather popular, but around 2010 the busybody
brigade at uOttawa sent me a demand to prove that the website satisfied
(name your jurisdiction, I think mine was Ontario provincial something)
accessibility requirements.

I figured my time to do this was worth $2-3K and simply went out and
bought service for about $100. It's now hosted on ottawaenglishdance.org.
Interestingly the main contributor to my site at the time was blind.
Go figure.

The point I'm getting at is that it may make people feel good to
legislate about accessibility, but my guess is the old adage of
catching more flies with honey than vinegar is illustrated here to a
horrifying degree. I'm afraid I've no practical advice on how to
satisfy the "rules".

Best of luck getting things available for as many folk as possible,
no matter their particular disabilities. It's something I support,
just not a lot of rules.

JN


On 2024-01-15 07:10, peter dalgaard wrote:

Yes,

Jonathon Godfrey, who wrote the r-devel/2022-December mail (and is himself 
blind), would be my standard go-to guy in matters relating to visual 
impairment, screen readers and all that.

Peter D.


On 13 Jan 2024, at 00:14 , Ben Bolker  wrote:

I would be very surprised if anyone had written up a VPAT 
 for R.

  It won't help you with the bureaucratic requirements, but R is in fact very 
accessible to visually impaired users: e.g. see

https://community.rstudio.com/t/accessibility-of-r-rstudio-compared-to-excel-for-student-that-is-legally-blind/103849/3

 From https://github.com/ajrgodfrey/BrailleR


R is perhaps the most blind-friendly statistical software option because all 
scripts can be written in plain text, using the text editor a user prefers, and 
all output can be saved in a wide range of file formats. The advent of R 
markdown and other reproducible research techniques can offer the blind user a 
degree of efficiency that is not offered in many other statistical software 
options. In addition, the processed Rmd files are usually HTML which are the 
best supported files in terms of screen reader development.


  (And there is continued attention to making sure R stays accessible in this 
way: https://stat.ethz.ch/pipermail/r-devel/2022-December/082180.html; 
https://stat.ethz.ch/pipermail/r-devel/2023-February/082313.html)

  R is also easy to use without a mouse, which should improve accessibility for 
users with neuromuscular conditions.

   cheers
Ben Bolker




On 2024-01-12 2:50 p.m., Hunter, Zayne via R-devel wrote:

Hello,
I am working with Ball State University to obtain a license of R. As part of 
our requirements for obtaining new software, we must review the VPAT for ADA 
compliance. Can you provide this information for me?
Thanks,
Zayne Hunter
Technology Advisor & Vendor Relations Manager
Ball State University
zayne.hun...@bsu.edu
(765)285-7853
[[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




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

Re: [Rd] Advice debugging M1Mac check errors

[J C Nash]

Slightly tangential: I had some woes with some vignettes in my
optimx and nlsr packages (actually in examples comparing to OTHER
packages) because the M? processors don't have 80 bit registers of
the old IEEE 754 arithmetic, so some existing "tolerances" are too
small when looking to see if is small enough to "converge", and one
gets "did not converge" type errors. There are workarounds,
but the discussion is beyond this post. However, worth awareness that
the code may be mostly correct except for appropriate tests of
smallness for these processors.

JN




On 2024-02-04 11:51, Dirk Eddelbuettel wrote:


On 4 February 2024 at 20:41, Holger Hoefling wrote:
| I wanted to ask if people have good advice on how to debug M1Mac package
| check errors when you don´t have a Mac? Is a cloud machine the best option
| or is there something else?

a) Use the 'mac builder' CRAN offers:
https://mac.r-project.org/macbuilder/submit.html

b) Use the newly added M1 runners at GitHub Actions,

https://github.blog/changelog/2024-01-30-github-actions-introducing-the-new-m1-macos-runner-available-to-open-source/

Option a) is pretty good as the machine is set up for CRAN and builds
fast. Option b) gives you more control should you need it.

Dirk



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

Re: [Rd] Advice debugging M1Mac check errors

[J C Nash]

80 bit registers (I don't have my original docs with me here in Victoria)
seem to have been part of the 1985 standard to which I was one of the 31 named
contributors. See

https://stackoverflow.com/questions/612507/what-are-the-applications-benefits-of-an-80-bit-extended-precision-data-type

or the Wikipedia item on IEEE 754.

It appears to have been omitted from 2008 and 2020 versions, but is still (I 
believe) part of
many processors. It's an internal precision for handling multiplications and 
accumulation,
and not one of the storage modes.

Most of the time this makes very little difference in results for R, since it 
is only some
operations where the extended precision gets activated. If we store quantities, 
we get the
regular precision. Thus very few situations using the M? chips give 
differences, but when
they do, it is a nuisance.

There is plenty of scope for debating the pros and cons of extended precision 
internally.
Not having it likely contributes to speed / bang for the buck of the M? chips. 
But we do
now have occasional differences in outcomes which will lead to confusion and 
extra work.

JN




On 2024-02-04 15:26, Duncan Murdoch wrote:

Hi John.

I don't think the 80 bit format was part of IEEE 754; I think it was an Intel invention for the 8087 chip (which I 
believe preceded that standard), and didn't make it into the standard.


The standard does talk about 64 bit and 128 bit floating point formats, but not 
80 bit.

Duncan Murdoch

On 04/02/2024 4:47 p.m., J C Nash wrote:

Slightly tangential: I had some woes with some vignettes in my
optimx and nlsr packages (actually in examples comparing to OTHER
packages) because the M? processors don't have 80 bit registers of
the old IEEE 754 arithmetic, so some existing "tolerances" are too
small when looking to see if is small enough to "converge", and one
gets "did not converge" type errors. There are workarounds,
but the discussion is beyond this post. However, worth awareness that
the code may be mostly correct except for appropriate tests of
smallness for these processors.

JN




On 2024-02-04 11:51, Dirk Eddelbuettel wrote:


On 4 February 2024 at 20:41, Holger Hoefling wrote:
| I wanted to ask if people have good advice on how to debug M1Mac package
| check errors when you don´t have a Mac? Is a cloud machine the best option
| or is there something else?

a) Use the 'mac builder' CRAN offers:
 https://mac.r-project.org/macbuilder/submit.html

b) Use the newly added M1 runners at GitHub Actions,
 
https://github.blog/changelog/2024-01-30-github-actions-introducing-the-new-m1-macos-runner-available-to-open-source/


Option a) is pretty good as the machine is set up for CRAN and builds
fast. Option b) gives you more control should you need it.

Dirk



__
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] Advice debugging M1Mac check errors

[J C Nash]

Simon's comments add another viewpoint to mine. My own knowledge of the
impact of "disable-long-double" does not include an understanding of
exactly what effect this has. One needs to spend a lot of time and effort
with excruciating details. Fortunately, we can usually get away with
64 bit FP arithmetic for almost all applications. I suspect applications
that need really long precision are likely best handled with special
hardware.

JN


On 2024-02-04 16:46, Simon Urbanek wrote:




On Feb 5, 2024, at 12:26 PM, Duncan Murdoch  wrote:

Hi John.

I don't think the 80 bit format was part of IEEE 754; I think it was an Intel 
invention for the 8087 chip (which I believe preceded that standard), and 
didn't make it into the standard.

The standard does talk about 64 bit and 128 bit floating point formats, but not 
80 bit.



Yes, the 80 bit was Intel-specific (motivated by internal operations, not as external 
format), but as it used to be most popular architecture, people didn't quite realize that 
tests relying on Intel results will be Intel-specific (PowerPC Macs had 128-bit floating 
point, but they were not popular enough to cause trouble in the same way). The IEEE 
standard allows "extended precision" formats, but doesn't prescribe their 
format or precision - and they are optional. Arm64 CPUs only support 64-bit double 
precision in hardware (true both on macOS and Windows), so only what is in the basic 
standard. There are 128-bit floating point solutions in software, but, obviously, they 
are a lot slower (several orders of magnitude). Apple has been asking for priorities in 
the scientific community and 128-bit floating number support was not something high on 
people's priority list. It is far from trivial, because there is a long list of 
operations (all variations of the math functions) so I wouldn't expect this to change 
anytime soon - in fact once Microsoft's glacial move is done we'll be likely seeing only 
64-bit everywhere.

That said even if you don't have a arm64 CPU, you can build R with 
--disable-long-double to get closer to the arm64 results if that is your worry.

Cheers,
Simon




On 04/02/2024 4:47 p.m., J C Nash wrote:

Slightly tangential: I had some woes with some vignettes in my
optimx and nlsr packages (actually in examples comparing to OTHER
packages) because the M? processors don't have 80 bit registers of
the old IEEE 754 arithmetic, so some existing "tolerances" are too
small when looking to see if is small enough to "converge", and one
gets "did not converge" type errors. There are workarounds,
but the discussion is beyond this post. However, worth awareness that
the code may be mostly correct except for appropriate tests of
smallness for these processors.
JN
On 2024-02-04 11:51, Dirk Eddelbuettel wrote:


On 4 February 2024 at 20:41, Holger Hoefling wrote:
| I wanted to ask if people have good advice on how to debug M1Mac package
| check errors when you don´t have a Mac? Is a cloud machine the best option
| or is there something else?

a) Use the 'mac builder' CRAN offers:
 https://mac.r-project.org/macbuilder/submit.html

b) Use the newly added M1 runners at GitHub Actions,
 
https://github.blog/changelog/2024-01-30-github-actions-introducing-the-new-m1-macos-runner-available-to-open-source/

Option a) is pretty good as the machine is set up for CRAN and builds
fast. Option b) gives you more control should you need it.

Dirk


__
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

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

Re: [Rd] Advice debugging M1Mac check errors

[J C Nash]

M1mac numerical issues should be rare, but when they do pop up they can be 
disconcerting.

The following little script reveals what happens with no extended precision. A 
few months
ago I built this into a "package" and used 
https://mac.r-project.org/macbuilder/submit.html
to run it, getting the indicated result of 0 for (sum(vv1) - 1.0e0), with 
non-zero on my
Ryzen 7 laptop.

JN

# FPExtendedTest.R   J C Nash
loopsum <- function(vec){
   n <- length(vec)
   vsum<-0.0
   for (i in 1:n) { vsum <- vsum + vec[i]}
   vsum
}
small<-.Machine$double.eps/4 # 1/4 of the machine precision
vsmall <- rep(small, 1e4) # a long vector of small numbers
vv1 <- c(1.0, vsmall) # 1 at the front of this vector
vv2 <- c(vsmall, 1.0) # 1 at the end
(sum(vv1) - 1.0e0) # Should be > 0 for extended precision, 0 otherwise
(sum(vv2) - 1.0e0) # Should be > 0
(loopsum(vv1) - 1.0e0) # should be zero
(loopsum(vv2) - 1.0e0) # should be greater than zero



On 2024-02-06 08:06, Prof Brian Ripley via R-devel wrote:



We were left to guess, but I doubt this has to do with the lack of 'extended precision' nor long doubles longer than 
doubles on arm64 macOS.  And issues with that are rather rare (much rarer than numerical issues for non-reference x86_64 
BLAS/LAPACKs).  Of the 20,300 CRAN packages just 18 have M1mac-specific errors, none obviously from numerical 
inaccuracy.  A quick look back suggests we get about 20 a year with M1mac numerical issues, about half of which were 
mirrored on the x86_64 'noLD' checks.




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

Re: [Rd] [External] Non-API updates

[J C Nash]

On 2024-06-25 12:25, Josiah Parry wrote:

"Stop using them" is pithy advice but far easier said than done!

With respect to NOTES and WARN on CRAN, these do not result in any
package maintainer notifications. The only notification that the developers
receive is the threatening one that states that the packages will be
removed
from CRAN with a very short timeline.

Is there a link you can provide regarding the "Moving into C compliance?" It
cannot be found at https://cran.r-project.org/doc/manuals/R-exts.html.

Given that there hasn't been a stable release with these changes, how are
we
to know these changes have been stabilized and begin working towards
remedying them? These WARNs come from R-devel which is, as its name
suggests,
a development version which is vacillating daily.




On Tue, Jun 25, 2024 at 12:10 PM  wrote:


On Tue, 25 Jun 2024, Josiah Parry wrote:


Hey folks,

I'm sure many of you all woke to the same message I did: "Please correct
before 2024-07-09 to safely retain your package on CRAN" caused by

Non-API

changes to CRAN.

This is quite unexpected as Luke Tierney's June 6th email writes

(emphasis

mine):

"An *experimental* *effort* is underway to add annotations to the WRE..."

"*Once things have gelled a bit more *I hope to turn this into a blog
post that will include some examples of moving non-API entry point
uses into compliance."

Since then there has not been any indication of stabilization of the
Non-API changes nor has there been a blog post outlining how to migrate.

As

things have been coming and going from the Non-API changes for quite some
time now, we (the royal we, here) have been waiting for an official
announcement from CRAN on the stabilizing changes.


I posted an update to this list a few days ago. If you missed it you
can find it in the archive.


*Can we extend this very short notice to handle the Non-API changes

before

removal from CRAN? *


Timing decisions are up to CRAN.


In the case of the 3 packages I have to fix within 2 weeks, these are all
using Rust. These changes require upstream changes to the extendr

library.

There are other packages that are also affected here. Making these

changes

is a delicate act and requires care and focus. All of the extendr
developers work full time and cannot make addressing these changes their
only priority for the next 2 weeks.


Using non-API entry points is a choice that comes with risks. The ones
leading to WARNINGs for your packages (PRSEEN and SYMVALUE)have been
receiving NOTEs in check results for several weeks. Using
tools:::checkPkgAPI you can see that your packages are referencing a
lot of non-API entry points. Some of these may be added to the API,
but most will not. This may be a good time to look into that.

To minimize disruption we have been adding entry points to the API as
long as it is safe to do so, in some cases against our better
judgment. But ones that are unsafe to use will not be
added. Eventually their declarations will be removed from public
header files and they will be hidden when that is possible. Packages
that have chosen to use these non-API entry points will have to adapt
if they want to pass R CMD check. For now, we will try to first have
use of these entry points result in NOTEs, and then WARNINGs. Once
their declarations are removed and they are hidden, packages using
them will fail to install.


Additionally, a blog post with "examples of moving non-API entry point

uses

into compliance" would be very helpful in this endeavor.


WRE now contains a section 'Moving into C API compliance'; that seems
a better option for the moment given that things are still very much
in flux. We will try to add to this section as needed. For the
specific entry points generating WARNINGs for your packages the advice
is simple: stop using them.

Best,

luke



   [[alternative HTML version deleted]]

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



--
Luke Tierney
Ralph E. Wareham Professor of Mathematical Sciences
University of Iowa  Phone: 319-335-3386
Department of Statistics andFax:   319-335-3017
 Actuarial Science
241 Schaeffer Hall  email:   luke-tier...@uiowa.edu
Iowa City, IA 52242 WWW:  http://www.stat.uiowa.edu/



[[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] [External] Non-API updates

[J C Nash]

It is probably important to note that the WRE with the new section on C API 
compliance is
in the R-Devel docs, not the current ones.

JN

On 2024-06-25 12:10, luke-tierney--- via R-devel wrote:

On Tue, 25 Jun 2024, Josiah Parry wrote:


Hey folks,

I'm sure many of you all woke to the same message I did: "Please correct
before 2024-07-09 to safely retain your package on CRAN" caused by Non-API
changes to CRAN.

This is quite unexpected as Luke Tierney's June 6th email writes (emphasis
mine):

"An *experimental* *effort* is underway to add annotations to the WRE..."

"*Once things have gelled a bit more *I hope to turn this into a blog
post that will include some examples of moving non-API entry point
uses into compliance."

Since then there has not been any indication of stabilization of the
Non-API changes nor has there been a blog post outlining how to migrate. As
things have been coming and going from the Non-API changes for quite some
time now, we (the royal we, here) have been waiting for an official
announcement from CRAN on the stabilizing changes.


I posted an update to this list a few days ago. If you missed it you
can find it in the archive.


*Can we extend this very short notice to handle the Non-API changes before
removal from CRAN? *


Timing decisions are up to CRAN.


In the case of the 3 packages I have to fix within 2 weeks, these are all
using Rust. These changes require upstream changes to the extendr library.
There are other packages that are also affected here. Making these changes
is a delicate act and requires care and focus. All of the extendr
developers work full time and cannot make addressing these changes their
only priority for the next 2 weeks.


Using non-API entry points is a choice that comes with risks. The ones
leading to WARNINGs for your packages (PRSEEN and SYMVALUE)have been
receiving NOTEs in check results for several weeks. Using
tools:::checkPkgAPI you can see that your packages are referencing a
lot of non-API entry points. Some of these may be added to the API,
but most will not. This may be a good time to look into that.

To minimize disruption we have been adding entry points to the API as
long as it is safe to do so, in some cases against our better
judgment. But ones that are unsafe to use will not be
added. Eventually their declarations will be removed from public
header files and they will be hidden when that is possible. Packages
that have chosen to use these non-API entry points will have to adapt
if they want to pass R CMD check. For now, we will try to first have
use of these entry points result in NOTEs, and then WARNINGs. Once
their declarations are removed and they are hidden, packages using
them will fail to install.


Additionally, a blog post with "examples of moving non-API entry point uses
into compliance" would be very helpful in this endeavor.


WRE now contains a section 'Moving into C API compliance'; that seems
a better option for the moment given that things are still very much
in flux. We will try to add to this section as needed. For the
specific entry points generating WARNINGs for your packages the advice
is simple: stop using them.

Best,

luke



[[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] request for discussion on lonely doc patch suggestion

[J C Nash]

For Linux users, meld is quite nice for side by side editing, though I've never 
tried using it for
display. Just checking now suggests it isn't obvious how to "print" side by 
side display.

I've made meld easier for my own use by creating an icon in Double Commander 
(DC allows
the user to create iconized links to scripts and programs). There are two panes 
in the DC
file manager. I highlight one file in each then click. This saves typing two 
full paths in
a command

   meld  path/to/file1 path/to/file2

I suspect the highlight and click makes my use of meld reasonably attractive. 
I'm not sure
I'd use it in the raw command line mode.

Like Duncan, I welcome suggestions for similar tools, especially if there's a 
display option.

John Nash

On 2025-03-24 15:21, Duncan Murdoch wrote:

I sent some comments directly to Ben.  I just want to reply publicly to this 
part:

On 2025-03-24 1:18 p.m., Ben Bolker wrote:


    The patch file is attached (also available at bugzilla, if it doesn't
get through to the list). I find the patch format a little hard to read,
so I'm reproducing just the *new* text below.


I agree absolutely about the lack of readability of patch files.  A side by side display is much nicer.  If anyone out 
there isn't using one, you should.


I really like the one I use ("Beyond Compare"), but it's not open source.  I've been using it for a very long time (20 
years or more, I think), and I suspect there are very good open source competitors out there now (and may have been for 
all the time I've been using BC). Suggestions?


Duncan Murdoch

__
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] Model forecasts with new factor levels - predict.warn

[Prof J C Nash]

predict.warn() -- a function to display factor levels in new data
for linear model prediction that do not exist in the
estimating data.

Date: 2005-8-24
From: John C. Nash (with thanks to Uwe Ligges for suggestions)
[EMAIL PROTECTED]

Motivation: In computing predictions from a linear model using factors,
it is possible to introduce new factor levels. This was encountered on
a practical level in forecasting retail sales where certain events
had not occurred during the estimation period but did occur during
the validation (i.e., trial forecast) period. predict() then gives
an error message that it cannot continue because the new data has
new levels of a factor. Indeed, it will give the new level of the
first such factor found. predict.warn is designed to give more
detailed information. for example, if there are other factor levels
that are also "new".

The underlying interest is to find a suitable method for introducing
estimates of the effects of such "new" factor levels into the model,
and I welcome exchanges on how that may be done. Indeed, one of my
main interests in submitting this tool is to learn how that may be
done appropriately and with suitable information to the user to avoid
egregious errors. For example, a new sales promotion may have many
similarities with existing promotions, so a rough estimate for the
effect may be possible, allowing for rudimentary validation testing
and for "working" values of forecasts.

Unfortunately, the lack of some facility to provide "guesstimates"
in this way leads users to ignore R or other tools and to work with
Excel. I have experience of this situation arising in practise, and
we may anticipate that these users will be difficult to convert to
using better tools. I believe we must find ways to allow users to
carry out such operations in as safe and informed manner possible.
Ideally, the follow-on to predict.warn should be a function called
something like predict.expand or predict.force that asks the user
to supply values (or possibly sets or distributions) of the "new"
factor level estimates.



Sample function:

predict.warn<-function (object, newdata,
  ...)
# this is function predict.warn.R
# Tries to warn if we do not have appropriate factor levels in 
estimation data
# that appear in the newdata.
# It does NOT yet test for non-factor variables, but that would also be 
nice.
# However, we usually assume that numerical variables can be 
interpolated and
# extrapolated.
{   message("predict.warn -- check for factor levels not present")
 message("in linear models with factors")
 message(" ")
 if (missing(newdata) || is.null(newdata)) {
# Without newdata, issue warning and return
message("You do not have any new data.")
 }
 else {
 xlev <- object$xlevels
 nn<-names(xlev)
dnn<-length(nn)
message("there are ",dnn," factor names ")
dostop<-FALSE # Ensure the return value is defined
for (i in seq(along=nn)) {
message("Factor ",nn[i]," in estimation data")
ofac <- object[["model"]][[nn[i]]]
print(ofac)  # diagnostic
oldlev<-levels(factor(ofac))
nfac <- newdata[[nn[i]]]  # ?? note need for double brackets. 
Why?
message("Factor ",nn[i], " in new data")
print(nfac) # diagnostic
message("about to try levels on newfac")
newlev<-levels(factor(nfac))
for (j in seq(along=newlev)) {
if (!(newlev[j] %in% oldlev)) {
message("New level ", newlev[j]," not found in 
estimation levels of 
factor ", nn[i])
dostop = TRUE
}
}
    } # end loop over names
 } # end part where we DO have new data
 return(dostop) # TRUE if we should stop
}

# end of predict.warn code

Test example:

# A simple test of predict.warn.R for 2 factors
# J C Nash 20050822
rm(list=ls())
source("predict.warn.R")
x<-c(1, 2, 3, 4, 5, 6, 7, 8)
fac1<-c("A", "A", "", "", "", "B", "")
fac1<-c(fac1, "A")
fac1
y<-c( 2, 4, 5, 8, 9.5, 12, 13, 17)
fac2<-c(1,1,1,2,2,1,2,3)
fac2<-as.factor(fac2) # Note: or else defaults to num
joe<-data.frame(x, y, fac1, fac2)
bill<-subset(joe, x<=5)
mary<-subset(joe, x>5)
message("bill")
str(bill)
print(bill)
message("mary")
str(mary)
print(mary)
est<-lm(y ~ fac1 + fac2 + x, bill)
print(est)
test<-predict.warn(est, newdata = mary)
message('predict.warn returns ',test)
message("\n")
message(&qu

Re: [Rd] if(!CRAN()){...} / Repository policies

[U30A J C Nash]

If CRAN were a passive repository, the discussion about its policies 
would not be relevant to this list e.g., SourceForge. However, the 
development of R and its packages are very intimately connected to the 
CRAN repository policy.


I doubt any of the players in building our current R ecosystem ever 
imagined it would become this big. The serious implications of 
repository and related policies for future development and health of 
this ecosystem and its community members deserve transparent debate. As 
with any movement based on volunteer contribution, R will adapt or be 
replaced.


As far as I'm aware -- and welcome correction -- we don't have a list to 
discuss policies by CRAN or other actors. This list seems appropriate 
because it generally involves those doing development for R. Whether 
R-core should or should not listen to such deliberations is also open to 
debate. My own view -- which could change depending on how things evolve 
-- is that a small central committee is fine as long as they are able to 
both support and be supported by the wider development community that 
provides the packages giving R its strength. Unfortunately, I sense that 
the time demands of what is now a major software system are such that 
communications (and tempers) are getting frayed. I urge attention to 
this, and also to ideas of organizational behaviour and evolution, so we 
keep R healthy. Those who know me will realize I'm willing to put in 
time to act as the glue in building communications.


I'm also a bit surprised nobody has mentioned the mirrors -- if the 
mirror repositories choose, say, the Alternative R Archive Network, 
which might also include sweater patterns, the effectiveness of the CRAN 
system would be crippled.


John Nash

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

Re: [Rd] What algorithm is R using to calculate mean?

[Prof J C Nash (U30A)]

There is quite a literature on related methods for variance. If anyone 
is interested, I did some work (and even published the code in the 
magazine Interface Age in 1981) on some of these. I could probably put 
together scans of relevant materials, some of which are not easily 
available. It would make a good Master's or senior undergrad project, 
and could be (possibly has been) extended to compute covariance matrices 
properly for very large data sets and streaming data.


JN

On 13-07-27 06:00 AM, r-devel-requ...@r-project.org wrote:

Message: 6
Date: Fri, 26 Jul 2013 12:58:28 +
From: Ravi Varadhan
To: Joshua Ulrich, Zach Harrington

Cc:"r-devel@r-project.org List"  
Subject: Re: [Rd] What algorithm is R using to calculate mean?
Message-ID:
<2f9ea67ef9ae1c48a147cb41be2e15c3491...@dom-eb-mail1.win.ad.jhu.edu>
Content-Type: text/plain; charset="us-ascii"

This uses the idea of Kahan's summation, if I am not mistaken.

http://en.wikipedia.org/wiki/Kahan_summation_algorithm

Ravi
__


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

Re: [Rd] Version of L-BFGS-B used in optim etc

[Prof J C Nash (U30A)]

This issue has been known for some time and I've had "why don't you fix 
this?" queries. However, I'm not one of the R-core folk who could do so, 
and don't code in C.  Moreover, as far as I can tell, the version of 
L-BFGS-B in R is not one of the standard releases from Morales and Nocedal.


As maintainer of optimx, I can state that I won't be including L-BFGS-B 
except through optim() unless someone implements the new code as a 
package, which I would welcome. I've decided to focus only on codes 
written all in R. For bounds constrained optimization, I've put together 
Rvmmin which is the "BFGS" method of optim() (a misnomer, but I'd better 
not get started...) and Rcgmin. I've got my brother's truncated Newton 
code running unconstrained, but still trying to find time to get the 
bounds debugged. I find codes in R more flexible in allowing 
improvements, and easier to avoid those nasty mismatches in coding style 
that happen across languages.


I'd be glad to collaborate with anyone on getting latest L-BFGS-B 
packaged, or on improvements to the TNmin etc. And, of course, I welcome 
examples where things don't work smoothly in my codes.


There are upgrades to optimx on R-forge that will not likely move to 
CRAN for a few months until TN is completed and added etc.


Best, JN


On 13-10-10 06:00 AM, r-devel-requ...@r-project.org wrote:

Message: 1
Date: Wed, 9 Oct 2013 15:46:40 +1000
From: David Duffy
To:
Subject: [Rd] Version of L-BFGS-B used in optim etc
Message-ID:
Content-Type: text/plain; format=flowed; charset="US-ASCII"

Hi.

I just noticed the paper by Morales and Nocedal

Remark on "Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale
Bound Constrained Optimization". TOMS 2011; 38(1): 7

http://www.ece.northwestern.edu/~morales/PSfiles/acm-remark.pdf

which describes a couple of improvements (speed and accuracy) to the
original Netlib code which AFAICT is that still used by optim()
via f2c.  Updated code is under

http://www.ece.northwestern.edu/~nocedal/lbfgsb.html

released under the New BSD License.  Has this already been made available
in R, perhaps in other packages such as optimx?

Cheers, David Duffy.


| David Duffy (MBBS PhD) ,-_|\
| email:dav...@qimr.edu.au   ph: INT+61+7+3362-0217 fax: -0101  / *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v


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

[Rd] Byte code compile not helpful in R3.0.2

[Prof J C Nash (U30A)]

I had a bunch of examples of byte code compiles in something I was
writing. Changed to 3.0.2 and the advantage of compiler disappears. I've
looked in the NEWS file but do not see anything that suggests that the
compile is now built-in. Possibly I've just happened on a bunch of
examples where it does not help, but experiences of a year ago do not
seem to remain valid now. Just wondering if my experience was consistent
with what is expected now in 3.0.2.

John Nash

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

Re: [Rd] Byte code compile not helpful in R3.0.2

[Prof J C Nash (U30A)]

My bad to not give details. I'm comparing (though not quite directly) to
results in the posting
http://rwiki.sciviews.org/doku.php?id=tips:rqcasestudy.

What prompted the query was a write up of "for" versus "while" loops,
where there was a speedup using compiler for one of these. I had the
example in a knitr file, and when I was reviewing the text before
sending it to an editor, I realized the timings no longer supported the
text. They were, as I recall, developed in R 2.15.2, and I just looked
through my VMs with different OS's to see if there is one with that
still extant, but except for a real Win7 case I have been too "good" and
updated to at least 3.0.1, where I'm getting no advantage. The Win7 case
is R 2.15.1, and there the compiler actually went slower on one run of
the code below. That may be due to antivirus running -- had not booted
that partition for quite a while.

Here is the for-while test code:

#  forwhiletime.R
library(microbenchmark)
require(compiler)

tfor <- function(n){
for (i in 1:n) {
   xx<-exp(sin(cos(as.double(i
}
xx
}

twhile <- function(n){
    i<-0
while (i On 13-11-03 2:15 PM, Prof J C Nash (U30A) wrote:
>> I had a bunch of examples of byte code compiles in something I was
>> writing. Changed to 3.0.2 and the advantage of compiler disappears. I've
>> looked in the NEWS file but do not see anything that suggests that the
>> compile is now built-in. Possibly I've just happened on a bunch of
>> examples where it does not help, but experiences of a year ago do not
>> seem to remain valid now. Just wondering if my experience was consistent
>> with what is expected now in 3.0.2.
> 
> Post some details, please.  Are the times in 3.0.2 like the times in
> 3.0.1 with or without compiling?  Or were you comparing to some other
> version?
> 
> Duncan Murdoch
> 
> 
>

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

Re: [Rd] Byte code compile (not) helpful in R3.0.2 -- Fixed.

[Prof J C Nash (U30A)]

Thanks. I should not try adjusting code after some hours of proofreading.

Making that change gave a suitable time difference.

Best, JN



On 13-11-03 03:46 PM, Henrik Bengtsson wrote:
> tfor <- cmpfun(tfor)
> twhile <- cmpfun(twhile)
> 
> /Henrik
> 
> 
> On Sun, Nov 3, 2013 at 11:55 AM, Prof J C Nash (U30A)  
> wrote:
>> My bad to not give details. I'm comparing (though not quite directly) to
>> results in the posting
>> http://rwiki.sciviews.org/doku.php?id=tips:rqcasestudy.
>>
>> What prompted the query was a write up of "for" versus "while" loops,
>> where there was a speedup using compiler for one of these. I had the
>> example in a knitr file, and when I was reviewing the text before
>> sending it to an editor, I realized the timings no longer supported the
>> text. They were, as I recall, developed in R 2.15.2, and I just looked
>> through my VMs with different OS's to see if there is one with that
>> still extant, but except for a real Win7 case I have been too "good" and
>> updated to at least 3.0.1, where I'm getting no advantage. The Win7 case
>> is R 2.15.1, and there the compiler actually went slower on one run of
>> the code below. That may be due to antivirus running -- had not booted
>> that partition for quite a while.
>>
>> Here is the for-while test code:
>>
>> #  forwhiletime.R
>> library(microbenchmark)
>> require(compiler)
>>
>> tfor <- function(n){
>> for (i in 1:n) {
>>xx<-exp(sin(cos(as.double(i
>> }
>> xx
>> }
>>
>> twhile <- function(n){
>> i<-0
>> while (i>i<-i+1
>>xx<-exp(sin(cos(as.double(i
>> }
>> xx
>> }
>> n<-1
>>
>> timfor<-microbenchmark(tfor(n))
>> timwhile<-microbenchmark(twhile(n))
>> timfor
>> timwhile
>> cmpfun(tfor)
>> cmpfun(twhile)
>> timforc<-microbenchmark(tfor(n))
>> timwhilec<-microbenchmark(twhile(n))
>> timforc
>> timwhilec
>> looptimes<-data.frame(timfor$time, timforc$time, timwhile$time,
>> timwhilec$time)
>> colMeans(looptimes)
>>
>>
>> Actually, I'm not greatly axious about all this. Mainly I want to make
>> sure that I get whatever advice is to be rendered so it is correct.
>>
>> Best,
>>
>> JN
>>
>>
>> On 13-11-03 02:22 PM, Duncan Murdoch wrote:
>>> On 13-11-03 2:15 PM, Prof J C Nash (U30A) wrote:
>>>> I had a bunch of examples of byte code compiles in something I was
>>>> writing. Changed to 3.0.2 and the advantage of compiler disappears. I've
>>>> looked in the NEWS file but do not see anything that suggests that the
>>>> compile is now built-in. Possibly I've just happened on a bunch of
>>>> examples where it does not help, but experiences of a year ago do not
>>>> seem to remain valid now. Just wondering if my experience was consistent
>>>> with what is expected now in 3.0.2.
>>>
>>> Post some details, please.  Are the times in 3.0.2 like the times in
>>> 3.0.1 with or without compiling?  Or were you comparing to some other
>>> version?
>>>
>>> Duncan Murdoch
>>>
>>>
>>>
>>
>> __
>> 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] cat with backspace and newline characters

[Prof J C Nash (U30A)]

Over the years, this has been useful to me (not just in R) for many
nonlinear optimization tasks. The alternatives often clutter the screen.

> On 13-11-06 06:00 AM, r-devel-requ...@r-project.org wrote:

> People do sometimes use this pattern for displaying progress (e.g. iteration 
> counts). 
>> 
>> 


As this looks like a "terminal" matter, the first level solution may be
to take a couple of examples and prepare a table of what happens on the
most used platforms, including in things like RStudio and ESS. Might be
a good issue for Rwiki, as I'm sure few of us have all the choices.

I can generally live with warning users when the iteration display may
display oddly. Perhaps others have more stringent requirements, but if
documentation is enough, we can avoid fixing something that is more or
less outside R and may be unnecessary work. And documented examples can
also show us if things are changing.

JN

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

Re: [Rd] R CMD check for the R code from vignettes -- thread fraying?

[Prof J C Nash (U30A)]

I noted Duncan's comment that an answer had been provided, and went to
the archives to find his earlier comment, which I am fairly sure I saw a
day or two ago. However, neither May nor June archives show Duncan in
the thread except for the msg below (edited for space). Possibly tech
failures are causing misunderstandings.

JN

On 14-06-02 06:00 AM, r-devel-requ...@r-project.org wrote:
> Message: 4 Date: Mon, 02 Jun 2014 14:06:28 +0900 From: Duncan Murdoch
>  To: Carl Boettiger ,
> Gabriel Becker  Cc: Henrik Bengtsson
> , R-devel  Subject: Re: [Rd]
> R CMD check for the R code from vignettes Message-ID:
> <538c0654.7050...@gmail.com> Content-Type: text/plain;
> charset=ISO-8859-1; format=flowed On 02/06/2014, 1:41 PM, Carl Boettiger
> wrote:
>> > Thanks both for the replies.
>> >


> You haven't been reading very carefully.  I saw several:  mine, Martin 
> Morgan's, Kasper Daniel Hansen's.
> 
> Duncan Murdoch



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

[Rd] a question about optim.R and optim.c in R

[Prof J C Nash (U30A)]

As you dig deeper you will find vmmin.c, cgmin.c and (I think) nmmin.c
etc. Those were, as I understand, converted by p2c from my Pascal codes
that you can find in the pascal library on netlib.org. These can be run
with the Free Pascal compiler.

Given how long ago these were developed (>30 years in all cases), they
are due for review. The packages Rvmmin and Rcgmin are all-R
replacements for the 1st two, and nmkb from dfoptim generally offers a
better version of the Nelder-Mead approach. All have bounds-constrained
variants, but for efficiency, there are direct calls of the
unconstrained methods, though I need to provide some nice examples of
when and how to call each. There is likely a place for some compiling of
sections to speed things up, but the R codes are not particularly sluggish.

Side comment: At UseR last week, Yihui Xie sat with me and we
implemented a Fortran language engine for knitr. It looks like a Pascal
one may also be possible, and maybe even a BASIC (though that may need
variants for different platforms). This would allow vignettes to
document some of the legacy code to be written, and that may be an
important matter as the expertise for such older tools moves into
retirement. Off-list communication about such ideas welcome.

John Nash


On 14-07-08 06:00 AM, r-devel-requ...@r-project.org wrote:
> Message: 2
> Date: Mon, 7 Jul 2014 16:34:59 -0400
> From: Zhiyuan Dong 
> To: r-devel@r-project.org
> Subject: [Rd] a question about optim.R and optim.c in R
> Message-ID:
>   
> Content-Type: text/plain
> 
> Hi, I am learning R by reading R source code. Here is one question I have
> about the optim function in R.
> 
> The context : In the optim.R, after all the prep steps, the main function
> call call is made via :
> 
> .External2(C_optim, par, fn1, gr1, method, con, lower, upper).
> 
> So, it seems to me, to follow what is going on from here, that I should
> read the optim function in \src\library\stats\src\optim.c
> 
> where it has this signature :
> 
> SEXP optim(SEXP call, SEXP op, SEXP args, SEXP rho)
> 
> I am not sure I follow here : In the .External2 call, we have 7 parameters
> :  par, fn1, gr1, method, con, lower, upper; This does not seem to match
> the signature of
> 
> SEXP optim(SEXP call, SEXP op, SEXP args, SEXP rho)
> 
> However, it seems (from the source code) that the 7 parameters are somehow
> embedded in the 'args' parameter. I am not sure what is going on...Am I
> missing something?
> 
> Thanks much!!!
> 
> Best,
> 
> Zhiyuan
> 
>   [[alternative HTML version deleted]]

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

Re: [Rd] lbfgsb from C/C++

[Prof J C Nash (U30A)]

I won't comment on the C/C++ option, as I'm not expert in that. However,
R users and developers should know that Nocedal et al. who developed
L-BFGS-B released an update to correct a fault in 2011. It was important
enough that an ACM TOMS article was used for the announcement.

I recently implemented a wrapper to the new Fortran code. As it is a
reverse communication code, the main task was getting rid of all the
Fortran WRITE() statements (not a happy experience!).

Maybe Dirk has ideas on how to make the wrapper more efficient and the
call from C or C++ direct for those who need that.

The code will be incorporated eventually in the optimx package, but for
now an experimental (NOTE - EXPERIMENTAL) version lbfgsb3 is up on
r-forge under the optimizer project. I'd be happy to exchange ideas
off-list on this.

Here's the reference to the TOMS announcement:

@Article{Morales:2011:FSL,
  author =   "Jos\'e Luis Morales and Jorge Nocedal",
  title ="Remark on ``{Algorithm} 778: {L-BFGS-B}: {Fortran}
Subroutines
  for Large-Scale Bound Constrained Optimization''",
  journal =  "{ACM} Transactions on Mathematical Software",
  accepted = "15 April 2011",
  volume =   "38",
  number =   "1",
  month =nov,
  URL =  "http://doi.acm.org/10.1145/2049662.2049669";,
  pages ="7:1--7:4",
  year = "2011",
  abstract = "
 This remark describes an improvement and a correction
 to Algorithm 778. It is shown that the performance of
 the algorithm can be improved significantly by making
 a relatively simple modification to the subspace
 minimization phase. The correction concerns an error
 caused by the use of routine dpmeps to estimate machine
 precision.",
}

John Nash

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

[Rd] How can I use R-function in My C++ project ? (using optim)

[Prof J C Nash (U30A)]

As the author of 3 of the 5 methods in optim, I think you may be wasting
your time if this is for performance. My reasons are given in

http://www.jstatsoft.org/v60/i02

Note that most of the speed benefits of compilation are found in the
objective and gradient function, with generally more minor improvements
in the method.

JN


On 14-10-22 06:00 AM, r-devel-requ...@r-project.org wrote:
> Date: Tue, 21 Oct 2014 20:14:03 +0800 (CST)
> From: ?? <2012111...@njau.edu.cn>
> To: R-devel@r-project.org
> Subject: [Rd] How can I use R-function in My C++ project ?
> Message-ID: <1f74c14.d1e7.14932a0ddbf.coremail.2012111...@njau.edu.cn>
> Content-Type: text/plain; charset="UTF-8"
> 
> Dear seniors:
>I am a student in Nanjing Agricultural University of China.
> 
> 
>I want to use the  function "optim"  of package stats in my C++ project. I 
> have got the R.dll , R.def and R.lib,
> 
> 
> but I can't find the function prototypes of "optim" in R.def. 
>   
>How  can  I  do ?  Is the Method I call R function in C++ with R.dll 
> feasible ?  I hope to get your help ! Thanks 
> 
> 
> for reading my question.
> 
> 
> Yours sincerely,
> 
> 
> Bo Huang
>   [[alternative HTML version deleted]]

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

[Rd] Help finding source of warnings

[Prof J C Nash (U30A)]

I've been implementing a wrapper to the 2011 Fortran version of 
L-BFGS-B. In optim(), R uses a C translation of a Fortran version (the 
version number does not appear to be documented by the original 
authors). The authors of the original Fortran code have updated it and 
published the reasons in ACM TOMS due to inefficiencies and a bug.


In running the checks on the resulting package (which is on R-forge 
under the optimizer project), I'm getting a number of warning messages 
of the type


Warning in file.copy(file.path(.Library, pkg, "DESCRIPTION"), pd) :
  problem copying /usr/lib/R/library/mgcv/DESCRIPTION to 
/tmp/Rtmp0kkeHo/RLIBS_1214765d1c5f/mgcv/DESCRIPTION: No such file or 
directory


which reference DESCRIPTIONs for a number of packages other than the one 
being checked -- here mgcv -- and which are not referenced in my package 
as far as I can determine.


Possibly unrelated, when I run the code on a problem, it works for one 
run, then gives a NAMESPACE error and failure on the second try. Apart 
from this, checks and unit tests appear to work correctly.


Does anyone have pointers where I might find some ideas on the origin of 
the issue(s)? I suspect the warning messages are not particularly 
indicative of the source of the warnings, but that I have some subtle 
glitch in the setup and call to the Fortran.


I suspect this is not platform dependent, but I'm running Linux Mint 
17.1 (ubuntu derivative), and  R 3.1.2.


Cheers, JN

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

Re: [Rd] Help finding source of warnings

[Prof J C Nash (U30A)]

Kurt pointed to the issue. Thanks. I did install r-recommended, but it 
seems something went wrong at some point. A reinstall got rid of the 
warnings.


Thanks to Dirk for his offer of help.

Now I'm still getting a namespace issue on the second run of an 
optimization problem. However, I think I need to do some more digging to 
narrow down where this issue is lurking. It may be some local matter, as 
with the r-recommended links failing.


Best, JN

On 15-01-18 09:27 AM, Kurt Hornik wrote:

Prof J C Nash (U30A) writes:



I've been implementing a wrapper to the 2011 Fortran version of
L-BFGS-B. In optim(), R uses a C translation of a Fortran version (the
version number does not appear to be documented by the original
authors). The authors of the original Fortran code have updated it and
published the reasons in ACM TOMS due to inefficiencies and a bug.



In running the checks on the resulting package (which is on R-forge
under the optimizer project), I'm getting a number of warning messages
of the type



Warning in file.copy(file.path(.Library, pkg, "DESCRIPTION"), pd) :
problem copying /usr/lib/R/library/mgcv/DESCRIPTION to
/tmp/Rtmp0kkeHo/RLIBS_1214765d1c5f/mgcv/DESCRIPTION: No such file or
directory



which reference DESCRIPTIONs for a number of packages other than the one
being checked -- here mgcv -- and which are not referenced in my package
as far as I can determine.



Possibly unrelated, when I run the code on a problem, it works for one
run, then gives a NAMESPACE error and failure on the second try. Apart
from this, checks and unit tests appear to work correctly.



Does anyone have pointers where I might find some ideas on the origin of
the issue(s)? I suspect the warning messages are not particularly
indicative of the source of the warnings, but that I have some subtle
glitch in the setup and call to the Fortran.



I suspect this is not platform dependent, but I'm running Linux Mint
17.1 (ubuntu derivative), and  R 3.1.2.


John: maybe you did not install the recommended packages?
(On Debian, the corresponding package would be r-recommended.)

Best
-k


Cheers, 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] Requirement for pandoc 1.12.3 in R 3.1.3

[Prof J C Nash (U30A)]

Are other developers finding R 3.1.3 problematic because vignette
building requires pandoc 1.12.3, while Linux Mint 17 / Ubuntu 14.04 have
1.12.2.1? R 3.1.2 seems to work fine.

I'd very much like to avoid having to build as large a Linux package as
pandoc, which has given me issues outside of R (it leaves out words,
sentences or paragraphs when converting Latex to epub in a novel I'm
working on, and does so without warning). Possibly concerns like this
are why R has moved to a later pandoc.

Suggestions welcome.

John Nash

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

Re: [Rd] Requirement for pandoc 1.12.3 in R 3.1.3

[Prof J C Nash (U30A)]

the derivatives table
>> Error in deriv.formula(~abs(x), "x") :
>>   Function 'abs' is not in the derivatives table
>> Error in D(expression(sign(x)), "x") :
>>   Function 'sign' is not in the derivatives table
>> Error in deriv.formula(~sign(x), "x") :
>>   Function 'sign' is not in the derivatives table
>> sh: 1: yacas: not found
>> Quitting from lines 629-641 (nls14tutorial.Rmd)
>> Error: processing vignette 'nls14tutorial.Rmd' failed with
>> diagnostics:
>> cannot open the connection
>> Execution halted
>>
>> For reference:
>>
>> john@john-J6-2015 ~/current/nls14work $ cat /etc/os-release
>> NAME="Ubuntu"
>> VERSION="14.04.2 LTS, Trusty Tahr"
>> ID=ubuntu
>> ID_LIKE=debian
>> PRETTY_NAME="Ubuntu 14.04.2 LTS"
>> VERSION_ID="14.04"
>> HOME_URL="http://www.ubuntu.com/";
>> SUPPORT_URL="http://help.ubuntu.com/";
>> BUG_REPORT_URL="http://bugs.launchpad.net/ubuntu/";
>>
>> john@john-J6-2015 ~/current/nls14work $ cat /etc/linuxmint/info
>> RELEASE=17.1
>> CODENAME=rebecca
>> EDITION="MATE 64-bit"
>> DESCRIPTION="Linux Mint 17.1 Rebecca"
>> DESKTOP=MATE
>> TOOLKIT=GTK
>> NEW_FEATURES_URL=http://www.linuxmint.com/rel_rebecca_mate_whatsnew.php
>> RELEASE_NOTES_URL=http://www.linuxmint.com/rel_rebecca_mate.php
>> USER_GUIDE_URL=help:linuxmint
>> GRUB_TITLE=Linux Mint 17.1 MATE 64-bit
>>
>>
>>
>> john@john-J6-2015 ~/current/nls14work $ dpkg -l | grep -i pandoc
>> ii  libghc-pandoc-citeproc-data 0.2-3build1
>> all  Pandoc support for Citation
>> Style Language - data files
>>
>> ii  pandoc  1.12.2.1-1build2
>> amd64general markup converter
>>
>> ii  pandoc-citeproc 0.2-3build1
>> amd64Pandoc support for Citation
>> Style Language - tools
>>
>> ii  pandoc-data 1.12.2.1-1build2
>> all  general markup converter -
>> data
>> files
>>
>> JN
>>
>>
>> On 15-03-12 10:21 AM, Prof Brian Ripley wrote:
>>> On 12/03/2015 13:51, Prof J C Nash (U30A) wrote:
>>>> Are other developers finding R 3.1.3 problematic because vignette
>>>> building requires pandoc 1.12.3, while Linux Mint 17 / Ubuntu
>>>> 14.04 have
>>>> 1.12.2.1? R 3.1.2 seems to work fine.
>>>
>>> R has no built-in support for non-Sweave vignettes, and there is no
>>> mention of pandoc in the R 3.1.3 sources except for the manual:
>>>
>>> 'Complete checking of a package which contains a file
>>> @file{README.md}
>>> needs @command{pandoc} installed: see
>>> @uref{http://johnmacfarlane.net/@/pandoc/@/installing.html}.'
>>>
>>> which is true (but is not done with R 3.1.3).
>>>
>>> I suspect you are confusing an R update with an update of whatever
>>> packages you use to process your vignettes: package rmarkdown has a
>>> pandoc version requirement of
>>>
>>> SystemRequirements: pandoc (>= 1.12.3) -
>>>
>>
>> __
>> 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] nls

[Prof J C Nash (U30A)]

nls() is using
1) only a Gauss-Newton code which is prone to some glitches
2) approximate derivatives

Package nlmrt uses symbolic derivatives for expressions (you have to
provide Jacobian code for R functions) and an aggressive Marquardt
method to try to reduce the sum of squares. It does return more
information about the problem (singular values of the final Jacobian
and gradient at the proposed solution) but does NOT return the nls
structured object. And it will usually take more time and computing
effort because it tries hard to reduce the SS.

A reproducible example would get you a more informed response.

John Nash


On 15-03-19 07:00 AM, r-devel-requ...@r-project.org wrote:
> Date: Wed, 18 Mar 2015 14:14:12 +0200
> From: Evans Otieno Ochiaga 
> To: r-devel@r-project.org
> Subject: [Rd] Help
> Message-ID:
>   
> Content-Type: text/plain; charset="UTF-8"
> 
> Hi to All,
> 
> I am fitting some models to a data using non linear least square, and
> whenever i run the command, parameters value have good convergence but I
> get the  error in red as shown below. Kindly how can I fix this problem.
> 
> 
> Convergence of parameter values
> 
> 0.2390121 :  0.1952981 0.975 1.000
> 0.03716107 :  0.1553976 0.910 1.000
> 0.009478433 :  0.2011017 0.798 1.000
> 0.004108196 :  0.2640111 0.693 1.000
> 0.003705189 :  0.2938360 0.652 1.000
> 0.003702546 :  0.2965745 0.650 1.000
> 0.003702546 :  0.2965898 0.650 1.000
> 0.003702546 :  0.2965898 0.650 1.000
> 0.003702546 :  0.2965898 0.650 1.000
> 
> Error in nls(Occupancy ~ 1 - (theta * beta^(2 * Resolution^(1/2)) *
> delta^Resolution),  :
>   step factor 0.000488281 reduced below 'minFactor' of 0.000976562
> 
> Regards,
> 
> 
> 
> 
> *Evans Ochiaga*

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

[Rd] Title case in DESCRIPTION for package where a word is a function name

[Prof J C Nash (U30A)]

I was preparing a fix for a minor glitch in my optimx package and R CMD
check gave an error that the title was not in title case. It is

A Replacement and Extension of the optim() Function

R CMD check suggests the incorrect form

A Replacement and Extension of the Optim() Function

'Writing R Extensions' suggests single quotes, i.e.,

A Replacement and Extension of the 'optim()' Function

which R CMD check still complains about.

I have found

A Replacement and Extension of the _optim()_ Function

does not get the complaint, but I'm not sure the underscore is allowed.

Given that I've obeyed the RTFM rule, I'm wondering what to do now.

On a related matter, I'm finding the reverse dependency check for optimx
takes a very long time and sometimes stalls for reasons I have not yet
sorted out. I run it in virtual machines for R3.2 and R-devel. Possibly
optimx needs so many packages I'm hitting memory or disk limits. Perhaps
off-list discussion could suggest a way to set up a reverse check
server. I'd be willing to help on such a project, which might be helpful
for many developers.

JN

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

Re: [Rd] Title case in DESCRIPTION for package where a word is a function name

[Prof J C Nash (U30A)]

Hendrik pointed out it was the parentheses that gave the complaint.
Single quotes and no parentheses seem to satisfy R CMD check. Perhaps
that needs to be in the WRE.

However, I have for some time used the parentheses to distinguish
functions from packages. optim() is a function, optimx a package.
Is this something CRAN should be thinking about? I would argue greater
benefit to users than title case.

JN


On 15-04-24 06:17 PM, Uwe Ligges wrote:
> 
> 
> On 24.04.2015 22:44, Ben Bolker wrote:
>> Prof J C Nash (U30A  uottawa.ca> writes:
>>
>>>
>>> I was preparing a fix for a minor glitch in my optimx package and R CMD
>>> check gave an error that the title was not in title case.
>>
>>[snip] to make Gmane happy ...
>>
>>> I have found
>>>
>>> A Replacement and Extension of the _optim()_ Function
>>>
>>> does not get the complaint, but I'm not sure the underscore is allowed.
>>>
>>> Given that I've obeyed the RTFM rule, I'm wondering what to do now.
>>
>>Presumably you should ask the CRAN maintainers?  That seems to
>> be the only possible answer -- I don't think anyone else can guess
>> very accurately ...
> 
> From WRE:
> 
> "Refer to other packages and external software in single quotes, and to
> book titles (and similar) in double quotes."
> 
> Other non-English usage (as documented for the Description field; this
> inlcudes function names) can also be used in single quotes.
> 
> Best,
> Uwe Ligges
> 
> 
>>
>>Ben Bolker
>>
>> __
>> 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] Title case in DESCRIPTION for package where a word is a function namei

[Prof J C Nash (U30A)]

How about allowing underscore? (I believe WRE is silent on this, and I
have not tried submitting a package with underscore in the title.) As I
pointed out in my OP, _optim()_ works. And we have the advantage that we
can distinguish package from function.

The purpose of consistent editing is surely to provide the affordances
that save us from needing extra documentation, as per Donald Norman's
excellent discussions on Design of Everyday Things, or Turn Signals are
the Facial Expressions of Automobiles. Changing the name of a function
in a case-sensitive computing language may not be a bug, but it is
asking for trouble.

JN

On 15-04-25 07:57 AM, peter dalgaard wrote:
> 
>> On 25 Apr 2015, at 13:11 , Prof J C Nash (U30A)  wrote:
>>
>> Hendrik pointed out it was the parentheses that gave the complaint.
>> Single quotes and no parentheses seem to satisfy R CMD check. Perhaps
>> that needs to be in the WRE.
> 
> Well, it is in ?toTitleCase:
> 
>  ...However, unknown
>  technical terms will be capitalized unless they are single words
>  enclosed in single quotes: names of packages and libraries should
>  be quoted in titles.
> 
> ..and it is the "single word" bit that gets you. AFAICT, the issue is that it 
> splits the text into words and then looks for words that begin and end with a 
> single quote. And parentheses count as word separators, so the quotes of 
> 'optim()' end up in two different words. 
> 
> It's one of those things that aren't easy to fix: Presumably you do want 
> capitalization within parentheses so we can't just not let them be 
> separators, and we can't just look for sets of single quotes with arbitrary 
> content because they get used inside ordinary text (e.g. the beginning of 
> this paragraph contains 's one of those things that aren'). So either we need 
> more heuristics, like only counting () as separators when preceded by or 
> preceding a space, or some sort of explicit escape mechanism, like BibTeX's 
> {foo}.
> 
>>
>> However, I have for some time used the parentheses to distinguish
>> functions from packages. optim() is a function, optimx a package.
>> Is this something CRAN should be thinking about? I would argue greater
>> benefit to users than title case.
>>
>> JN
>>
>>
>> On 15-04-24 06:17 PM, Uwe Ligges wrote:
>>>
>>>
>>> On 24.04.2015 22:44, Ben Bolker wrote:
>>>> Prof J C Nash (U30A  uottawa.ca> writes:
>>>>
>>>>>
>>>>> I was preparing a fix for a minor glitch in my optimx package and R CMD
>>>>> check gave an error that the title was not in title case.
>>>>
>>>>   [snip] to make Gmane happy ...
>>>>
>>>>> I have found
>>>>>
>>>>> A Replacement and Extension of the _optim()_ Function
>>>>>
>>>>> does not get the complaint, but I'm not sure the underscore is allowed.
>>>>>
>>>>> Given that I've obeyed the RTFM rule, I'm wondering what to do now.
>>>>
>>>>   Presumably you should ask the CRAN maintainers?  That seems to
>>>> be the only possible answer -- I don't think anyone else can guess
>>>> very accurately ...
>>>
>>> From WRE:
>>>
>>> "Refer to other packages and external software in single quotes, and to
>>> book titles (and similar) in double quotes."
>>>
>>> Other non-English usage (as documented for the Description field; this
>>> inlcudes function names) can also be used in single quotes.
>>>
>>> Best,
>>> Uwe Ligges
>>>
>>>
>>>>
>>>>   Ben Bolker
>>>>
>>>> __
>>>> 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

Re: [Rd] Speeding up R (was Using multicores in R)

[Prof J C Nash (U30A)]

For info, I put a little study I did about the byte code compiler and 
other speedup approaches (but not multicore) on the Rwiki at


http://rwiki.sciviews.org/doku.php?id=tips:rqcasestudy

which looks at a specific problem, so may not be relevant to everyone.
However, one of my reasons for doing it was to document the "how to" a 
little.


JN



   2.  Have you tried the "compiler" package?  If I understand
correctly, R is a two-stage interpreter, first translating what we know
as R into byte code, which is then interpreted by a byte code
interpreter.  If my memory is correct, this approach can cut the compute
time by a factor of 100.





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

[Rd] Lessons from LibreOffice project

[Prof J C Nash (U30A)]

The message below came to me from the Getting Open Source Logic INto 
Government list. I'm passing it on to the devel list as the infoworld 
article may have some ideas of relevance to the R project, mainly 
concerning build and test issues and tracking changes in the code base. 
While the LibreOffice project is very different from R, there may 
nevertheless be some tips we can borrow in the best open source 
tradition. FWIW, the article is quite short.


John Nash

 Original Message 
Subject: [OTT-GOSLING] Interesting article about LibreOffice project
Date: Wed, 06 Mar 2013 04:10:51 +
From: gabriel.cosse...@gmail.com
Reply-To: GOSLING members in Ottawa 


To: ottawa-gosl...@list.goslingcommunity.org

Hi everyone!

Here's a very interesting article about how the LibreOffice project is
evolving:

What you can learn from the monster LibreOffice project
http://www.infoworld.com/print/212908

Have a great day!


Gabriel Cossette
Coordonnateur de communautés sur les logiciels libres | Open Source
Software Communities Coordinator
Technologies de l'information | Information Technology
Services partagés Canada | Parcs Canada
Shared Services Canada | Parks Canada
3, passage du Chien-d'Or, Québec (Québec)  G1R 4V7
gabriel.cosse...@spc-ssc.gc.ca
Tél. | Tel.: 418 649-8163
Cell. | Cell.: 418 254-8558
Gouvernement du Canada | Government of Canada
___
Ottawa-gosling mailing list
ottawa-gosl...@list.goslingcommunity.org
http://list.goslingcommunity.org/mailman/listinfo/ottawa-gosling

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

Re: [Rd] R 3.0, Rtools3.0,l Windows7 64-bit, and permission agony

[Prof J C Nash (U30A)]

While as a Linux user who has not so far been banished to Winland I have 
not experienced this problem, it seems to be the type of issue where a 
"how to", for example, on the R Wiki, would be helpful. Moreover, surely 
this is a name conflict on different platforms, so possibly a list of 
these in a related location would also be useful. I have been struggling 
unsuccessfully to remember one that bit me recently when I did have to 
help a Windows user. It was not a major disaster, but more of a large 
sucking sound of time being wasted figuring out something that was just 
a shade more than trivial.


If there are folk who can get this started, do include me in an off-list 
copy and I'll be willing to test/edit wiki material.


JN

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