Re: [Rd] Bug in nlm, found using sem; failure in several flavors (PR#13881)
Hi Jeff,
As mentioned in my message, I *did* replicate on another platform.
One platform was running linux in a VM, the other running MacOSX, on two
different physical machines. The problem did not replicate on a PowerPC
box, so if it is a bug with intel chips, I think it's also appropriate to
report here given that the problem seems to cross the OS boundary. I do not
have access to a Windows machine.
Also, I have spent a decent amount of time running memory testers on
the first machine's memory, to ensure that physical memory is not a
problem...the classic "memcheck" program finds no trouble.
The CSV and .ram files both came through as base-64-encoded
attachments, which you actually quoted back to me in your response,
below...but in case your mail client can't read those and/or you are
interested in this bug, I have pasted them at the end of this message (since
they're quite small).
And thanks for the shout-out for bacon. I bought this domain name in
undergrad many years ago, and I must say, I've never regretted it!
--Adam
On Fri, 7 Aug 2009, Jeff Ryan wrote:
Adam,
It seems that your attachment didn't make it through.
That aside, my experience with strange errors like those (random type
not implemented ones) has been that you may be looking at a memory
problem on you machine. Given that you can't replicate on another
platform (and the .csv file didn't come through), I would think it
wise to start there.
My 2c. And I love bacon too :)
Jeff
On Fri, Aug 7, 2009 at 1:10 PM, wrote:
This message is in MIME format. The first part should be readable text,
while the remaining parts are likely unreadable without MIME-aware tools.
--1660387551-1458482416-1249639718=:2997
Content-Type: TEXT/PLAIN; CHARSET=US-ASCII; FORMAT=flowed
Content-ID:
Hello,
There appears to be a bug in the nlm function, which I discovered
when trying to run an SEM. The SEM is not bizarre; the covariance matrix is
solve-able and the eigenvalues are greater than zero, and I do not believe
the "sem" package per se to be the problem (as nlm keeps being the part that
fails, though I can't replicate this with other nlm tasks). I apologize if
I have put too many much information in this message; I'm not a programmer
by trade so I don't really know what's going on here, which hampers my
ability to write consise bug reports.
Here is the code I use:
library(sem)
ice.S <- read.csv("iceS.csv") # attached
rownames(ice.S) <- ice.S[,1]
ice.S[,1] <- NULL
ice.S <- as.matrix(ice.S)
ice.ram <- specify.model("ice.ram") # attached
ice.N <- 342
ice.sem <- sem(ram=ice.ram, S=ice.S, N=ice.N)
...at this point, any number of problems could occur. I have seen the
following.
1) Simple lack of convergence. (might be my model's fault.)
2) Error in nlm(if (analytic.gradient) objective.2 else objective.1, start,
:
type 31 is unimplemented in 'type2char'
3) *** caught segfault ***
address 0xc019c87b, cause 'memory not mapped'
Traceback:
1: nlm(if (analytic.gradient) objective.2 else objective.1, start,
hessian = TRUE, iterlim = maxiter, print.level = if (debug) 2 else 0,
typsize = typsize, ...)
2: sem.default(ram = ram, S = S, N = N, param.names = pars, var.names =
vars, fixed.x = fixed.x, debug = debug, ...)
3: sem(ram = ram, S = S, N = N, param.names = pars, var.names = vars,
fixed.x = fixed.x, debug = debug, ...)
4: sem.mod(ram = ice.ram, S = ice.S, N = ice.N)
5: sem(ram = ice.ram, S = ice.S, N = ice.N)
Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection: 1
aborting ...
Segmentation fault
swiss:data$
(no core was dumped).
Trying with debug mode provides other interesting errors:
ice.sem <- sem(ram=ice.ram, S=ice.S, N=ice.N, debug=TRUE)
...gets up to some iteration (not always 15), and then R dies ungracefully,
and exits to the shell:
iteration = 15
Step:
[1] 1.253132e-02 1.183343e-02 -7.651342e-03 -2.490800e-03 2.278938e-03
[6] 3.197431e-04 6.137849e-04 -2.496882e-03 -1.065829e-03 -2.118179e-03
[11] 2.942936e-03 -1.335936e-03 -3.665618e-03 3.090566e-03 8.534956e-04
[16] -1.440421e-03 -5.230877e-04 -1.053686e-04 -9.771005e-04 -4.269216e-04
[21] 7.261694e-05 -1.039668e-03 -8.409151e-03 -3.497456e-03 2.216899e-03
[26] -4.633520e-03 -4.052130e-03 -4.746537e-03 -1.589622e-03 -2.523766e-04
Parameter:
[1] -0.76604614 -1.19639662 0.83456888 0.72540413 0.08482452 0.56180393
[7] 0.50615814 0.55728015 0.83796696 0.88371335 -0.70465116 0.85251098
[13] -0.18346956 0.66857787 0.57012481 0.39116561 0.91237990 0.63951482
[19] 0.26601566 0.29240836 0.44710919 0.94734056 6.52039015 0.02524762
[25] -0.01614603 2.88198219 0.03442452 3.52272237 1.44698423 -0.72964745
Function Value
[1] -15175.94
Gradient:
[1] -2.085412e+07 -3.819717e+07 3.883989e+07 1.352594e+00 -4.283329e+00
[6] -1.437250e+00 -6.558913e-01 1.358276e+00 7.930865e+05 -1.293853e+06
[11] -5.816627e+03 5.90819
Re: [Rd] Bug in nlm, found using sem; failure in several flavors (PR#13883)
This message is in MIME format. The first part should be readable text,
while the remaining parts are likely unreadable without MIME-aware tools.
--1660387551-150661043-1249684349=:2997
Content-Type: TEXT/PLAIN; charset=iso-8859-1; format=flowed
Content-Transfer-Encoding: QUOTED-PRINTABLE
Hi Jeff,
=09As mentioned in my message, I *did* replicate on another platform.=20
One platform was running linux in a VM, the other running MacOSX, on two
different physical machines. The problem did not replicate on a PowerPC
box, so if it is a bug with intel chips, I think it's also appropriate to
report here given that the problem seems to cross the OS boundary. I do no=
t
have access to a Windows machine.
=09Also, I have spent a decent amount of time running memory testers on
the first machine's memory, to ensure that physical memory is not a
problem...the classic "memcheck" program finds no trouble.
=09The CSV and .ram files both came through as base-64-encoded
attachments, which you actually quoted back to me in your response,
below...but in case your mail client can't read those and/or you are
interested in this bug, I have pasted them at the end of this message (sinc=
e
they're quite small).
=09And thanks for the shout-out for bacon. I bought this domain name in
undergrad many years ago, and I must say, I've never regretted it!
--Adam
On Fri, 7 Aug 2009, Jeff Ryan wrote:
> Adam,
>
> It seems that your attachment didn't make it through.
>
> That aside, my experience with strange errors like those (random type
> not implemented ones) has been that you may be looking at a memory
> problem on you machine. Given that you can't replicate on another
> platform (and the .csv file didn't come through), I would think it
> wise to start there.
>
> My 2c. And I love bacon too :)
>
> Jeff
>
> On Fri, Aug 7, 2009 at 1:10 PM, wrote:
>> =A0This message is in MIME format. =A0The first part should be readable =
text,
>> =A0while the remaining parts are likely unreadable without MIME-aware to=
ols.
>>
>> --1660387551-1458482416-1249639718=3D:2997
>> Content-Type: TEXT/PLAIN; CHARSET=3DUS-ASCII; FORMAT=3Dflowed
>> Content-ID:
>>
>> Hello,
>>
>> =A0 =A0 =A0 =A0There appears to be a bug in the nlm function, which I di=
scovered
>> when trying to run an SEM. =A0The SEM is not bizarre; the covariance mat=
rix is
>> solve-able and the eigenvalues are greater than zero, and I do not belie=
ve
>> the "sem" package per se to be the problem (as nlm keeps being the part =
that
>> fails, though I can't replicate this with other nlm tasks). =A0I apologi=
ze if
>> I have put too many much information in this message; I'm not a programm=
er
>> by trade so I don't really know what's going on here, which hampers my
>> ability to write consise bug reports.
>>
>> Here is the code I use:
>>
>> library(sem)
>> ice.S <- read.csv("iceS.csv") # attached
>> rownames(ice.S) <- ice.S[,1]
>> ice.S[,1] <- NULL
>> ice.S <- as.matrix(ice.S)
>> ice.ram <- specify.model("ice.ram") # attached
>> ice.N <- 342
>> ice.sem <- sem(ram=3Dice.ram, S=3Dice.S, N=3Dice.N)
>>
>> ...at this point, any number of problems could occur. I have seen the
>> following.
>>
>> 1) Simple lack of convergence. (might be my model's fault.)
>> 2) Error in nlm(if (analytic.gradient) objective.2 else objective.1, sta=
rt,
>> :
>> =A0 type 31 is unimplemented in 'type2char'
>> 3) =A0*** caught segfault ***
>> address 0xc019c87b, cause 'memory not mapped'
>>
>> Traceback:
>> =A01: nlm(if (analytic.gradient) objective.2 else objective.1, start,
>> hessian =3D TRUE, iterlim =3D maxiter, print.level =3D if (debug) 2 else=
0,
>> typsize =3D typsize, ...)
>> =A02: sem.default(ram =3D ram, S =3D S, N =3D N, param.names =3D pars, v=
ar.names =3D
>> vars, =A0 =A0 fixed.x =3D fixed.x, debug =3D debug, ...)
>> =A03: sem(ram =3D ram, S =3D S, N =3D N, param.names =3D pars, var.names=
=3D vars,
>> fixed.x =3D fixed.x, debug =3D debug, ...)
>> =A04: sem.mod(ram =3D ice.ram, S =3D ice.S, N =3D ice.N)
>> =A05: sem(ram =3D ice.ram, S =3D ice.S, N =3D ice.N)
>>
>> Possible actions:
>> 1: abort (with core dump, if enabled)
>> 2: normal R exit
>> 3: exit R without saving workspace
>> 4: exit R saving workspace
>> Selection: 1
>> aborting ...
>> Segmentation fault
>> swiss:data$
>>
>> (no core was dumped).
>>
>> Trying with debug mode provides other interesting errors:
>>
>>> ice.sem <- sem(ram=3Dice.ram, S=3Dice.S, N=3Dice.N, debug=3DTRUE)
>>
>> ...gets up to some iteration (not always 15), and then R dies ungraceful=
ly,
>> and exits to the shell:
>>
>> iteration =3D 15
>> Step:
>> =A0[1] =A01.253132e-02 =A01.183343e-02 -7.651342e-03 -2.490800e-03 =A02.=
278938e-03
>> =A0[6] =A03.197431e-04 =A06.137849e-04 -2.496882e-03 -1.065829e-03 -2.11=
8179e-03
>> [11] =A02.942936e-03 -1.335936e-03 -3.665618e-03 =A03.090566e-03 =A08.53=
4956e-04
>> [16] -1.440421e-03 -5.230877e-04 -1.053686e-04 -9.771005e-04 -4.269216e-=
04
>> [21] =A07.261694e-05 -1.039668e-03 -8.409151e-03 -3.497456e-03 =A02.2168=
Re: [Rd] identical(0, -0)
> "DM" == Duncan Murdoch
> on Fri, 07 Aug 2009 12:55:50 -0400 writes:
DM> On 8/7/2009 11:41 AM, Martin Maechler wrote:
>>> "DM" == Duncan Murdoch
>>> on Fri, 07 Aug 2009 11:25:11 -0400 writes:
>>
DM> On 8/7/2009 10:46 AM, Martin Maechler wrote:
>> >>> "TH" == Ted Harding
>> >>> on Fri, 07 Aug 2009 14:49:54 +0100 (BST) writes:
>> >>
TH> On 07-Aug-09 11:07:08, Duncan Murdoch wrote:
>> >> >> Martin Maechler wrote:
>> >> William Dunlap
>> >> on Thu, 6 Aug 2009 15:06:08 -0700 writes:
>> >> >>> >> -Original Message- From:
>> >> >>> >> r-help-boun...@r-project.org
>> >> >>> >> [mailto:r-help-boun...@r-project.org] On Behalf Of
>> >> >>> >> Giovanni Petris Sent: Thursday, August 06, 2009 3:00 PM
>> >> >>> >> To: milton.ru...@gmail.com Cc: r-h...@r-project.org;
>> >> >>> >> daniel.gerl...@geodecapital.com Subject: Re: [R] Why is 0
>> >> >>> >> not an integer?
>> >> >>> >>
>> >> >>> >>
>> >> >>> >> I ran an instant experiment...
>> >> >>> >>
>> >> >>> >> > typeof(0) [1] "double" > typeof(-0) [1] "double" >
>> >> >>> >> identical(0, -0) [1] TRUE
>> >> >>> >>
>> >> >>> >> Best, Giovanni
>> >> >>>
>> >> >>> > But 0.0 and -0.0 have different reciprocals
>> >> >>>
>> >> >>> >> 1.0/0.0
>> >> >>> >[1] Inf
>> >> >>> >> 1.0/-0.0
>> >> >>> >[1] -Inf
>> >> >>>
>> >> >>> > Bill Dunlap TIBCO Software Inc - Spotfire Division wdunlap
>> >> >>> > tibco.com
>> >> >>>
>> >> >>> yes. {finally something interesting in this boring thread !}
---> diverting to R-devel
>> >> >>>
>> >> >>> In April, I've had a private e-mail communication with John
>> >> >>> Chambers [father of S, notably S4, which also brought identical()]
>> >> >>> and Bill about the topic,
>> >> >>> where I had started suggesting that R should be changed such
>> >> >>> that
>> >> >>> identical(-0. , +0.)
>> >> >>> would return FALSE.
>> >> >>> Bill did mention that it does so for (newish versions of) S+
>> >> >>> and that he'd prefer that, too,
>> >> >>> and John said
>> >> >>>
>> >> >>> >> I agree on having a preference for a bitwise comparison for
>> >> >>> >> identical()---that's what the name means after all. But since
>> >> >>> >> someone implemented the numerical case as the C == it's
probably
>> >> >>> >> going to be more hassle than it's worth to change it. But we
>> >> >>> >> should make the implementation clear in the documentation.
>> >> >>>
>> >> >>> so in principle, we all agreed that R's identical() should be
>> >> >>> changed here, namely by using something like memcmp() instead
>> >> >>> of simple '==' , however we haven't bothered to actually
>> >> >>> *implement* this change.
>> >> >>>
>> >> >>> I am currently testing a patch which would lead to
>> >> >>> identical(0, -0) return FALSE.
>> >> >>>
>> >> >> I don't think that would be a good idea. Other expressions besides
>> >> >> "-0"
>> >> >> calculate the zero with the negative sign bit, e.g. the following
>> >> >> sequence:
>> >> >>
>> >> >> pos <- 1
>> >> >> neg <- -1
>> >> >> zero <- 0
>> >> >> y <- zero*pos
>> >> >> z <- zero*neg
>> >> >> identical(y, z)
>> >> >>
>> >> >> I think most R users would expect the last expression there to be
>> >> >> TRUE based on the previous two lines, given that pos and neg both
>> >> >> have finite values. In a simple case like this y == z would be a
>> >> >> better test to use, but if those were components of a larger
>> >> >> structure, identical() is all we've got, and people would waste a
>> >> >> lot of time tracking down why structures differing only in the
>> >> >> sign of zero were not identical, even though every element tested
>> >> >> equal.
>> >>
>> >> identical() *is* not the same as '==' even if you think of a
>> >> generalized '==',
>> >> and your example is not convincing to me.
>>
DM> Fair enough, but after your change, how would one do what
DM> identical(list(pos, neg, zero, y), list(pos, neg, zero, z)) does now?
DM> That seems to me to be a more useful comparison than one that declares
DM> those to be unequal because the signs of y and z differ.
>>
>> Maybe something like
>>
>> all(mapply(`==`, list(pos, neg, zero, y), list(pos, neg, zero, z)))
>>
>> ## or even
>>
>> isTRUE(all.equal( list(pos, neg, zero, y), list(pos, neg, zero, z),
>> tol = 0))
DM> I think I didn't make my point clearly. I'm not particularly worried
DM> about lists of numbers, I'm worried about signed zeros buried in
complex
DM> structures. identical(struc1, struc2) works nicely now for that sort
of
DM> comparison; indeed the man page for i
Re: [Rd] identical(0, -0)
I'll save space and not include previous messages. My 2 cents: At the very least the documentation needs a fix. If it is easy to do, then Ted Harding's suggestion of a switch (default OFF) to check for sign difference would be sensible. I would urge inclusion in the documentation of the +0, -0 example(s) if there is NOT a way in R to distinguish these. There are occasions where it is useful to be able to detect things like this (and NaN and Inf and -Inf etc.). They are usually not of interest to users, but sometimes are needed for developers to check edge effects. For those cases it may be time to consider a package FPIEEE754 or some similar name to allow testing and possibly setting of flags for some of the fancier features. Likely used by just a few of us in extreme situations. Unfortunately, some of the nice exception handling that was suggested in the standard is apparently rarely implemented in compilers. For info, I was actually a voting member of IEEE 754 because I found a nice "feature" in the IBM Fortran G arithmetic, though I never sat in on any of the meetings. JN __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] identical(0, -0)
On Sat, Aug 8, 2009 at 10:39 AM, Prof. John C Nash wrote: > I would urge inclusion in the documentation of the +0, -0 example(s) if > there is NOT a way in R to distinguish these. There are occasions where it For single numbers try this: x <- +0 y <- -0 identical(x, y) && identical(1/x, 1/y) # FALSE __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Problem using model.frame with argument subset in own function
Dear List,
I am writing a formula method for a function in a package I maintain. I
want the method to return a data.frame that potentially only contains
some of the variables in 'data', as specified by the formula.
The problem I am having is in writing the function and wrapping it
around model.frame. Consider the following data frame:
dat <- data.frame(A = runif(10), B = runif(10), C = runif(10))
And the wrapper function:
foo <- function(formula, data = NULL, ..., subset = NULL,
na.action = na.pass) {
mt <- terms(formula, data = data, simplify = TRUE)
mf <- model.frame(formula(mt), data = data, subset = subset,
na.action = na.action)
## real function would do more stuff here and pass mf on to
## other functions
mf
}
This is how I envisage the function being called. The real world use
would have a data.frame with tens or hundreds of components where only a
few need to be excluded. Hence wanting formulas of the form below to
work.
foo(~ . - B, data = dat)
The aim is to return only columns A and C in an object returned by
model.frame. However, when I run the above, I get the following error:
> foo(~ A + B, data = dat)
Error in xj[i] : invalid subscript type 'closure'
I've tracked this down to the line in model.frame.default
subset <- eval(substitute(subset), data, env)
After evaluating this line, subset contains:
Browse[1]> subset
function (x, ...)
UseMethod("subset")
Not NULL, and hence the error later on when calling the internal
model.frame code.
So the question is, what am I doing wrong?
If I leave the subset argument out of the definition of foo and rely
upon the default in model.frame.default, the function works as
expected.
Perhaps the question should be, how do I modify foo() to allow it to
have a formal subset argument, passed to model.frame?
Any other suggestions gratefully accepted.
Thanks in advance,
G
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
