Re: [Rd] Replace isnan and lgamma in Fortran subroutine in R package

2014-09-23 Thread Martyn Plummer
On Tue, 2014-09-23 at 07:43 +0200, Berend Hasselman wrote:
> On 23-09-2014, at 00:33, Wang, Zhu  wrote:
> 
> > Hello,
> > 
> > I submitted a package which used Fortran functions isnan and lgamma. 
> > However, I was told that:
> > 
> > isnan and lgamma are not Fortran 95 functions.
> > 
> > I was asked to write 'cross-platform portable code' and so should not be 
> > writing GNU extensions to Fortran.
> > 
> > See http://cran.r-project.org/web/checks/check_results_mpath.html, which 
> > will shortly show installation failures under Solaris.
> > 
> > I will appreciate advice on how to replace these two functions to avoid 
> > failure on some platforms.
> > 
> 
> I don’t know about lgamma.
> 
> Instead of isnan you could use Lapack’s logical function disnan to test for 
> NaN (it’s in lapack 3.4.2; I don’t know about earlier versions).
> 
> Another way would be to write a C function using functions provided by R to 
> test for NaN. 
> That function should be declared with F77_NAME so that it can be called from 
> a Fortran routine.
> 
> I haven’t tried this so I don’t know if this would be foolproof or is the 
> best way to do it..
> 
> Berend
> 

As described by Berend (See also section 6.6 of Writing R Extensions)
you can define Fortran-callable wrappers around ISNAN and lgammafn, both
of which are provided with the C interface to R (lgammafn is declared in
Rmath.h).  Then declare these wrappers as external subroutines in your
Fortran code and call them in place of isnan and lgamma.

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


Re: [Rd] Replace isnan and lgamma in Fortran subroutine in R package

2014-09-23 Thread Martyn Plummer
Try this patch.
Martyn

On Mon, 2014-09-22 at 22:33 +, Wang, Zhu wrote:
> Hello,
> 
> I submitted a package which used Fortran functions isnan and lgamma. However, 
> I was told that:
> 
> isnan and lgamma are not Fortran 95 functions.
> 
> I was asked to write 'cross-platform portable code' and so should not be 
> writing GNU extensions to Fortran.
> 
> See http://cran.r-project.org/web/checks/check_results_mpath.html, which will 
> shortly show installation failures under Solaris.
> 
> I will appreciate advice on how to replace these two functions to avoid 
> failure on some platforms.
> 
> Thanks in advance,
> 
> Zhu Wang
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel


---
This message and its attachments are strictly confidential. If you are
not the intended recipient of this message, please immediately notify
the sender and delete it. Since its integrity cannot be guaranteed,
its content cannot involve the sender's responsibility. Any misuse,
any disclosure or publication of its content, either whole or partial,
is prohibited, exception made of formally approved use
---
diff -uNr mpath/src/cfuns.c mpath-fixed/src/cfuns.c
--- mpath/src/cfuns.c	1970-01-01 01:00:00.0 +0100
+++ mpath-fixed/src/cfuns.c	2014-09-23 09:56:24.303347209 +0200
@@ -0,0 +1,5 @@
+#include 
+#include 
+
+int F77_SUB(cisnan)(double *x) { return ISNAN(*x); }
+double F77_SUB(rlgamma)(double *x) { return lgammafn(*x); }
diff -uNr mpath/src/midloop.f mpath-fixed/src/midloop.f
--- mpath/src/midloop.f	2014-09-22 19:32:25.0 +0200
+++ mpath-fixed/src/midloop.f	2014-09-23 09:56:55.455497773 +0200
@@ -255,6 +255,8 @@
 
   integer n, family, i
   double precision dev, y(n), mu(n), theta, weights(n),tmp
+  integer cisnan
+  external cisnan
  
   dev = 0.0D0
   do 10 i=1, n
