pintos-os.org Git - pspp/blob - src/math/covariance.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009, 2010, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include "math/covariance.h"
  20
  21 #include <gsl/gsl_matrix.h>
  22
  23 #include "data/case.h"
  24 #include "data/variable.h"
  25 #include "libpspp/assertion.h"
  26 #include "libpspp/misc.h"
  27 #include "math/categoricals.h"
  28 #include "math/interaction.h"
  29 #include "math/moments.h"
  30
  31 #include "gl/xalloc.h"
  32
  33 #define n_MOMENTS (MOMENT_VARIANCE + 1)
  34
  35
  36 /* Create a new matrix of NEW_SIZE x NEW_SIZE and copy the elements of
  37    matrix IN into it.  IN must be a square matrix, and in normal usage
  38    it will be smaller than NEW_SIZE.
  39    IN is destroyed by this function.  The return value must be destroyed
  40    when no longer required.
  41 */
  42 static gsl_matrix *
  43 resize_matrix (gsl_matrix *in, size_t new_size)
  44 {
  45   size_t i, j;
  46
  47   gsl_matrix *out = NULL;
  48
  49   assert (in->size1 == in->size2);
  50
  51   if (new_size <= in->size1)
  52     return in;
  53
  54   out = gsl_matrix_calloc (new_size, new_size);
  55
  56   for (i = 0; i < in->size1; ++i)
  57     {
  58       for (j = 0; j < in->size2; ++j)
  59         {
  60           double x = gsl_matrix_get (in, i, j);
  61
  62           gsl_matrix_set (out, i, j, x);
  63         }
  64     }
  65
  66   gsl_matrix_free (in);
  67
  68   return out;
  69 }
  70
  71 struct covariance
  72 {
  73   /* True if the covariances are centerered. (ie Real covariances) */
  74   bool centered;
  75
  76   /* The variables for which the covariance matrix is to be calculated. */
  77   size_t n_vars;
  78   const struct variable *const *vars;
  79
  80   /* Categorical variables. */
  81   struct categoricals *categoricals;
  82
  83   /* Array containing number of categories per categorical variable. */
  84   size_t *n_categories;
  85
  86   /* Dimension of the covariance matrix. */
  87   size_t dim;
  88
  89   /* The weight variable (or NULL if none) */
  90   const struct variable *wv;
  91
  92   /* A set of matrices containing the 0th, 1st and 2nd moments */
  93   gsl_matrix **moments;
  94
  95   /* The class of missing values to exclude */
  96   enum mv_class exclude;
  97
  98   /* An array of doubles representing the covariance matrix.
  99      Only the top triangle is included, and no diagonals */
 100   double *cm;
 101   int n_cm;
 102
 103   /* 1 for single pass algorithm;
 104      2 for double pass algorithm
 105   */
 106   short passes;
 107
 108   /*
 109     0 : No pass has  been made
 110     1 : First pass has been started
 111     2 : Second pass has been
 112
 113     IE: How many passes have been (partially) made. */
 114   short state;
 115
 116   /* Flags indicating that the first case has been seen */
 117   bool pass_one_first_case_seen;
 118   bool pass_two_first_case_seen;
 119
 120   gsl_matrix *unnormalised;
 121 };
 122
 123
 124
 125 /* Return a matrix containing the M th moments.
 126    The matrix is of size  NxN where N is the number of variables.
 127    Each row represents the moments of a variable.
 128    In the absence of missing values, the columns of this matrix will
 129    be identical.  If missing values are involved, then element (i,j)
 130    is the moment of the i th variable, when paired with the j th variable.
 131  */
 132 gsl_matrix *
 133 covariance_moments (const struct covariance *cov, int m)
 134 {
 135   return cov->moments[m];
 136 }
 137
 138
 139
 140 /* Create a covariance struct.
 141  */
 142 struct covariance *
 143 covariance_1pass_create (size_t n_vars, const struct variable *const *vars,
 144                          const struct variable *weight, enum mv_class exclude,
 145                          bool centered)
 146 {
 147   size_t i;
 148   struct covariance *cov = XZALLOC (struct covariance);
 149
 150   cov->centered = centered;
 151   cov->passes = 1;
 152   cov->state = 0;
 153   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 154
 155   cov->vars = vars;
 156
 157   cov->wv = weight;
 158   cov->n_vars = n_vars;
 159   cov->dim = n_vars;
 160
 161   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 162
 163   for (i = 0; i < n_MOMENTS; ++i)
 164     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 165
 166   cov->exclude = exclude;
 167
 168   cov->n_cm = (n_vars * (n_vars - 1)) / 2;
 169
 170
 171   cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
 172   cov->categoricals = NULL;
 173
 174   return cov;
 175 }
 176
 177 /*
 178   Create a covariance struct for a two-pass algorithm. If categorical
 179   variables are involed, the dimension cannot be know until after the
 180   first data pass, so the actual covariances will not be allocated
 181   until then.
 182  */
 183 struct covariance *
 184 covariance_2pass_create (size_t n_vars, const struct variable *const *vars,
 185                          struct categoricals *cats,
 186                          const struct variable *wv, enum mv_class exclude,
 187                          bool centered)
 188 {
 189   size_t i;
 190   struct covariance *cov = xmalloc (sizeof *cov);
 191
 192   cov->centered = centered;
 193   cov->passes = 2;
 194   cov->state = 0;
 195   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 196
 197   cov->vars = vars;
 198
 199   cov->wv = wv;
 200   cov->n_vars = n_vars;
 201   cov->dim = n_vars;
 202
 203   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 204
 205   for (i = 0; i < n_MOMENTS; ++i)
 206     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 207
 208   cov->exclude = exclude;
 209
 210   cov->n_cm = -1;
 211   cov->cm = NULL;
 212
 213   cov->categoricals = cats;
 214   cov->unnormalised = NULL;
 215
 216   return cov;
 217 }
 218
 219 /* Return an integer, which can be used to index
 220    into COV->cm, to obtain the I, J th element
 221    of the covariance matrix.  If COV->cm does not
 222    contain that element, then a negative value
 223    will be returned.
 224 */
 225 static int
 226 cm_idx (const struct covariance *cov, int i, int j)
 227 {
 228   int as;
 229   const int n2j = cov->dim - 2 - j;
 230   const int nj = cov->dim - 2 ;
 231
 232   assert (i >= 0);
 233   assert (j < cov->dim);
 234
 235   if (i == 0)
 236     return -1;
 237
 238   if (j >= cov->dim - 1)
 239     return -1;
 240
 241   if (i <= j)
 242     return -1 ;
 243
 244   as = nj * (nj + 1) ;
 245   as -= n2j * (n2j + 1) ;
 246   as /= 2;
 247
 248   return i - 1 + as;
 249 }
 250
 251
 252 /*
 253   Returns true iff the variable corresponding to the Ith element of the covariance matrix
 254    has a missing value for case C
 255 */
 256 static bool
 257 is_missing (const struct covariance *cov, int i, const struct ccase *c)
 258 {
 259   const struct variable *var = i < cov->n_vars ?
 260     cov->vars[i] :
 261     categoricals_get_interaction_by_subscript (cov->categoricals, i - cov->n_vars)->vars[0];
 262
 263   const union value *val = case_data (c, var);
 264
 265   return (var_is_value_missing (var, val) & cov->exclude) != 0;
 266 }
 267
 268
 269 static double
 270 get_val (const struct covariance *cov, int i, const struct ccase *c)
 271 {
 272   if (i < cov->n_vars)
 273     {
 274       const struct variable *var = cov->vars[i];
 275
 276       const union value *val = case_data (c, var);
 277
 278       return val->f;
 279     }
 280
 281   return categoricals_get_effects_code_for_case (cov->categoricals, i - cov->n_vars, c);
 282 }
 283
 284 #if 0
 285 void
 286 dump_matrix (const gsl_matrix *m)
 287 {
 288   size_t i, j;
 289
 290   for (i = 0 ; i < m->size1; ++i)
 291     {
 292       for (j = 0 ; j < m->size2; ++j)
 293         printf ("%02f ", gsl_matrix_get (m, i, j));
 294       printf ("\n");
 295     }
 296 }
 297 #endif
 298
 299 /* Call this function for every case in the data set */
 300 void
 301 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 302 {
 303   size_t i, j, m;
 304   const double weight = cov->wv ? case_num (c, cov->wv) : 1.0;
 305
 306   assert (cov->passes == 2);
 307   if (!cov->pass_one_first_case_seen)
 308     {
 309       assert (cov->state == 0);
 310       cov->state = 1;
 311     }
 312
 313   if (cov->categoricals)
 314     categoricals_update (cov->categoricals, c);
 315
 316   for (i = 0 ; i < cov->dim; ++i)
 317     {
 318       double v1 = get_val (cov, i, c);
 319
 320       if (is_missing (cov, i, c))
 321         continue;
 322
 323       for (j = 0 ; j < cov->dim; ++j)
 324         {
 325           double pwr = 1.0;
 326
 327           if (is_missing (cov, j, c))
 328             continue;
 329
 330           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 331             {
 332               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 333
 334               *x += pwr * weight;
 335               pwr *= v1;
 336             }
 337         }
 338     }
 339
 340   cov->pass_one_first_case_seen = true;
 341 }
 342
 343
 344 /* Call this function for every case in the data set */
 345 void
 346 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 347 {
 348   size_t i, j;
 349   const double weight = cov->wv ? case_num (c, cov->wv) : 1.0;
 350
 351   assert (cov->passes == 2);
 352   assert (cov->state >= 1);
 353
 354   if (! cov->pass_two_first_case_seen)
 355     {
 356       size_t m;
 357       assert (cov->state == 1);
 358       cov->state = 2;
 359
 360       if (cov->categoricals)
 361         categoricals_done (cov->categoricals);
 362
 363       cov->dim = cov->n_vars;
 364
 365       if (cov->categoricals)
 366         cov->dim += categoricals_df_total (cov->categoricals);
 367
 368       cov->n_cm = (cov->dim * (cov->dim - 1)) / 2;
 369       cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
 370
 371       /* Grow the moment matrices so that they're large enough to accommodate the
 372          categorical elements */
 373       for (i = 0; i < n_MOMENTS; ++i)
 374         {
 375           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 376         }
 377
 378       /* Populate the moments matrices with the categorical value elements */
 379       for (i = cov->n_vars; i < cov->dim; ++i)
 380         {
 381           for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
 382             {
 383               double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
 384
 385               gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
 386
 387               w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
 388
 389               gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
 390             }
 391         }
 392
 393       /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
 394        now it assumes there are none */
 395       for (m = 0 ; m < n_MOMENTS; ++m)
 396         {
 397           for (i = 0 ; i < cov->dim ; ++i)
 398             {
 399               double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
 400               for (j = cov->n_vars; j < cov->dim; ++j)
 401                 {
 402                   gsl_matrix_set (cov->moments[m], i, j, x);
 403                 }
 404             }
 405         }
 406
 407       /* Divide the means by the number of samples */
 408       for (i = 0; i < cov->dim; ++i)
 409         {
 410           for (j = 0; j < cov->dim; ++j)
 411             {
 412               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 413               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 414             }
 415         }
 416     }
 417
 418   for (i = 0 ; i < cov->dim; ++i)
 419     {
 420       double v1 = get_val (cov, i, c);
 421
 422       if (is_missing (cov, i, c))
 423         continue;
 424
 425       for (j = 0 ; j < cov->dim; ++j)
 426         {
 427           int idx;
 428           double ss ;
 429           double v2 = get_val (cov, j, c);
 430
 431           const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 432
 433           if (is_missing (cov, j, c))
 434             continue;
 435
 436           {
 437             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 438             *x += s;
 439           }
 440
 441           ss =
 442             (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 443             *
 444             (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 445             * weight
 446             ;
 447
 448           idx = cm_idx (cov, i, j);
 449           if (idx >= 0)
 450             {
 451               cov->cm [idx] += ss;
 452             }
 453         }
 454     }
 455
 456   cov->pass_two_first_case_seen = true;
 457 }
 458
 459
 460 /* Call this function for every case in the data set.
 461    After all cases have been passed, call covariance_calculate
 462  */
 463 void
 464 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 465 {
 466   size_t i, j, m;
 467   const double weight = cov->wv ? case_num (c, cov->wv) : 1.0;
 468
 469   assert (cov->passes == 1);
 470
 471   if (!cov->pass_one_first_case_seen)
 472     {
 473       assert (cov->state == 0);
 474       cov->state = 1;
 475     }
 476
 477   for (i = 0 ; i < cov->dim; ++i)
 478     {
 479       const union value *val1 = case_data (c, cov->vars[i]);
 480
 481       if (is_missing (cov, i, c))
 482         continue;
 483
 484       for (j = 0 ; j < cov->dim; ++j)
 485         {
 486           double pwr = 1.0;
 487           int idx;
 488           const union value *val2 = case_data (c, cov->vars[j]);
 489
 490           if (is_missing (cov, j, c))
 491             continue;
 492
 493           idx = cm_idx (cov, i, j);
 494           if (idx >= 0)
 495             {
 496               cov->cm [idx] += val1->f * val2->f * weight;
 497             }
 498
 499           for (m = 0 ; m < n_MOMENTS; ++m)
 500             {
 501               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 502
 503               *x += pwr * weight;
 504               pwr *= val1->f;
 505             }
 506         }
 507     }
 508
 509   cov->pass_one_first_case_seen = true;
 510 }
 511
 512
 513 /*
 514    Allocate and return a gsl_matrix containing the covariances of the
 515    data.
 516 */
 517 static gsl_matrix *
 518 cm_to_gsl (struct covariance *cov)
 519 {
 520   int i, j;
 521   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
 522
 523   /* Copy the non-diagonal elements from cov->cm */
 524   for (j = 0 ; j < cov->dim - 1; ++j)
 525     {
 526       for (i = j+1 ; i < cov->dim; ++i)
 527         {
 528           double x = cov->cm [cm_idx (cov, i, j)];
 529           gsl_matrix_set (m, i, j, x);
 530           gsl_matrix_set (m, j, i, x);
 531         }
 532     }
 533
 534   /* Copy the diagonal elements from cov->moments[2] */
 535   for (j = 0 ; j < cov->dim ; ++j)
 536     {
 537       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 538       gsl_matrix_set (m, j, j, sigma);
 539     }
 540
 541   return m;
 542 }
 543
 544
 545 static gsl_matrix *
 546 covariance_calculate_double_pass (struct covariance *cov)
 547 {
 548   size_t i, j;
 549   for (i = 0 ; i < cov->dim; ++i)
 550     {
 551       for (j = 0 ; j < cov->dim; ++j)
 552         {
 553           int idx;
 554           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 555           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 556
 557           idx = cm_idx (cov, i, j);
 558           if (idx >= 0)
 559             {
 560               x = &cov->cm [idx];
 561               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 562             }
 563         }
 564     }
 565
 566   return  cm_to_gsl (cov);
 567 }
 568
 569 static gsl_matrix *
 570 covariance_calculate_single_pass (struct covariance *cov)
 571 {
 572   size_t i, j;
 573   size_t m;
 574
 575   for (m = 0; m < n_MOMENTS; ++m)
 576     {
 577       /* Divide the moments by the number of samples */
 578       if (m > 0)
 579         {
 580           for (i = 0 ; i < cov->dim; ++i)
 581             {
 582               for (j = 0 ; j < cov->dim; ++j)
 583                 {
 584                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 585                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 586
 587                   if (m == MOMENT_VARIANCE)
 588                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 589                 }
 590             }
 591         }
 592     }
 593
 594   if (cov->centered)
 595     {
 596       /* Centre the moments */
 597       for (j = 0 ; j < cov->dim - 1; ++j)
 598         {
 599           for (i = j + 1 ; i < cov->dim; ++i)
 600             {
 601               double *x = &cov->cm [cm_idx (cov, i, j)];
 602
 603               *x /= gsl_matrix_get (cov->moments[0], i, j);
 604
 605               *x -=
 606                 gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 607                 *
 608                 gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 609             }
 610         }
 611     }
 612
 613   return cm_to_gsl (cov);
 614 }
 615
 616
 617 /* Return a pointer to gsl_matrix containing the pairwise covariances.  The
 618    caller owns the returned matrix and must free it when it is no longer
 619    needed.
 620
 621    Call this function only after all data have been accumulated.  */
 622 gsl_matrix *
 623 covariance_calculate (struct covariance *cov)
 624 {
 625   if (cov->state <= 0)
 626     return NULL;
 627
 628   switch (cov->passes)
 629     {
 630     case 1:
 631       return covariance_calculate_single_pass (cov);
 632       break;
 633     case 2:
 634       return covariance_calculate_double_pass (cov);
 635       break;
 636     default:
 637       NOT_REACHED ();
 638     }
 639 }
 640
 641 /*
 642   Covariance computed without dividing by the sample size.
 643  */
 644 static gsl_matrix *
 645 covariance_calculate_double_pass_unnormalized (struct covariance *cov)
 646 {
 647   return  cm_to_gsl (cov);
 648 }
 649
 650 static gsl_matrix *
 651 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
 652 {
 653   size_t i, j;
 654
 655   if (cov->centered)
 656     {
 657       for (i = 0 ; i < cov->dim; ++i)
 658         {
 659           for (j = 0 ; j < cov->dim; ++j)
 660             {
 661               double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 662               *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 663                 / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 664             }
 665         }
 666
 667       for (j = 0 ; j < cov->dim - 1; ++j)
 668         {
 669           for (i = j + 1 ; i < cov->dim; ++i)
 670             {
 671               double *x = &cov->cm [cm_idx (cov, i, j)];
 672
 673               *x -=
 674                 gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 675                 *
 676                 gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i)
 677                 / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 678             }
 679         }
 680     }
 681
 682   return cm_to_gsl (cov);
 683 }
 684
 685
 686 /* Return a pointer to gsl_matrix containing the pairwise covariances.  The
 687    returned matrix is owned by the structure, and must not be freed.
 688
 689    Call this function only after all data have been accumulated.  */
 690 const gsl_matrix *
 691 covariance_calculate_unnormalized (struct covariance *cov)
 692 {
 693   if (cov->state <= 0)
 694     return NULL;
 695
 696   if (cov->unnormalised != NULL)
 697     return cov->unnormalised;
 698
 699   switch (cov->passes)
 700     {
 701     case 1:
 702       cov->unnormalised =  covariance_calculate_single_pass_unnormalized (cov);
 703       break;
 704     case 2:
 705       cov->unnormalised =  covariance_calculate_double_pass_unnormalized (cov);
 706       break;
 707     default:
 708       NOT_REACHED ();
 709     }
 710
 711   return cov->unnormalised;
 712 }
 713
 714 /* Function to access the categoricals used by COV
 715    The return value is owned by the COV
 716 */
 717 const struct categoricals *
 718 covariance_get_categoricals (const struct covariance *cov)
 719 {
 720   return cov->categoricals;
 721 }
 722
 723
 724 /* Destroy the COV object */
 725 void
 726 covariance_destroy (struct covariance *cov)
 727 {
 728   size_t i;
 729
 730   categoricals_destroy (cov->categoricals);
 731
 732   for (i = 0; i < n_MOMENTS; ++i)
 733     gsl_matrix_free (cov->moments[i]);
 734
 735   gsl_matrix_free (cov->unnormalised);
 736   free (cov->moments);
 737   free (cov->cm);
 738   free (cov);
 739 }
 740
 741 size_t
 742 covariance_dim (const struct covariance * cov)
 743 {
 744   return (cov->dim);
 745 }
 746
 747 \f
 748
 749 /*
 750   Routines to assist debugging.
 751   The following are not thoroughly tested and in certain respects
 752   unreliable.  They should only be
 753   used for aids to development. Not as user accessible code.
 754 */
 755
 756 #include "libpspp/str.h"
 757 #include "output/pivot-table.h"
 758 #include "data/format.h"
 759
 760
 761 /* Create a table which can be populated with the encodings for
 762    the covariance matrix COV */
 763 struct pivot_table *
 764 covariance_dump_enc_header (const struct covariance *cov)
 765 {
 766   struct pivot_table *table = pivot_table_create ("Covariance Encoding");
 767
 768   struct pivot_dimension *factors = pivot_dimension_create (
 769     table, PIVOT_AXIS_COLUMN, "Factor");
 770   for (size_t i = 0 ; i < cov->n_vars; ++i)
 771     pivot_category_create_leaf (factors->root,
 772                                 pivot_value_new_variable (cov->vars[i]));
 773   for (size_t i = 0, n = 0; i < cov->dim - cov->n_vars; n++)
 774     {
 775       const struct interaction *iact =
 776         categoricals_get_interaction_by_subscript (cov->categoricals, i);
 777
 778       struct string str = DS_EMPTY_INITIALIZER;
 779       interaction_to_string (iact, &str);
 780       struct pivot_category *group = pivot_category_create_group__ (
 781         factors->root,
 782         pivot_value_new_user_text_nocopy (ds_steal_cstr (&str)));
 783
 784       int df = categoricals_df (cov->categoricals, n);
 785       for (int j = 0; j < df; j++)
 786         pivot_category_create_leaf_rc (group, pivot_value_new_integer (j),
 787                                        PIVOT_RC_INTEGER);
 788
 789       i += df;
 790     }
 791
 792   struct pivot_dimension *matrix = pivot_dimension_create (
 793     table, PIVOT_AXIS_ROW, "Matrix", "Matrix");
 794   matrix->hide_all_labels = true;
 795
 796   return table;
 797 }
 798
 799
 800 /*
 801   Append table T, which should have been returned by covariance_dump_enc_header
 802   with an entry corresponding to case C for the covariance matrix COV
 803  */
 804 void
 805 covariance_dump_enc (const struct covariance *cov, const struct ccase *c,
 806                      struct pivot_table *table)
 807 {
 808   int row = pivot_category_create_leaf (
 809     table->dimensions[1]->root,
 810     pivot_value_new_integer (table->dimensions[1]->n_leaves));
 811
 812   for (int i = 0 ; i < cov->dim; ++i)
 813     pivot_table_put2 (
 814       table, i, row, pivot_value_new_number (get_val (cov, i, c)));
 815 }