* Changed dict_get_case_weight() to accept an additional int * flag to complain about...
[pspp] / src / crosstabs.q
index d880ceb03ae3f5e673044a74895181664248081b..4a58d10c96aa9979395f4eba9e6ccd8ff97ce1e7 100644 (file)
 */
 
 #include <config.h>
-#include <assert.h>
+#include "error.h"
 #include <ctype.h>
 #include <stdlib.h>
 #include <stdio.h>
+#include <gsl/gsl_cdf.h>
 #include "algorithm.h"
 #include "alloc.h"
 #include "hash.h"
@@ -43,7 +44,6 @@
 #include "error.h"
 #include "magic.h"
 #include "misc.h"
-#include "stats.h"
 #include "output.h"
 #include "tab.h"
 #include "value-labels.h"
@@ -64,7 +64,7 @@
             tabl:!tables/notables,
             box:!box/nobox,
             pivot:!pivot/nopivot;
-     +cells[cl_]=count,none,row,column,total,expected,residual,sresidual,
+     +cells[cl_]=count,none,expected,row,column,total,residual,sresidual,
                 asresidual,all;
      +statistics[st_]=chisq,phi,cc,lambda,uc,none,btau,ctau,risk,gamma,d,
                      kappa,eta,corr,all.
@@ -178,6 +178,8 @@ cmd_crosstabs (void)
 static int
 internal_cmd_crosstabs (void)
 {
+  int i;
+
   variables = NULL;
   variables_cnt = 0;
   xtab = NULL;
@@ -185,7 +187,6 @@ internal_cmd_crosstabs (void)
   pl_tc = pool_create ();
   pl_col = pool_create ();
 
-  lex_match_id ("CROSSTABS");
   if (!parse_crosstabs (&cmd))
     return CMD_FAILURE;
 
@@ -200,11 +201,9 @@ internal_cmd_crosstabs (void)
   if (!cmd.sbc_cells)
     {
       cmd.a_cells[CRS_CL_COUNT] = 1;
-      num_cells = 1;
     }
   else 
     {
-      int i;
       int count = 0;
 
       for (i = 0; i < CRS_CL_count; i++)
@@ -224,10 +223,10 @@ internal_cmd_crosstabs (void)
          cmd.a_cells[CRS_CL_ALL] = 0;
        }
       cmd.a_cells[CRS_CL_NONE] = 0;
-      for (num_cells = i = 0; i < CRS_CL_count; i++)
-       if (cmd.a_cells[i])
-          cmd.a_cells[num_cells++] = i;
     }
+  for (num_cells = i = 0; i < CRS_CL_count; i++)
+    if (cmd.a_cells[i])
+      cells[num_cells++] = i;
 
   /* STATISTICS. */
   if (cmd.sbc_statistics)
@@ -274,8 +273,9 @@ internal_cmd_crosstabs (void)
   else
     write = CRS_WR_NONE;
 
-  procedure (precalc, mode == GENERAL ? calc_general : calc_integer, postcalc,
-             NULL);
+  procedure_with_splits (precalc,
+                         mode == GENERAL ? calc_general : calc_integer,
+                         postcalc, NULL);
 
   return CMD_SUCCESS;
 }
