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