bugPSPP - Bugs: bug #25466, Compile problem on wilcoxon.c

 
 

bug #25466: Compile problem on wilcoxon.c

Submitter:  Michel Boaventura <michelboaventura>
Submitted:  Mon 02 Feb 2009 06:35:48 PM UTC
   
 
Category:  Compilation/Portability Severity:  5 - Average
Status:  Fixed Assigned to:  None
Open/Closed:  Closed Release:  None
Effort:  0.00
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 06 Feb 2009 06:28:32 AM UTC, comment #17: 


>You need also to remove the references to SIGALRM in main.c and
>psppire.c


Thanks for noticing!

With that change, I pushed these changes, and so I'm closing the bug.  Michel, if there is still a problem, please do re-open it (or file a new one).

Ben Pfaff <blp>
Group administrator
Fri 06 Feb 2009 05:53:24 AM UTC, comment #16: 

You need also to remove the references to SIGALRM in main.c and psppire.c

Other than that, I'm happy with these patches.

John Darrington <jmd>
Group administrator
Fri 06 Feb 2009 05:32:13 AM UTC, comment #15: 


>...less than 15 lines isn't "legally significant"...


OK.  Here is take 2, then.

I'm also attaching an actual fix for this bug.

(file #17410, file #17411)

Ben Pfaff <blp>
Group administrator
Thu 05 Feb 2009 10:45:12 PM UTC, comment #14: 

The guidelines at http://www.gnu.org/prep/maintain/maintain.html#Legally-Significant say that less than 15 lines isn't "legally significant", so I don't think we need to worry about these 4.

In the interests of proper attribution, perhaps we should retain a comment indicating the original author.

John Darrington <jmd>
Group administrator
Thu 05 Feb 2009 03:27:42 PM UTC, comment #13: 


>As there's no longer any 3rd party code involved, there's no
>reason for it to stay in lib/misc. It can go in src/math.


I was assuming that the first few lines of LevelOfSignificanceWXMPSR() came from Rob van Son, e.g. this code:


  /* Determine Wmax, i.e., work with the largest Rank Sum */
  MaximalW = N*(N+1)/2;
  if(Winput < MaximalW/2)Winput = MaximalW - Winput;
  W = Winput;    /* Convert to long int */
  if(W != Winput)++W;  /* Increase to next full integer */


Given that, I was staying on the safe side by leaving it in lib/misc.

Am I wrong?

Ben Pfaff <blp>
Group administrator
Thu 05 Feb 2009 07:30:52 AM UTC, comment #12: 

As there's no longer any 3rd party code involved, there's no reason for it to stay in lib/misc.  It can go in src/math.

John Darrington <jmd>
Group administrator
Thu 05 Feb 2009 06:56:23 AM UTC, comment #11: 

Here is a proposed patch to fix the performance problem.  I have not yet written up a patch to remove or change the timer code, but pushing this commit would enable that.

(file #17404)

Ben Pfaff <blp>
Group administrator
Thu 05 Feb 2009 06:27:31 AM UTC, comment #10: 


>a) if we end up using routines from a third party library, then
>this won't be available to us


longjmp() from a signal handler out of arbitrary third-party code seems to be courting disaster.  At the very least it invites memory leaks.

>b) Executing the test condition would have slowed the system even
>more.


If you do it right then you don't put the test in your innermost loop and the test is only a pair of test-and-jump instructions that the CPU branch predictor always gets right.

The model checking code contains something like this.

Ben Pfaff <blp>
Group administrator
Thu 05 Feb 2009 05:34:07 AM UTC, comment #9: 


> Do we really need to use setjmp()/longjmp() to break out of the algorithm though? If we avoid that, by inserting occasional code sequences of the form "if (timed_out()) return ERROR_CODE;", then it would be easy to make it portable.


I did consider that when I first wrote this command.  I decided not to because a) if we end up using routines from a third party library, then this won't be available to us; b) Executing the test condition would have slowed the system even more.

However, it might be the best option now that we have (possibly) a fast algorithm of our own.