@@ -584,8 +584,10 @@ precalc (void *aux UNUSED)
 static int
 calc_general (struct ccase *c, void *aux UNUSED)
 {
+  int bad_warn = 1;
+
   /* Case weight. */
-  double weight = dict_get_case_weight (default_dict, c);
+  double weight = dict_get_case_weight (default_dict, c, &bad_warn);
 
   /* Flattened current table index. */
   int t;
@@ -656,8 +658,10 @@ calc_general (struct ccase *c, void *aux UNUSED)
 static int
 calc_integer (struct ccase *c, void *aux UNUSED)
 {
+  int bad_warn = 1;
+
   /* Case weight. */
-  double weight = dict_get_case_weight (default_dict, c);
+  double weight = dict_get_case_weight (default_dict, c, &bad_warn);
   
   /* Flattened current table index. */
   int t;
@@ -1766,22 +1770,31 @@ display_dimensions (struct tab_table *table, int first_difference, struct table_
                         x->vars[first_difference]);
 }
 
-/* Put value V into cell (C,R) of TABLE, suffixed with letter M. */
+/* Put VALUE into cell (C,R) of TABLE, suffixed with character
+   SUFFIX if nonzero.  If MARK_MISSING is nonzero the entry is
+   additionally suffixed with a letter `M'. */
 static void
-float_M_suffix (struct tab_table *table, int c, int r, double v)
+format_cell_entry (struct tab_table *table, int c, int r, double value,
+                   char suffix, int mark_missing)
 {
-  static const struct fmt_spec f = {FMT_F, 8, 0};
+  const struct fmt_spec f = {FMT_F, 10, 1};
+  union value v;
   struct len_string s;
-
-  s.length = 9;
-  s.string = tab_alloc (table, 9);
-  s.string[8] = 'M';
-  format_short (s.string, &f, (union value *) &v);
+  
+  s.length = 10;
+  s.string = tab_alloc (table, 16);
+  v.f = value;
+  data_out (s.string, &f, &v);
   while (*s.string == ' ')
     {
       s.length--;
       s.string++;
     }
+  if (suffix != 0)
+    s.string[s.length++] = suffix;
+  if (mark_missing)
+    s.string[s.length++] = 'M';
+
   tab_raw (table, c, r, TAB_RIGHT, &s);
 }
 
@@ -1811,10 +1824,16 @@ display_crosstabulation (void)
          tab_hline (table, TAL_1, -1, n_cols, 0);
        for (c = 0; c < n_cols; c++)
          {
-           double expected_value = row_tot[r] * col_tot[c] / W;
+            int mark_missing = 0;
+            double expected_value = row_tot[r] * col_tot[c] / W;
+            if (cmd.miss == CRS_REPORT
+                && (is_num_user_missing (cols[c].f, x->vars[COL_VAR])
+                    || is_num_user_missing (rows[r].f, x->vars[ROW_VAR])))
+              mark_missing = 1;
            for (i = 0; i < num_cells; i++)
              {
                double v;
+                int suffix = 0;
 
                switch (cells[i])
                  {
@@ -1823,12 +1842,15 @@ display_crosstabulation (void)
                    break;
                  case CRS_CL_ROW:
                    v = *mp / row_tot[r] * 100.;
+                    suffix = '%';
                    break;
                  case CRS_CL_COLUMN:
                    v = *mp / col_tot[c] * 100.;
+                    suffix = '%';
                    break;
                  case CRS_CL_TOTAL:
                    v = *mp / W * 100.;
+                    suffix = '%';
                    break;
                  case CRS_CL_EXPECTED:
                    v = expected_value;
@@ -1847,14 +1869,10 @@ display_crosstabulation (void)
                    break;
                  default:
                    assert (0);
+                    abort ();
                  }
 
-               if (cmd.miss == CRS_REPORT
-                   && (is_num_user_missing (cols[c].f, x->vars[COL_VAR])
-                       || is_num_user_missing (rows[r].f, x->vars[ROW_VAR])))
-                 float_M_suffix (table, c, i, v);
-               else if (v != 0.)
-                 tab_float (table, c, i, TAB_RIGHT, v, 8, 0);
+                format_cell_entry (table, c, i, v, suffix, mark_missing);
              }
 
            mp++;
@@ -1869,48 +1887,56 @@ display_crosstabulation (void)
     int r, i;
 
     tab_offset (table, -1, tab_row (table) - num_cells * n_rows);
-    for (r = 0; r < n_rows; r++)
-      for (i = 0; i < num_cells; i++)
-       {
-         double v;
-
-         switch (cells[i])
-           {
-           case CRS_CL_COUNT:
-             v = row_tot[r];
-             break;
-           case CRS_CL_ROW:
-             v = 100.;
-             break;
-           case CRS_CL_COLUMN:
-             v = row_tot[r] / W * 100.;
-             break;
-           case CRS_CL_TOTAL:
-             v = row_tot[r] / W * 100.;
-             break;
-           case CRS_CL_EXPECTED:
-           case CRS_CL_RESIDUAL:
-           case CRS_CL_SRESIDUAL:
-           case CRS_CL_ASRESIDUAL:
-             v = 0.;
-             break;
-           default:
-             assert (0);
-           }
-
-         if (cmd.miss == CRS_REPORT
-             && is_num_user_missing (rows[r].f, x->vars[ROW_VAR]))
-           float_M_suffix (table, n_cols, 0, v);
-         else if (v != 0.)
-           tab_float (table, n_cols, 0, TAB_RIGHT, v, 8, 0);
-
-         tab_next_row (table);
-       }
+    for (r = 0; r < n_rows; r++) 
+      {
+        char suffix = 0;
+        int mark_missing = 0;
+
+        if (cmd.miss == CRS_REPORT
+            && is_num_user_missing (rows[r].f, x->vars[ROW_VAR]))
+          mark_missing = 1;
+
+        for (i = 0; i < num_cells; i++)
+          {
+            double v;
+
+            switch (cells[i])
+              {
+              case CRS_CL_COUNT:
+                v = row_tot[r];
+                break;
+              case CRS_CL_ROW:
+                v = 100.;
+                suffix = '%';
+                break;
+              case CRS_CL_COLUMN:
+                v = row_tot[r] / W * 100.;
+                suffix = '%';
+                break;
+              case CRS_CL_TOTAL:
+                v = row_tot[r] / W * 100.;
+                suffix = '%';
+                break;
+              case CRS_CL_EXPECTED:
+              case CRS_CL_RESIDUAL:
+              case CRS_CL_SRESIDUAL:
+              case CRS_CL_ASRESIDUAL:
+                v = 0.;
+                break;
+              default:
+                assert (0);
+                abort ();
+              }
+
+            format_cell_entry (table, n_cols, 0, v, suffix, mark_missing);
+            tab_next_row (table);
+          } 
+      }
   }
 
   /* Column totals, grand total. */
   {
-    int c, j;
+    int c;
     int last_row = 0;
 
     if (num_cells > 1)
@@ -1918,9 +1944,15 @@ display_crosstabulation (void)
     for (c = 0; c <= n_cols; c++)
       {
        double ct = c < n_cols ? col_tot[c] : W;
-       int i;
+        int mark_missing = 0;
+        char suffix = 0;
+        int i;
            
-       for (i = j = 0; i < num_cells; i++)
+        if (cmd.miss == CRS_REPORT && c < n_cols 
+            && is_num_user_missing (cols[c].f, x->vars[COL_VAR]))
+          mark_missing = 1;
+
+        for (i = 0; i < num_cells; i++)
          {
            double v;
 
@@ -1928,15 +1960,19 @@ display_crosstabulation (void)
              {
              case CRS_CL_COUNT:
                v = ct;
+                suffix = '%';
                break;
              case CRS_CL_ROW:
                v = ct / W * 100.;
+                suffix = '%';
                break;
              case CRS_CL_COLUMN:
                v = 100.;
+                suffix = '%';
                break;
              case CRS_CL_TOTAL:
                v = ct / W * 100.;
+                suffix = '%';
                break;
              case CRS_CL_EXPECTED:
              case CRS_CL_RESIDUAL:
@@ -1945,17 +1981,12 @@ display_crosstabulation (void)
                continue;
              default:
                assert (0);
+                abort ();
              }
 
-           if (cmd.miss == CRS_REPORT && c < n_cols 
-               && is_num_user_missing (cols[c].f, x->vars[COL_VAR]))
-             float_M_suffix (table, c, j, v);
-           else if (v != 0.)
-             tab_float (table, c, j, TAB_RIGHT, v, 8, 0);
-
-           j++;
+            format_cell_entry (table, c, i, v, suffix, mark_missing);
          }
-        last_row = j;
+        last_row = i;
       }
 
     tab_offset (table, -1, tab_row (table) + last_row);
@@ -2003,7 +2034,7 @@ display_chisq (void)
          tab_float (chisq, 1, 0, TAB_RIGHT, chisq_v[i], 8, 3);
          tab_float (chisq, 2, 0, TAB_RIGHT, df[i], 8, 0);
          tab_float (chisq, 3, 0, TAB_RIGHT,
-                    chisq_sig (chisq_v[i], df[i]), 8, 3);
+                    gsl_cdf_chisq_Q (chisq_v[i], df[i]), 8, 3);
        }
       else
        {
@@ -2668,10 +2699,10 @@ calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
 
                if (cmd.a_statistics[CRS_ST_D])
                  {
-                   d_yx_cum += fij * sqr (Dr * (Cij - Dij)
-                                          - (P - Q) * (W - row_tot[i]));
-                   d_xy_cum += fij * sqr (Dc * (Dij - Cij)
-                                          - (Q - P) * (W - col_tot[j]));
+                   d_yx_cum += fij * pow2 (Dr * (Cij - Dij)
+                                            - (P - Q) * (W - row_tot[i]));
+                   d_xy_cum += fij * pow2 (Dc * (Dij - Cij)
+                                            - (Q - P) * (W - col_tot[j]));
                  }
                
                if (++j == n_cols)
@@ -2691,8 +2722,8 @@ calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
       }
 
       btau_var = ((btau_cum
-                  - (W * sqr (W * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
-                 / sqr (Dr * Dc));
+                  - (W * pow2 (W * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
+                 / pow2 (Dr * Dc));
       if (cmd.a_statistics[CRS_ST_BTAU])
        {
          ase[3] = sqrt (btau_var);
@@ -2717,17 +2748,17 @@ calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
          somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
          somers_d_t[0] = (somers_d_v[0]
                           / (4 / (Dc + Dr)
-                             * sqrt (ctau_cum - sqr (P - Q) / W)));
+                             * sqrt (ctau_cum - pow2 (P - Q) / W)));
          somers_d_v[1] = (P - Q) / Dc;
-         somers_d_ase[1] = 2. / sqr (Dc) * sqrt (d_xy_cum);
+         somers_d_ase[1] = 2. / pow2 (Dc) * sqrt (d_xy_cum);
          somers_d_t[1] = (somers_d_v[1]
                           / (2. / Dc
-                             * sqrt (ctau_cum - sqr (P - Q) / W)));
+                             * sqrt (ctau_cum - pow2 (P - Q) / W)));
          somers_d_v[2] = (P - Q) / Dr;
-         somers_d_ase[2] = 2. / sqr (Dr) * sqrt (d_yx_cum);
+         somers_d_ase[2] = 2. / pow2 (Dr) * sqrt (d_yx_cum);
          somers_d_t[2] = (somers_d_v[2]
                           / (2. / Dr
-                             * sqrt (ctau_cum - sqr (P - Q) / W)));
+                             * sqrt (ctau_cum - pow2 (P - Q) / W)));
        }
 
       free (cum);
