chisqsim broken on alpha - solved + next bug in fft found

[Albrecht Gebhardt]

On Wed, 15 Sep 1999 ripley@stats.ox.ac.uk wrote:

> On Tue, 14 Sep 1999, Albrecht Gebhardt wrote:
> 
> [...]
> 
> > chisqsim (called from chisq.test() crashes. First I found that the crash
> > is at 
> > chisqsim.c:167
> > 167             fact[i] = x;
> 
> [...]
> 
> chisqsum.c has its interface in terms of long *, called by .C that supplies
> int *. I think those are different on an Alpha.  Could you please
> edit long to int in that file and try again.
> 
> Brian
> 

Yes this helps.
The patch is:
--- ./src/appl/chisqsim.c.chisqsim-patch        Wed Sep 15 09:17:56 1999
+++ ./src/appl/chisqsim.c       Wed Sep 15 14:30:56 1999
@@ -11,10 +11,15 @@
 #include "S.h"
 
 static void
+#ifdef __osf__
+rcont2(int *nrow, int *ncol, int *nrowt, int *ncolt, int *ntotal,
+       double *fact, int *jwork, int *matrix)
+#else
 rcont2(long *nrow, long *ncol, long *nrowt, long *ncolt, long *ntotal,
        double *fact, long *jwork, long *matrix) 
+#endif
 {
-    long nlmp, j, l, m, ia, ib, ic, jc, id, ie, ii, nrowtl, iap, idp,
+    int nlmp, j, l, m, ia, ib, ic, jc, id, ie, ii, nrowtl, iap, idp,
        igp, ihp, iip, nll, nlm, nrowm, ncolm, lsm, lsp;
     double x, y, dummy, sumprb;
 
@@ -151,12 +156,18 @@
    */
 
 void
+#ifdef __osf__
+chisqsim(int *nrow, int *ncol, int *nrowt, int *ncolt, int *n,
+        int *b, double *expected, int *observed, double *fact,
+        int *jwork, double *results)
+#else
 chisqsim(long *nrow, long *ncol, long *nrowt, long *ncolt, long *n,
         long *b, double *expected, long *observed, double *fact,
         long *jwork, double *results)
+#endif
 {
     /* Local variables */
-    long i, j, iter;
+    int i, j, iter;
     double chi, e, o, x;
 
     /* Calculate log-factorials */


And I get to the next bug:

> data(faithful)
> x <- xx <- faithful$eruptions
> x[i.out <- sample(length(x), 10)] <- NA
  I modifed this to 
  x[i.out <- c(1,10,20,50,100,200)] <- NA
> doR <- density(x, bw=0.15, na.rm = TRUE) # works
> doN <- density(x, bw=0.15, na.rm = FALSE)
segfault at approx.c: 

    for(i=0 ; i<*nxy ; i++)
	if(ISNA(x[i]) || ISNA(y[i]))
	    error("approx(): attempted to interpolate NA values\n");

Both x and y are NULL pointers. 
The error is somewhere before, comparing with R on SuSE 5.3:
somewhere in density():

debug: kords[(n + 2):(2 * n)] <- -kords[n:2]
linux and osf both get:
kords
[1]  0.00000000  0.01094819  0.02189638  0.03284457  0.04379277 ...
[1021] -0.04379277 -0.03284457 -0.02189638 -0.01094819

then in the following if-statement we go to
debug: kords <- convolve(y, kords)[1:n]
with linux result
kords
  [1] 2.523589e-13 4.199167e-13 6.950994e-13 1.144613e-12 ...
[511] 5.032797e-13 3.017881e-13
and alpha:
kords
  [1] NA NA NA NA NA NA NA NA ...
[501] NA NA NA NA NA NA NA NA NA NA NA NA
... now on the alpha x and y are set to numeric(0) and approx is called
with this arguments, bang.

Now in convolve(x,y) on entry x=y and y=kords on linux and osf are
identical.
x:
   [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00  ...
 [151] 4.866728e-03 6.204044e-03 3.258272e-03 3.552390e-03 4.218750e-03
[1021] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00  ...
y:
   [1]  2.659615e+00  2.652540e+00  2.631429e+00  2.596616e+00 ...
[1021]  2.548649e+00  2.596616e+00  2.631429e+00  2.652540e+00 

now the error is in line
debug: x <- fft(fft(x) * (if (conj) Conj(fft(y)) else fft(y)), inv = TRUE)
of convolve() which returns only NAs on the alpha.

... so tomorrow I will have a look into fft.c



Albrecht

-------------------------------------------------------------------------------
Albrecht Gebhardt                   email   : albrecht.gebhardt@uni-klu.ac.at
Institut fuer Mathematik            Tel.    : (++43 463) 2700/837
Universitaet Klagenfurt             Fax     : (++43 463) 2700/834
Villacher Str. 161
A-9020 Klagenfurt, Austria
-------------------------------------------------------------------------------


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._