Speed up Wilcoxon matched-pairs signed-ranks significance routine.

[pspp-builds.git] / src / language / stats / wilcoxon.c

diff --git a/src/language/stats/wilcoxon.c b/src/language/stats/wilcoxon.c
index 7552d7ece1f84e37acd189de6aaeebf4ce11ba2a..a485f3506fd8b226e956ca88b3b126dc5297bec3 100644 (file)
--- a/src/language/stats/wilcoxon.c
+++ b/src/language/stats/wilcoxon.c
@@ -1,5 +1,5 @@
 /* Pspp - a program for statistical analysis.
-   Copyright (C) 2008 Free Software Foundation, Inc.
+   Copyright (C) 2008, 2009 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -21,14 +21,14 @@
 #include <data/variable.h>
 #include <data/casereader.h>
 #include <data/casewriter.h>
-#include <data/case-ordering.h>
+#include <data/subcase.h>
 #include <math/sort.h>
 #include <libpspp/message.h>
 #include <xalloc.h>
 #include <output/table.h>
 #include <data/procedure.h>
 #include <data/dictionary.h>
-#include <misc/wx-mp-sr.h>
+#include <math/wilcoxon-sig.h>
 #include <gsl/gsl_cdf.h>
 #include <unistd.h>
 #include <signal.h>
@@ -88,8 +88,8 @@ wilcoxon_execute (const struct dataset *ds,
     {
       struct casereader *r = casereader_clone (input);
       struct casewriter *writer;
-      struct ccase c;
-      struct case_ordering *ordering = case_ordering_create ();
+      struct ccase *c;
+      struct subcase ordering;
       variable_pair *vp = &t2s->pairs[i];
 
       const int reader_width = weight ? 3 : 2;
@@ -97,49 +97,45 @@ wilcoxon_execute (const struct dataset *ds,
       ws[i].sign = var_create_internal (0);
       ws[i].absdiff = var_create_internal (1);
 
-      case_ordering_add_var (ordering, ws[i].absdiff, SRT_ASCEND);
-
-
       r = casereader_create_filter_missing (r, *vp, 2,
                                            exclude,
                                            NULL, NULL);
 
-      writer = sort_create_writer (ordering, reader_width);
-      while (casereader_read (r, &c))
-       {
-         struct ccase output;
-         double d = append_difference (&c, 0, vp);
+      subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
+      writer = sort_create_writer (&ordering, reader_width);
+      subcase_destroy (&ordering);
 
-         case_create (&output, reader_width);
+      for (; (c = casereader_read (r)) != NULL; case_unref (c))
+       {
+         struct ccase *output = case_create (reader_width);
+         double d = append_difference (c, 0, vp);
 
          if (d > 0)
            {
-             case_data_rw (&output, ws[i].sign)->f = 1.0;
+             case_data_rw (output, ws[i].sign)->f = 1.0;
 
            }
          else if (d < 0)
            {
-             case_data_rw (&output, ws[i].sign)->f = -1.0;
+             case_data_rw (output, ws[i].sign)->f = -1.0;
            }
          else
            {
              double w = 1.0;
              if (weight)
-               w = case_data (&c, weight)->f;
+               w = case_data (c, weight)->f;
 
              /* Central point values should be dropped */
              ws[i].n_zeros += w;
-             case_destroy (&c);
-             continue;
+              continue;
            }
 
-         case_data_rw (&output, ws[i].absdiff)->f = fabs (d);
+         case_data_rw (output, ws[i].absdiff)->f = fabs (d);
 
          if (weight)
-          case_data_rw (&output, weightx)->f = case_data (&c, weight)->f;
+          case_data_rw (output, weightx)->f = case_data (c, weight)->f;
 
-         casewriter_write (writer, &output);
-         case_destroy (&c);
+         casewriter_write (writer, output);
        }
       casereader_destroy (r);
       ws[i].reader = casewriter_make_reader (writer);
@@ -148,7 +144,7 @@ wilcoxon_execute (const struct dataset *ds,
   for (i = 0 ; i < t2s->n_pairs; ++i )
     {
       struct casereader *rr ;
-      struct ccase c;
+      struct ccase *c;
       enum rank_error err = 0;
 
       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
@@ -156,13 +152,13 @@ wilcoxon_execute (const struct dataset *ds,
                                          distinct_callback, &ws[i]
                                          );
 
-      while (casereader_read (rr, &c))
+      for (; (c = casereader_read (rr)) != NULL; case_unref (c))
        {
-         double sign = case_data (&c, ws[i].sign)->f;
-         double rank = case_data_idx (&c, weight ? 3 : 2)->f;
+         double sign = case_data (c, ws[i].sign)->f;
+         double rank = case_data_idx (c, weight ? 3 : 2)->f;
          double w = 1.0;
          if (weight)
-           w = case_data (&c, weightx)->f;
+           w = case_data (c, weightx)->f;
 
          if ( sign > 0 )
            {
@@ -176,8 +172,6 @@ wilcoxon_execute (const struct dataset *ds,
            }
          else
            NOT_REACHED ();
-
-         case_destroy (&c);
        }
 
       casereader_destroy (rr);