@@ -2820,12 +2851,12 @@ calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
                     / (W * (W * W - sum_rici) * (W * W - sum_rici)));
 #if 0
       t[8] = v[8] / sqrt (W * (((sum_fii * (W - sum_fii))
-                               / sqr (W * W - sum_rici))
+                               / pow2 (W * W - sum_rici))
                               + ((2. * (W - sum_fii)
                                   * (2. * sum_fii * sum_rici
                                      - W * sum_fiiri_ci))
                                  / cube (W * W - sum_rici))
-                              + (sqr (W - sum_fii)
+                              + (pow2 (W - sum_fii)
                                  * (W * sum_fijri_ci2 - 4.
                                     * sum_rici * sum_rici)
                                  / pow4 (W * W - sum_rici))));
@@ -2988,7 +3019,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
            {
              const int deltaj = j == cm_index;
              accum += (mat[j + i * n_cols]
-                       * sqr ((j == fim_index[i])
+                       * pow2 ((j == fim_index[i])
                               - deltaj
                               + v[0] * deltaj));
            }
@@ -3004,7 +3035,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
          if (cm_index != fim_index[i])
            accum += (mat[i * n_cols + fim_index[i]]
                      + mat[i * n_cols + cm_index]);
-       t[2] = v[2] / (sqrt (accum - sqr (sum_fim - cm) / W) / (W - cm));
+       t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / W) / (W - cm));
       }
 
       /* ASE1 for X given Y. */
