Re: [Rd] From .Fortran to .Call?

2020-12-24 Thread Koenker, Roger W
Yes, dotcall64 looks interesting.  There is a paper about it here:

https://www.sciencedirect.com/science/article/pii/S2352711018300785?via%3Dihub

and the R spam package contains many examples of its use for fortran source.

> On Dec 24, 2020, at 12:39 AM, Balasubramanian Narasimhan  
> wrote:
> 
> Also, just came to know about dotcall64::.C64() (on CRAN) which allows for 
> Fortran to be called using .Call().
> 
> -Naras
> 
> On 12/23/20 8:34 AM, Balasubramanian Narasimhan wrote:
>> I think it should be pretty easy to fix up SUtools to use the .Call instead 
>> of .Fortran following along the lines of
>> 
>> https://urldefense.com/v3/__https://github.com/wrathematics/Romp__;!!DZ3fjg!r3_sswU4ZHCe3huoGUy2boX-Vr7aUS-RaExyeh_Rsv8gvGiABcqzvOOKZinG4kC7RtA$
>>  
>> I too deal with a lot of f77 and so I will most likely finish it before the 
>> new year, if not earlier. (Would welcome testers besides myself.)
>> 
>> Incidentally, any idea of what the performance hit is, quantitatively? I 
>> confess I never paid attention to it myself as most Fortran code I use seems 
>> pretty fast, i.e. glmnet.
>> 
>> -Naras
>> 
>> 
>> On 12/23/20 3:57 AM, Koenker, Roger W wrote:
>>> Thanks to all and best wishes for a better 2021.
>>> 
>>> Unfortunately I remain somewhat confused:
>>> 
>>> o  Bill reveals an elegant way to get from my rudimentary registration 
>>> setup to
>>> one that would explicitly type the C interface functions,
>>> 
>>> o Ivan seems to suggest that there would be no performance gain from 
>>> doing this.
>>> 
>>> o  Naras’s pcLasso package does use the explicit C typing, but then 
>>> uses .Fortran
>>> not .Call.
>>> 
>>> o  Avi uses .Call and cites the Romp package 
>>> https://urldefense.com/v3/__https://github.com/wrathematics/Romp__;!!DZ3fjg!r3_sswU4ZHCe3huoGUy2boX-Vr7aUS-RaExyeh_Rsv8gvGiABcqzvOOKZinG4kC7RtA$
>>>  where it is asserted that "there is a (nearly) deprecated interface 
>>> .Fortran() which you
>>> should not use due to its large performance overhead.”
>>> 
>>> As the proverbial naive R (ab)user I’m left wondering:
>>> 
>>> o  if I updated my quantreg_init.c file in accordance with Bill’s 
>>> suggestion could I
>>> then simply change my .Fortran calls to .Call?
>>> 
>>> o  and if so, do I need to include ALL the fortran subroutines in my 
>>> src directory
>>> or only the ones called from R?
>>> 
>>> o  and in either case could I really expect to see a significant 
>>> performance gain?
>>> 
>>> Finally, perhaps I should stipulate that my fortran is strictly f77, so no 
>>> modern features
>>> are in play, indeed most of the code is originally written in ratfor, Brian 
>>> Kernighan’s
>>> dialect from ancient times at Bell Labs.
>>> 
>>> Again,  thanks to all for any advice,
>>> Roger
>>> 
>>> 
 On Dec 23, 2020, at 1:11 AM, Avraham Adler  wrote:
 
 Hello, Ivan.
 
 I used .Call instead of .Fortran in the Delaporte package [1]. What
 helped me out a lot was Drew Schmidt's Romp examples and descriptions
 [2]. If you are more comfortable with the older Fortran interface,
 Tomasz Kalinowski has a package which uses Fortran 2018 more
 efficiently [3]. I haven't tried this last in practice, however.
 
 Hope that helps,
 
 Avi
 
 [1] 
 https://urldefense.com/v3/__https://CRAN.R-project.org/package=Delaporte__;!!DZ3fjg!s1-ihrZ9DPUtXpxdIpJPA1VedpZFt12Ahmn4CycOmile_uSahFZnJPn_5KPITBN5NK8$
 [2] 
 https://urldefense.com/v3/__https://github.com/wrathematics/Romp__;!!DZ3fjg!s1-ihrZ9DPUtXpxdIpJPA1VedpZFt12Ahmn4CycOmile_uSahFZnJPn_5KPISF5aCYs$
 [3] 
 https://urldefense.com/v3/__https://github.com/t-kalinowski/RFI__;!!DZ3fjg!s1-ihrZ9DPUtXpxdIpJPA1VedpZFt12Ahmn4CycOmile_uSahFZnJPn_5KPIbwXmXqY$
 
 Tomasz Kalinowski
 
 
 
 On Tue, Dec 22, 2020 at 7:24 PM Balasubramanian Narasimhan
  wrote:
> To deal with such Fortran issues in several packages I deal with, I
> wrote the SUtools package 
> (https://urldefense.com/v3/__https://github.com/bnaras/SUtools__;!!DZ3fjg!s1-ihrZ9DPUtXpxdIpJPA1VedpZFt12Ahmn4CycOmile_uSahFZnJPn_5KPIJ5BbDwA$
>  ) that you
> can try.  The current version generates the registration assuming
> implicit Fortran naming conventions though. (I've been meaning to
> upgrade it to use the gfortran -fc-prototypes-external flag which should
> be easy; I might just finish that during these holidays.)
> 
> There's a vignette as well:
> https://urldefense.com/v3/__https://bnaras.github.io/SUtools/articles/SUtools.html__;!!DZ3fjg!s1-ihrZ9DPUtXpxdIpJPA1VedpZFt12Ahmn4CycOmile_uSahFZnJPn_5KPITq9-Quc$
>  
> 
> -Naras
> 
> 
> On 12/19/20 9:53 AM, Ivan Krylov wrote:
>> On Sat, 19 Dec 2020 17:04:59 +
>> "Koenker, Roger W"  wrote:
>> 
>>> There are comments in various places, including R-extensions §5.4
>>> suggesting that .Fortran is (nearly) deprecated and hint

Re: [Rd] Silent failure with NA results in fligner.test()

2020-12-24 Thread Martin Maechler
Not sure
If all of the variances are zero,  they are homogenous in that sense,
and I would give a  p-value of 1  ..
if only *some* of the variances are zero... it's less easy.

I still would try to *not* give an error in such cases  and even
prefer  NA  statistic and p-value..  because yes, these are "not
available" for such data.
But it is not strictly an error to try such a test on data of the
correct format...   Consequently, personally I would even try to not
give the current error ... but rather return NA values here:
>  if (all(x == 0))
>  stop("data are essentially constant")

On Mon, Dec 21, 2020 at 12:22 PM Kurt Hornik  wrote:
>
> > Karolis K writes:
>
> Any preferences?
>
> Best
> -k
>
> > Hello,
> > In certain cases fligner.test() returns NaN statistic and NA p-value.
> > The issue happens when, after centering with the median, all absolute 
> > values become constant, which ten leads to identical ranks.
>
> > Below are a few examples:
>
> > # 2 groups, 2 values each
> > # issue is caused by residual values after centering (-0.5, 0.5, -0.5, 0.5)
> > # then, after taking the absolute value, all the ranks become identical.
> >> fligner.test(c(2,3,4,5), gl(2,2))
>
> > Fligner-Killeen test of homogeneity of variances
>
> > data:  c(2, 3, 4, 5) and gl(2, 2)
> > Fligner-Killeen:med chi-squared = NaN, df = 1, p-value = NA
>
>
> > # similar situation with more observations and 3 groups
> >> fligner.test(c(2,3,2,3,4,4,5,5,8,9,9,8), gl(3,4))
>
> > Fligner-Killeen test of homogeneity of variances
>
> > data:  c(2, 3, 2, 3, 4, 4, 5, 5, 8, 9, 9, 8) and gl(3, 4)
> > Fligner-Killeen:med chi-squared = NaN, df = 2, p-value = NA
>
>
> > Two simple patches are proposed below. One returns an error, and another 
> > returns a p-value of 1.
> > Not sure which one is more appropriate, so submitting both.
>
> > Warm regards,
> > Karolis Koncevičius
>
> > ---
>
> > Index: fligner.test.R
> > ===
> > --- fligner.test.R(revision 79650)
> > +++ fligner.test.R(working copy)
> > @@ -59,8 +59,13 @@
> >  stop("data are essentially constant")
>
> >  a <- qnorm((1 + rank(abs(x)) / (n + 1)) / 2)
> > -STATISTIC <- sum(tapply(a, g, "sum")^2 / tapply(a, g, "length"))
> > -STATISTIC <- (STATISTIC - n * mean(a)^2) / var(a)
> > +if (var(a) > 0) {
> > +STATISTIC <- sum(tapply(a, g, "sum")^2 / tapply(a, g, "length"))
> > +STATISTIC <- (STATISTIC - n * mean(a)^2) / var(a)
> > +}
> > +else {
> > +STATISTIC <- 0
> > +}
> >  PARAMETER <- k - 1
> >  PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
> >  names(STATISTIC) <- "Fligner-Killeen:med chi-squared”
>
> > ---
>
> > Index: fligner.test.R
> > ===
> > --- fligner.test.R(revision 79650)
> > +++ fligner.test.R(working copy)
> > @@ -57,6 +57,8 @@
> >  x <- x - tapply(x,g,median)[g]
> >  if (all(x == 0))
> >  stop("data are essentially constant")
> > +if (var(abs(x)) == 0)
> > +stop("absolute residuals from the median are essentially constant")
>
> >  a <- qnorm((1 + rank(abs(x)) / (n + 1)) / 2)
> >  STATISTIC <- sum(tapply(a, g, "sum")^2 / tapply(a, g, "length"))
>
> > __
> > 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
>


-- 
Martinhttp://stat.ethz.ch/~maechler
Seminar für Statistik, ETH Zürich HG G 16   Rämistrasse 101
CH-8092 Zurich, SWITZERLAND   ☎ +41 44 632 3408<><

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