John Darrington <jmd>
Group administrator
Wed 04 Feb 2009 03:27:00 PM UTC, comment #8: 


>Even if that's true for this case, there will be many algorithms
>in other (yet to be implemented) subcommands in NPAR TESTS which
>are intractable. We'll need the timer for those (and to be keep
>spss refugees happy).


Definitely.

Do we really need to use setjmp()/longjmp() to break out of the algorithm though?  If we avoid that, by inserting occasional code sequences of the form "if (timed_out()) return ERROR_CODE;", then it would be easy to make it portable.

Ben Pfaff <blp>
Group administrator
Wed 04 Feb 2009 03:24:18 PM UTC, comment #7: 

Actually, here's a solution that is even better.  A test program using it calculates (correctly) and prints all the values for 1 <= N <= 30, 1 <= W <= N*(N+1)/2, in under .2 seconds total.

The runtime of this algorithm is O(N*W) and it uses O(W) space.


int ranksum6(int N, int W)
{
  int *array;
  int max;
  int total;

  assert (N >= 0);
  if (N == 0)
    return 0;
  else if (W <= 0)
    return 1 << N;
  else if (W > N * (N + 1) / 2)
    return 0;
  else if (N == 1)
    return 1;

  array = calloc (sizeof *array, W + 1);
  array[W] = 1;

  max = W;
  total = 0;
  for (; N > 1; N--)
    {
      int max_sum = N * (N + 1) / 2;
      int i;

      if (max_sum < max)
        max = max_sum;

      for (i = 1; i <= max; i++)
        if (array[i] != 0)
          {
            int new_W = i - N;
            if (new_W <= 0)
              total += array[i] * (1 << (N - 1));
            else
              array[new_W] += array[i];
          }
    }
  total += array[1];
  free (array);
  return total;
}


Ben Pfaff <blp>
Group administrator
Tue 03 Feb 2009 10:22:17 PM UTC, comment #6: 


>  there's a lot of redundant calculation in ranksum2(). But that's still a big improvement, and good enough to get rid of the timer entirely I think.



Even if that's true for this case, there will be many algorithms in other (yet to be implemented) subcommands in NPAR TESTS which are intractable.  We'll need the timer for those (and to be keep spss refugees happy).

John Darrington <jmd>
Group administrator
Tue 03 Feb 2009 10:18:15 PM UTC, comment #5: 

This is certainly very interesting.   I spent quite a lot of time searching the literature for a more efficient algorithm, and concluded that there wasn't one.

0.  Can you demonstrate that the two algorithms are identical?

1.  What is the complexity of your new algorithm?  It's not clear to me.

2.  I'd like to see a statistician comment on this.  The original algorithm's author asked for improvements to be sent to him. His home page is http://www.fon.hum.uva.nl/rob/  ; perhaps it would be appropriate to ask for his comments.

3.  If this works, you should publish a paper on it.

John Darrington <jmd>
Group administrator
Tue 03 Feb 2009 07:59:45 AM UTC, comment #4: 


>Unless anyone has any better ideas, I suggest that we just compile
>out support for the TIMER subcommand on platforms which don't have
>a SIGALRM.


The real problem, I think, is that LevelOfSignificanceWXMPSR is astoundingly inefficient.  The core of it is essentially this:


int ranksum1(int N, unsigned long W)
{
  unsigned long NumberOfPossibilities = 1 << N;
  int CountLarger = 0;
  unsigned long i;

  /* Generate all distributions of sign over ranks as bit-patterns (i). */
  for(i=0; i < NumberOfPossibilities; ++i)
  {
    unsigned long RankSum = 0;
    int j;

    /*
       Shift "sign" bits out of i to determine the Sum of Ranks (j).
    */
    for(j=0; j < N; ++j)
    {
      if((i >> j) & 1)RankSum += j + 1;
    };
    /*
    * Count the number of "samples" that have a Sum of Ranks larger than
    * or equal to the one found (i.e., >= W).
    */
    if(RankSum >= W)++CountLarger;
  }
  return CountLarger;
}


But that can be reimplemented in a much more efficient fashion as:


int ranksum2(int N, int W)
{
  if (N > 0 && W > 0)
    {
      unsigned long int total = 0;
      while (N >= W && N > 0)
        total += 1 << --N;
      return total + ranksum2(N - 1, W - N) + ranksum2(N - 1, W);
    }
  else
    return 0;
}


For N=30, W=120 ranksum1() takes 6 minutes 9 seconds on my system, ranksum2() takes 0.4 seconds.  And they both return the same result.

(Am I missing anything here?  This seems too easy.)

With dynamic programming I think I can do even better; there's a lot of redundant calculation in ranksum2().  But that's still a big improvement, and good enough to get rid of the timer entirely I think.

I'm attaching a test program that can be used to verify that both functions return the same result and to benchmark.

(file #17392)

Ben Pfaff <blp>
Group administrator
Tue 03 Feb 2009 01:03:01 AM UTC, comment #3: 

A bigger problem is the lack of SIGALRM on mingw.  Sadly there doesn't seem to be anything in gnulib which can help us either.

Unless anyone has any better ideas, I suggest that we just compile out support for the TIMER subcommand on platforms which don't have a SIGALRM.

John Darrington <jmd>
Group administrator
Mon 02 Feb 2009 06:50:45 PM UTC, comment #2: 

We should probably fall back to using setjmp instead of sigsetjmp on mingw.  It ought to work fine.

Ben Pfaff <blp>
Group administrator
Mon 02 Feb 2009 06:44:31 PM UTC, comment #1: 

Im looking closer to this bug, and seems like mingw doesn't have the implementation of sigjmp_buf. I will search if it will be done in the future, but meanwhile we should find a workaround, because this problem doesn't let new versions of pspp be compiled to windows.

Michel Boaventura <michelboaventura>
Group Member
Mon 02 Feb 2009 06:35:48 PM UTC, original submission:  

OS: Gentoo using mingw to cross-compile
PSPP version: git master branch

Get this compile error:

src/language/stats/wilcoxon.c:373: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or ‘__attribute__’ before ‘env’
src/language/stats/wilcoxon.c: In function ‘give_up_callback’:
src/language/stats/wilcoxon.c:378: warning: implicit declaration of function ‘siglongjmp’
src/language/stats/wilcoxon.c:378: error: ‘env’ undeclared (first use in this function)
src/language/stats/wilcoxon.c:378: error: (Each undeclared identifier is reported only once
src/language/stats/wilcoxon.c:378: error: for each function it appears in.)
src/language/stats/wilcoxon.c: In function ‘timed_wilcoxon_significance’:
src/language/stats/wilcoxon.c:401: warning: implicit declaration of function ‘sigsetjmp’
src/language/stats/wilcoxon.c:401: error: ‘env’ undeclared (first use in this function)
src/language/stats/wilcoxon.c:403: error: ‘SIGALRM’ undeclared (first use in this function)
src/language/stats/wilcoxon.c:404: warning: implicit declaration of function ‘alarm’
make[2]: ** [src/language/stats/wilcoxon.lo] Erro 1
make[2]: Saindo do diretório `/home/michel/pspp'
make[1]: ** [all-recursive] Erro 1
make[1]: Saindo do diretório `/home/michel/pspp'
make: ** [all] Erro 2

Michel Boaventura <michelboaventura>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #17410:  wilcoxon-speedup.patch added by blp (12KiB - text/x-diff)
file #17411:  bug-25466.patch added by blp (2KiB - text/x-diff)
file #17404:  foo.patch added by blp (7KiB - text/x-diff)
file #17392:  ranksum.c added by blp (2KiB - text/x-csrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jmd (Posted a comment)
  • -email is unavailable- added by blp (Posted a comment)
  • -email is unavailable- added by michelboaventura (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2009-02-06 blp StatusNone Fixed
        Open/ClosedOpen Closed
    2009-02-06 blp Attached File- Added wilcoxon-speedup.patch, #17410
        Attached File- Added bug-25466.patch, #17411
    2009-02-05 blp Attached File- Added foo.patch, #17404
    2009-02-03 blp Attached File- Added ranksum.c, #17392

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code