@@ -282,7 +284,7 @@
 dev=dev+2*(weights(i)*(y(i)*dlog(max(1.0D0,y(i))/mu(i))-
  +  (y(i)+theta)*dlog((y(i)+theta)/(mu(i) + theta 
   endif
-  if(isnan(dev)) then
+  if(cisnan(dev).NE.0) then
   call intpr("i=", -1, i, 1)
   call dblepr("y(i)=", -1, y(i), 1)
   call dblepr("mu(i)=", -1, mu(i), 1)
@@ -300,6 +302,8 @@
 
   integer n, family, i, y0
   double precision ll, y(n), mu(n), theta, w(n)
+  double precision rlgamma
+  external rlgamma
 
   ll = 0
   do 10 i=1,n
@@ -309,8 +313,8 @@
   else 
   y0=0
   endif
-   ll=ll + w(i) * (lgamma(theta+y(i))-lgamma(theta)-
- +   lgamma(y(i)+1)+
+   ll=ll + w(i) * (rlgamma(theta+y(i))-rlgamma(theta)-
+ +   rlgamma(y(i)+1)+
  +   theta*log(theta) + y(i)*log(mu(i)+y0) - (theta + y(i))*
  +   log(theta +  mu(i)))
   else if(family .EQ. 1)then !gaussian
@@ -320,7 +324,7 @@
 ll=ll + w(i)*(y(i)*log(mu(i)/(1-mu(i)))+log(1-mu(i)))
 endif
 else if(family .EQ. 3) then !poisson
-ll=ll + w(i)*(-mu(i) + y(i)*log(mu(i))-lgamma(y(i)+1))
+ll=ll + w(i)*(-mu(i) + y(i)*log(mu(i))-rlgamma(y(i)+1))
   endif
  10   continue
   return
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


Re: [Rd] Replace isnan and lgamma in Fortran subroutine in R package

2014-09-23 Thread Wang, Zhu
That patch works flawlessly.

Thanks Martyn.

Zhu

-Original Message-
From: Martyn Plummer [mailto:plumm...@iarc.fr] 
Sent: Tuesday, September 23, 2014 8:39 AM
To: Wang, Zhu
Cc: r-devel@r-project.org
Subject: Re: [Rd] Replace isnan and lgamma in Fortran subroutine in R package

Try this patch.
Martyn

On Mon, 2014-09-22 at 22:33 +, Wang, Zhu wrote:
> Hello,
> 
> I submitted a package which used Fortran functions isnan and lgamma. However, 
> I was told that:
> 
> isnan and lgamma are not Fortran 95 functions.
> 
> I was asked to write 'cross-platform portable code' and so should not be 
> writing GNU extensions to Fortran.
> 
> See http://cran.r-project.org/web/checks/check_results_mpath.html, which will 
> shortly show installation failures under Solaris.
> 
> I will appreciate advice on how to replace these two functions to avoid 
> failure on some platforms.
> 
> Thanks in advance,
> 
> Zhu Wang
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel


---
This message and its attachments are strictly confidenti...{{dropped:8}}

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


Re: [Rd] Bug in new behaviour for all.equal and environments?

2014-09-23 Thread Martin Maechler
> Duncan Murdoch 
> on Sun, 21 Sep 2014 15:57:00 -0400 writes:

> On 21/09/2014, 1:38 PM, Rui Barradas wrote:
>> Hello,
>> 
>> In R 3.1.1 on Windows 7 it's ok.
>> 
>> > all.equal(baseenv(), baseenv()) [1] TRUE >
>> sessionInfo() R version 3.1.1 (2014-07-10) Platform:
>> x86_64-w64-mingw32/x64 (64-bit)

> I'm not sure if that's really "ok", since it returns TRUE
> for any pair of environments, e.g.

> all.equal(baseenv(), emptyenv())

> The current R-devel behaviour is temporary and certain to
> change.

> Duncan Murdoch

Indeed.  A bit more than an hour ago, I have committed several
changes to R-devel's behavior in dealing with objects such as 
 
 (rp <- getClass("refClass")@prototype)
 (ner <- new("envRefClass"))

both of which would not even correctly *print* in current R.
They now do in R-devel, also work with str() there, and

  all.equal(baseenv(), baseenv())

and many more cases of all.equal() with environments now work
there as well.

I still expect that creative useRs will be able to construct
cases where the new  all.equal() methods will fail, possibly
even end in an infinite loop (which is caught and leads to an
error).

Please use R-devel  svn rev >= 8  for such experiments.

Thanking you for testing ...,
Martin Maechler
ETH Zurich



>> []

>> Rui Barradas
>> 
>> Em 21-09-2014 18:06, Kevin Ushey escreveu:
>>> Hi R-devel,
>>> 
>>> The following code:
>>> 
>>> all.equal(baseenv(), baseenv())
>>> 
>>> gives the error when run in a clean R session with
>>> latest R-devel (r66650):
>>> 
>>> kevin:~$ R --vanilla --slave -e "all.equal(baseenv(),
>>> baseenv())" Error in all.equal.envRefClass(target[[i]],
>>> current[[i]], check.attributes = check.attributes, :
>>> attempt to apply non-function Calls: all.equal
>>> ... all.equal.list -> all.equal -> all.equal.envRefClass
>>> Execution halted
>>> 
>>> Although I don't know if it's helpful -- it appears that
>>> packages that include some S4 machinery will effect the
>>> outcome of this error, e.g.  if we load 'Biobase' first:
>>> 
>>> kevin:~$ R --vanilla --slave -e
>>> "suppressPackageStartupMessages(library(Biobase));
>>> all.equal(baseenv(), baseenv())" Error in
>>> target$getClass() : object '.refClassDef' not found
>>> Calls: all.equal ... all.equal.list -> all.equal ->
>>> all.equal.envRefClass ->  Execution halted
>>> 
>>> We were bumping into an error with this in Rcpp -- we
>>> used all.equal (through RUnit) to confirm that
>>> baseenv(), and a function returning the base
>>> environment, would return TRUE.
>>> 
>>> For completeness:
>>> 
>>> kevin:~$ R --vanilla --slave -e "sessionInfo()" R Under
>>> development (unstable) (2014-09-20 r66650) Platform:
>>> x86_64-apple-darwin13.3.0 (64-bit)
>>> 
>>> locale: [1]
>>> en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
>>> 
>>> attached base packages: [1] stats graphics grDevices
>>> utils datasets methods base
>>> 
>>> Thanks, Kevin

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


[Rd] Patch for R to fix some buffer overruns and add a missing PROTECT().

2014-09-23 Thread Karl Millar
This patch is against current svn and contains three classes of fix:
   - Ensure the result is properly terminated after calls to strncpy()
   - Replace calls of sprintf() with snprintf()
   - Added a PROTECT() call in do_while which could cause memory
errors if evaluating the condition results in a warning.

Thanks,

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


Re: [Rd] Patch for R to fix some buffer overruns and add a missing PROTECT().

2014-09-23 Thread Duncan Murdoch

On 23/09/2014 3:20 PM, Karl Millar wrote:

This patch is against current svn and contains three classes of fix:
- Ensure the result is properly terminated after calls to strncpy()
- Replace calls of sprintf() with snprintf()
- Added a PROTECT() call in do_while which could cause memory
errors if evaluating the condition results in a warning.


Nothing was attached.

Generally fixes like this are best sent to bugs.r-project.org, and they 
receive highest priority if accompanied by code demonstrating why they 
are needed, i.e. crashes or incorrect results in current R.  Those will 
likely be incorporated as regression tests.


Duncan Murdoch

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


Re: [Rd] Patch for R to fix some buffer overruns and add a missing PROTECT().

2014-09-23 Thread Karl Millar
Bug submitted.  Thanks.

On Tue, Sep 23, 2014 at 12:42 PM, Duncan Murdoch
 wrote:
> On 23/09/2014 3:20 PM, Karl Millar wrote:
>>
>> This patch is against current svn and contains three classes of fix:
>> - Ensure the result is properly terminated after calls to strncpy()
>> - Replace calls of sprintf() with snprintf()
>> - Added a PROTECT() call in do_while which could cause memory
>> errors if evaluating the condition results in a warning.
>
>
> Nothing was attached.
>
> Generally fixes like this are best sent to bugs.r-project.org, and they
> receive highest priority if accompanied by code demonstrating why they are
> needed, i.e. crashes or incorrect results in current R.  Those will likely
> be incorporated as regression tests.
>
> Duncan Murdoch

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


Re: [Rd] Pointer ownership with GECreateDevDesc/GEDestroyDevDesc

2014-09-23 Thread Paul Murrell


This seems like the best solution (free and NULL the dev pointer in 
dev_Close(), BEFORE GEdestroyDevDesc() has a go).


I think this should be safe because dev_Close() and GEdestroyDevDesc() 
are only called in one place (so your free() should always happen before 
GEdestroyDevDesc() and your free() shouldn't happen in any other, 
potentially inappropriate, context).


Paul

On 09/23/14 03:15, Ritch Melton wrote:

I perused the source for RStudio and noticed that they had a similar
problem. I didn't think that the GEDevDesc structure was exported, but
apparently this is a worthwhile workaround to the defect.

void GD_Close(pDevDesc dev)
{
...elided...
   // explicitly free and then null out the dev pointer of the GEDevDesc
   // This is to avoid incompatabilities between the heap we are
compiled with
   // and the heap R is compiled with (we observed this to a problem with
   // 64-bit R)
   std::free(s_pGEDevDesc->dev);
   s_pGEDevDesc->dev = NULL;

   // set GDDevDesc to NULL so we don't reference it again
   s_pGEDevDesc = NULL;
}
}

I will

On Fri, Sep 19, 2014 at 10:37 AM, Ritch Melton  wrote:


According to the "R Internals" document, for a custom device, I should
create a pDevDesc structure that gets passed to GECreateDevDesc.

..elided...
pDevDesc dev;
/* Allocate and initialize the device driver data */
if (!(dev = (pDevDesc) calloc(1, sizeof(DevDesc return 0;
/* or error() */
/* set up device driver or free ’dev’ and error() */
gdd = GEcreateDevDesc(dev); GEaddDevice2(gdd, "dev_name");
...elided...

which indicates to me that the calling code owns the pDevDesc structure
and should be responsible for freeing it.
However, in GEdestroyDevDesc, that particular code calls free(dd->dev),
freeing the pointer.

In my case, I have a device implemented and created in a language that is
using a different allocator. The pDevDesc I pass in is not owned by R and
causes a crash on free with "dev.off()".

I thought I could mitigate the issue, by making a call to the allocator R
is using via an exported R function, such as R_alloc, but I did not see
such a function available to me.

I would like to know if I should be creating the pDevDesc via an imported
alloc routine somehow? or is there a way to prevent this code from freeing
the structure I pass in.

Blue Skies,
Ritch



[[alternative HTML version deleted]]

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



--
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/

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