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