pintos-os.org Git - pspp/blob - src/language/data-io/mconvert.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2021 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.h>
  20
  21 #include "data/any-reader.h"
  22 #include "data/any-writer.h"
  23 #include "data/casereader.h"
  24 #include "data/casewriter.h"
  25 #include "data/dataset.h"
  26 #include "data/dictionary.h"
  27 #include "language/data-io/file-handle.h"
  28 #include "language/data-io/matrix-reader.h"
  29 #include "language/lexer/lexer.h"
  30 #include "language/command.h"
  31
  32 #include "gettext.h"
  33 #define _(msgid) gettext (msgid)
  34
  35 int
  36 cmd_mconvert (struct lexer *lexer, struct dataset *ds)
  37 {
  38   bool append = false;
  39   struct file_handle *in = NULL;
  40   struct file_handle *out = NULL;
  41   while (lex_token (lexer) != T_ENDCMD)
  42     {
  43       lex_match (lexer, T_SLASH);
  44
  45       if (lex_match_id (lexer, "APPEND"))
  46         append = true;
  47       else if (lex_match_id (lexer, "REPLACE"))
  48         append = false;
  49       else
  50         {
  51           if (lex_match_id (lexer, "MATRIX"))
  52             lex_match (lexer, T_EQUALS);
  53
  54           struct file_handle **fhp = (lex_match_id (lexer, "IN") ? &in
  55                                       : lex_match_id (lexer, "OUT") ? &out
  56                                       : NULL);
  57           if (!fhp)
  58             {
  59               lex_error_expecting (lexer, "IN", "OUT", "APPEND", "REPLACE");
  60               goto error;
  61             }
  62           if (!lex_force_match (lexer, T_LPAREN))
  63             goto error;
  64
  65           fh_unref (*fhp);
  66           if (lex_match (lexer, T_ASTERISK))
  67             *fhp = NULL;
  68           else
  69             {
  70               *fhp = fh_parse (lexer, FH_REF_FILE, dataset_session (ds));
  71               if (!*fhp)
  72                 goto error;
  73             }
  74
  75           if (!lex_force_match (lexer, T_RPAREN))
  76             goto error;
  77         }
  78     }
  79
  80   if (!in && !dataset_has_source (ds))
  81     {
  82       msg (SE, _("No active file is defined and no external file is "
  83                  "specified on MATRIX=IN."));
  84       goto error;
  85     }
  86
  87   struct dictionary *d;
  88   struct casereader *cr;
  89   if (in)
  90     {
  91       cr = any_reader_open_and_decode (in, NULL, &d, NULL);
  92       if (!cr)
  93         goto error;
  94     }
  95   else
  96     {
  97       d = dict_clone (dataset_dict (ds));
  98       cr = proc_open (ds);
  99     }
 100
 101   struct matrix_reader *mr = matrix_reader_create (d, cr);
 102   if (!mr)
 103     {
 104       casereader_destroy (cr);
 105       dict_unref (d);
 106       if (!in)
 107         proc_commit (ds);
 108       goto error;
 109     }
 110
 111   struct casewriter *cw;
 112   if (out)
 113     {
 114       cw = any_writer_open (out, d);
 115       if (!cw)
 116         {
 117           matrix_reader_destroy (mr);
 118           casereader_destroy (cr);
 119           dict_unref (d);
 120           if (!in)
 121             proc_commit (ds);
 122           goto error;
 123         }
 124     }
 125   else
 126     cw = autopaging_writer_create (dict_get_proto (d));
 127
 128   for (;;)
 129     {
 130       struct matrix_material mm;
 131       struct casereader *group;
 132       if (!matrix_reader_next (&mm, mr, &group))
 133         break;
 134
 135       bool add_corr = mm.cov && !mm.corr;
 136       bool add_cov = mm.corr && !mm.cov && mm.var_matrix;
 137       bool add_stddev = add_corr && !mm.var_matrix;
 138       bool remove_corr = add_cov && !append;
 139       bool remove_cov = add_corr && !append;
 140
 141       struct ccase *model = casereader_peek (group, 0);
 142       for (size_t i = 0; i < mr->n_fvars; i++)
 143         *case_num_rw (model, mr->fvars[i]) = SYSMIS;
 144
 145       for (;;)
 146         {
 147           struct ccase *c = casereader_read (group);
 148           if (!c)
 149             break;
 150
 151           struct substring rowtype = matrix_reader_get_string (c, mr->rowtype);
 152           if ((remove_cov && ss_equals_case (rowtype, ss_cstr ("COV")))
 153               || (remove_corr && ss_equals_case (rowtype, ss_cstr ("CORR"))))
 154             case_unref (c);
 155           else
 156             casewriter_write (cw, c);
 157         }
 158       casereader_destroy (group);
 159
 160       if (add_corr)
 161         {
 162           for (size_t y = 0; y < mr->n_cvars; y++)
 163             {
 164               struct ccase *c = case_clone (model);
 165               for (size_t x = 0; x < mr->n_cvars; x++)
 166                 {
 167                   double d1 = gsl_matrix_get (mm.cov, x, x);
 168                   double d2 = gsl_matrix_get (mm.cov, y, y);
 169                   double cov = gsl_matrix_get (mm.cov, y, x);
 170                   *case_num_rw (c, mr->cvars[x]) = cov / sqrt (d1 * d2);
 171                 }
 172               matrix_reader_set_string (c, mr->rowtype, ss_cstr ("CORR"));
 173               matrix_reader_set_string (c, mr->varname,
 174                                         ss_cstr (var_get_name (mr->cvars[y])));
 175               casewriter_write (cw, c);
 176             }
 177         }
 178
 179       if (add_stddev)
 180         {
 181           struct ccase *c = case_clone (model);
 182           for (size_t x = 0; x < mr->n_cvars; x++)
 183             {
 184               double variance = gsl_matrix_get (mm.cov, x, x);
 185               *case_num_rw (c, mr->cvars[x]) = sqrt (variance);
 186             }
 187           matrix_reader_set_string (c, mr->rowtype, ss_cstr ("STDDEV"));
 188           matrix_reader_set_string (c, mr->varname, ss_empty ());
 189           casewriter_write (cw, c);
 190         }
 191
 192       if (add_cov)
 193         {
 194           for (size_t y = 0; y < mr->n_cvars; y++)
 195             {
 196               struct ccase *c = case_clone (model);
 197               for (size_t x = 0; x < mr->n_cvars; x++)
 198                 {
 199                   double d1 = gsl_matrix_get (mm.var_matrix, x, x);
 200                   double d2 = gsl_matrix_get (mm.var_matrix, y, y);
 201                   double corr = gsl_matrix_get (mm.corr, y, x);
 202                   *case_num_rw (c, mr->cvars[x]) = corr * sqrt (d1 * d2);
 203                 }
 204               matrix_reader_set_string (c, mr->rowtype, ss_cstr ("COV"));
 205               matrix_reader_set_string (c, mr->varname,
 206                                         ss_cstr (var_get_name (mr->cvars[y])));
 207               casewriter_write (cw, c);
 208             }
 209         }
 210
 211       case_unref (model);
 212       matrix_material_uninit (&mm);
 213     }
 214
 215   matrix_reader_destroy (mr);
 216   if (!in)
 217     proc_commit (ds);
 218   if (out)
 219     casewriter_destroy (cw);
 220   else
 221     {
 222       dataset_set_dict (ds, dict_ref (d));
 223       dataset_set_source (ds, casewriter_make_reader (cw));
 224     }
 225
 226   fh_unref (in);
 227   fh_unref (out);
 228   dict_unref (d);
 229   return CMD_SUCCESS;
 230
 231 error:
 232   fh_unref (in);
 233   fh_unref (out);
 234   return CMD_FAILURE;
 235 }
 236