refactoring
[pspp] / 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