@@ -3016,7 +3047,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
            {
              const int deltaj = i == rm_index;
              accum += (mat[j + i * n_cols]
-                       * sqr ((i == fmj_index[j])
+                       * pow2 ((i == fmj_index[j])
                               - deltaj
                               + v[0] * deltaj));
            }
@@ -3032,7 +3063,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
          if (rm_index != fmj_index[j])
            accum += (mat[j + n_cols * fmj_index[j]]
                      + mat[j + n_cols * rm_index]);
-       t[1] = v[1] / (sqrt (accum - sqr (sum_fmj - rm) / W) / (W - rm));
+       t[1] = v[1] / (sqrt (accum - pow2 (sum_fmj - rm) / W) / (W - rm));
       }
 
       /* Symmetric ASE0 and ASE1. */
@@ -3045,12 +3076,12 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
            {
              int temp0 = (fmj_index[j] == i) + (fim_index[i] == j);
              int temp1 = (i == rm_index) + (j == cm_index);
-             accum0 += mat[j + i * n_cols] * sqr (temp0 - temp1);
+             accum0 += mat[j + i * n_cols] * pow2 (temp0 - temp1);
              accum1 += (mat[j + i * n_cols]
-                        * sqr (temp0 + (v[0] - 1.) * temp1));
+                        * pow2 (temp0 + (v[0] - 1.) * temp1));
            }
        ase[0] = sqrt (accum1 - 4. * W * v[0] * v[0]) / (2. * W - rm - cm);
-       t[0] = v[0] / (sqrt (accum0 - sqr ((sum_fim + sum_fmj - cm - rm) / W))
+       t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / W))
                       / (2. * W - rm - cm));
       }
 
@@ -3066,7 +3097,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
        for (sum_fij2_ri = sum_fij2_ci = 0., i = 0; i < n_rows; i++)
          for (j = 0; j < n_cols; j++)
            {
-             double temp = sqr (mat[j + i * n_cols]);
+             double temp = pow2 (mat[j + i * n_cols]);
              sum_fij2_ri += temp / row_tot[i];
              sum_fij2_ci += temp / col_tot[j];
            }
@@ -3104,7 +3135,7 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
            if (entry <= 0.)
              continue;
            
-           P += entry * sqr (log (col_tot[j] * row_tot[i] / (W * entry)));
+           P += entry * pow2 (log (col_tot[j] * row_tot[i] / (W * entry)));
            UXY -= entry / W * log (entry / W);
          }
 
@@ -3116,27 +3147,27 @@ calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
            if (entry <= 0.)
              continue;
            
-           ase1_yx += entry * sqr (UY * log (entry / row_tot[i])
+           ase1_yx += entry * pow2 (UY * log (entry / row_tot[i])
                                    + (UX - UXY) * log (col_tot[j] / W));
-           ase1_xy += entry * sqr (UX * log (entry / col_tot[j])
+           ase1_xy += entry * pow2 (UX * log (entry / col_tot[j])
                                    + (UY - UXY) * log (row_tot[i] / W));
-           ase1_sym += entry * sqr ((UXY
+           ase1_sym += entry * pow2 ((UXY
                                      * log (row_tot[i] * col_tot[j] / (W * W)))
                                     - (UX + UY) * log (entry / W));
          }
       
       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
-      ase[5] = (2. / (W * sqr (UX + UY))) * sqrt (ase1_sym);
+      ase[5] = (2. / (W * pow2 (UX + UY))) * sqrt (ase1_sym);
       t[5] = v[5] / ((2. / (W * (UX + UY)))
-                    * sqrt (P - sqr (UX + UY - UXY) / W));
+                    * sqrt (P - pow2 (UX + UY - UXY) / W));
                    
       v[6] = (UX + UY - UXY) / UX;
       ase[6] = sqrt (ase1_xy) / (W * UX * UX);
-      t[6] = v[6] / (sqrt (P - W * sqr (UX + UY - UXY)) / (W * UX));
+      t[6] = v[6] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UX));
       
       v[7] = (UX + UY - UXY) / UY;
       ase[7] = sqrt (ase1_yx) / (W * UY * UY);
-      t[7] = v[7] / (sqrt (P - W * sqr (UX + UY - UXY)) / (W * UY));
+      t[7] = v[7] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UY));
     }
 
   /* Somers' D. */