ML to FIML conversion now happens in front-end.
[openmx:openmx.git] / src / omxMatrix.c
1 /*
2  *  Copyright 2007-2009 The OpenMx Project
3  *
4  *  Licensed under the Apache License, Version 2.0 (the "License");
5  *  you may not use this file except in compliance with the License.
6  *  You may obtain a copy of the License at
7  *
8  *       http://www.apache.org/licenses/LICENSE-2.0
9  *
10  *  Unless required by applicable law or agreed to in writing, software
11  *  distributed under the License is distributed on an "AS IS" BASIS,
12  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  *  See the License for the specific language governing permissions and
14  *  limitations under the License.
15  */
16
17 /***********************************************************
18 *
19 *  omxMatrix.cc
20 *
21 *  Created: Timothy R. Brick    Date: 2008-11-13 12:33:06
22 *
23 *       Contains code for the omxMatrix class
24 *   omxDataMatrices hold necessary information to simplify
25 *       dealings between the OpenMX back end and BLAS.
26 *
27 **********************************************************/
28 #include "omxMatrix.h"
29
30 const char omxMatrixMajorityList[3] = "Tn";             // BLAS Column Majority.
31
32 void omxPrintMatrix(omxMatrix *source, char* header) {
33         int j, k;
34
35         Rprintf("%s: (%d x %d) [%s-major]\n", header, source->rows, source->cols, (source->colMajor?"col":"row"));
36         if(OMX_DEBUG) {Rprintf("Matrix Printing is at %0x\n", source);}
37
38         if(source->colMajor) {
39                 for(j = 0; j < source->rows; j++) {
40                         for(k = 0; k < source->cols; k++) {
41                                 Rprintf("\t%3.6f", source->data[k*source->rows+j]);
42                         }
43                         Rprintf("\n");
44                 }
45         } else {
46                 for(j = 0; j < source->cols; j++) {
47                         for(k = 0; k < source->rows; k++) {
48                                 Rprintf("\t%3.6f", source->data[k*source->cols+j]);
49                         }
50                         Rprintf("\n");
51                 }
52         }
53 }
54
55 omxMatrix* omxInitMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
56
57         if(om == NULL) om = (omxMatrix*) R_alloc(1, sizeof(omxMatrix));
58         if(OMX_DEBUG) { Rprintf("Initializing 0x%0x to (%d, %d) with state at 0x%x.\n", om, nrows, ncols, os); }
59
60         om->rows = nrows;
61         om->cols = ncols;
62         om->colMajor = (isColMajor?1:0);
63
64         om->originalRows = om->rows;
65         om->originalCols = om->cols;
66         om->originalColMajor=om->colMajor;
67
68         if(om->rows == 0 || om->cols == 0) {
69                 om->data = NULL;
70                 om->localData = FALSE;
71         } else {
72                 om->data = (double*) Calloc(nrows * ncols, double);
73                 om->localData = TRUE;
74         }
75
76         om->populateFrom = NULL;
77         om->populateToCol = NULL;
78         om->populateToRow = NULL;
79         om->numPopulateLocations = 0;
80
81         om->aliasedPtr = NULL;
82         om->algebra = NULL;
83         om->objective = NULL;
84
85         om->currentState = os;
86         om->lastCompute = -1;
87         om->lastRow = -1;
88
89         omxComputeMatrix(om);
90
91         return om;
92
93 }
94
95 void omxCopyMatrix(omxMatrix *dest, omxMatrix *orig) {
96         /* Duplicate a matrix.  NOTE: Matrix maintains its algebra bindings. */
97
98         if(OMX_DEBUG) { Rprintf("omxCopyMatrix"); }
99
100         omxFreeMatrixData(dest);
101
102         dest->rows = orig->rows;
103         dest->cols = orig->cols;
104         dest->colMajor = orig->colMajor;
105         dest->originalRows = dest->rows;
106         dest->originalCols = dest->cols;
107         dest->originalColMajor = dest->colMajor;
108         dest->currentState = orig->currentState;
109         dest->lastCompute = orig->lastCompute;
110         dest->lastRow = orig->lastRow;
111
112         if(dest->rows == 0 || dest->cols == 0) {
113                 dest->data = NULL;
114                 dest->localData=FALSE;
115         } else {
116                 dest->data = (double*) Calloc(dest->rows * dest->cols, double);
117                 memcpy(dest->data, orig->data, dest->rows * dest->cols * sizeof(double));
118                 dest->localData = TRUE;
119         }
120
121         dest->aliasedPtr = NULL;
122
123         omxComputeMatrix(dest);
124
125 }
126
127 void omxAliasMatrix(omxMatrix *dest, omxMatrix *src) {
128         omxCopyMatrix(dest, src);
129         dest->aliasedPtr = src;                                 // Alias now follows back matrix precisely.
130         dest->algebra = NULL;                                   // Have to look at how this effect interacts with populating
131         dest->objective = NULL;                                 //              matrix values to other locations.
132 }
133
134 void omxFreeMatrixData(omxMatrix * om) {
135
136         if(om->localData && om->data != NULL) {
137                 if(OMX_DEBUG) { Rprintf("Freeing 0x%0x. Localdata = %d.\n", om->data, om->localData); }
138                 Free(om->data);
139                 om->data = NULL;
140                 om->localData = FALSE;
141         }
142
143 }
144
145 void omxFreeAllMatrixData(omxMatrix *om) {
146
147         if(OMX_DEBUG) { Rprintf("Freeing 0x%0x with data = %0x and algebra %0x.\n", om, om->data, om->algebra); }
148
149         if(om->localData && om->data != NULL) {
150                 Free(om->data);
151                 om->data = NULL;
152                 om->localData = FALSE;
153         }
154
155         if(om->algebra != NULL) {
156                 omxFreeAlgebraArgs(om->algebra);
157                 om->algebra = NULL;
158         }
159
160         if(om->objective != NULL) {
161                 omxFreeObjectiveArgs(om->objective);
162                 om->objective = NULL;
163         }
164
165 }
166
167 void omxResizeMatrix(omxMatrix *om, int nrows, int ncols, unsigned short keepMemory) {
168         // Always Recompute() before you Resize().
169         if(OMX_DEBUG) { Rprintf("Resizing matrix from (%d, %d) to (%d, %d) (keepMemory: %d)", om->rows, om->cols, nrows, ncols, keepMemory); }
170         if(keepMemory == FALSE) {
171                 if(OMX_DEBUG) { Rprintf(" and regenerating memory to do it"); }
172                 omxFreeMatrixData(om);
173                 om->data = (double*) Calloc(nrows * ncols, double);
174                 om->localData = TRUE;
175         } else if(om->originalRows * om->originalCols < nrows * ncols) {
176                 warning("Upsizing an existing matrix may cause undefined behavior.\n"); // TODO: Define this behavior?
177         }
178
179         if(OMX_DEBUG) { Rprintf(".\n"); }
180         om->rows = nrows;
181         om->cols = ncols;
182         if(keepMemory == FALSE) {
183                 om->originalRows = om->rows;
184                 om->originalCols = om->cols;
185         }
186
187         omxComputeMatrix(om);
188 }
189
190 void omxResetAliasedMatrix(omxMatrix *om) {
191         om->rows = om->originalRows;
192         om->cols = om->originalCols;
193         if(om->aliasedPtr != NULL) {
194                 memcpy(om->data, om->aliasedPtr->data, om->rows*om->cols*sizeof(double));
195                 om->colMajor = om->aliasedPtr->colMajor;
196         }
197         omxComputeMatrix(om);
198 }
199
200 void omxComputeMatrix(omxMatrix *om) {
201
202 //      if(OMX_DEBUG) { Rprintf("Matrix compute: 0x%0x, 0x%0x, %d.\n", om, om->currentState, om->colMajor); }
203         om->majority = &(omxMatrixMajorityList[(om->colMajor?1:0)]);
204         om->minority = &(omxMatrixMajorityList[(om->colMajor?0:1)]);
205         om->leading = (om->colMajor?om->rows:om->cols);
206         om->lagging = (om->colMajor?om->cols:om->rows);
207
208         for(int i = 0; i < om->numPopulateLocations; i++) {
209                 omxRecompute(om->populateFrom[i]);                              // Make sure it's up to date
210                 omxSetMatrixElement(om, om->populateToRow[i], om->populateToCol[i], om->populateFrom[i]->data[0]);
211                 // And then fill in the details.  Use the Helper here in case of transposition/downsampling.
212         }
213
214         om->isDirty = FALSE;
215         om->lastCompute = om->currentState->computeCount;
216         om->lastRow = om->currentState->currentRow;
217
218 }
219
220 double* omxLocationOfMatrixElement(omxMatrix *om, int row, int col) {
221         int index = 0;
222         if(om->colMajor) {
223                 index = col * om->rows + row;
224         } else {
225                 index = row * om->cols + col;
226         }
227         return om->data + index;
228 }
229
230 double omxVectorElement(omxMatrix *om, int index) {
231         if(index < om->rows * om->cols) {
232                 return om->data[index];
233         }
234 }
235
236 double omxAliasedMatrixElement(omxMatrix *om, int row, int col) {
237         int index = 0;
238         if(row >= om->originalRows || col >= om->originalCols) {
239                 char errstr[250];
240                 sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.", row, col, om->originalRows, om->originalCols);
241                 error(errstr);
242         }
243         if(om->colMajor) {
244                 index = col * om->originalRows + row;
245         } else {
246                 index = row * om->originalCols + col;
247         }
248         return om->data[index];
249
250 }
251
252 double omxMatrixElement(omxMatrix *om, int row, int col) {
253         int index = 0;
254         if(row >= om->rows || col >= om->cols) {
255                 char errstr[250];
256                 sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.", row, col, om->rows, om->cols);
257                 error(errstr);
258         }
259         if(om->colMajor) {
260                 index = col * om->rows + row;
261         } else {
262                 index = row * om->cols + col;
263         }
264         return om->data[index];
265 }
266
267 void omxSetMatrixElement(omxMatrix *om, int row, int col, double value) {
268         int index = 0;
269         if(om->colMajor) {
270                 index = col * om->rows + row;
271         } else {
272                 index = row * om->cols + col;
273         }
274         om->data[index] = value;
275 }
276
277 void omxMarkDirty(omxMatrix *om) { om->isDirty = TRUE; }
278
279 unsigned short omxMatrixNeedsUpdate(omxMatrix *om) {
280
281         for(int i = 0; i < om->numPopulateLocations; i++) {
282                 if(omxNeedsUpdate(om->populateFrom[i])) return TRUE;    // Make sure it's up to date
283         }
284
285 };
286
287 omxMatrix* omxNewMatrixFromMxMatrix(SEXP mxMatrix, omxState* state) {
288 /* Populates the fields of a omxMatrix with details from an R Matrix. */
289
290         SEXP className;
291         SEXP matrixDims;
292         SEXP matrix = mxMatrix;
293         int* dimList;
294         unsigned short int isMxMatrix = FALSE;
295
296         omxMatrix *om = NULL;
297         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
298
299         if(OMX_DEBUG) { Rprintf("Filling omxMatrix from R matrix.\n"); }
300
301         if(!isMatrix(mxMatrix) && !isVector(mxMatrix)) { // Sanity Check
302                 if(OMX_DEBUG) { Rprintf("R matrix is an object of some sort.\n"); }
303                 if(inherits(mxMatrix, "MxMatrix")) {
304                         if(OMX_DEBUG) { Rprintf("R matrix is Mx Matrix.  Processing.\n"); }
305                         PROTECT(matrix = GET_SLOT(mxMatrix,  install("values")));
306                         isMxMatrix = TRUE; // So we remember to unprotect.
307                 } else {
308                         error("Recieved unknown matrix type.");
309                 }
310         }
311
312         om->data = REAL(matrix);        // TODO: Class-check first?
313
314         if(isMatrix(matrix)) {
315                 PROTECT(matrixDims = getAttrib(matrix, R_DimSymbol));
316                 dimList = INTEGER(matrixDims);
317                 om->rows = dimList[0];
318                 om->cols = dimList[1];
319                 UNPROTECT(1);   // MatrixDims
320         } else if (isVector(matrix)) {          // If it's a vector, assume it's a row vector. BLAS doesn't care.
321                 if(OMX_DEBUG) { Rprintf("Vector discovered.  Assuming rowity.\n"); }
322                 om->rows = 1;
323                 om->cols = length(matrix);
324         }
325         if(OMX_DEBUG) { Rprintf("Matrix connected to (%d, %d) mxMatrix.\n", om->rows, om->cols); }
326
327         om->localData = FALSE;
328         om->colMajor = TRUE;
329         om->originalRows = om->rows;
330         om->originalCols = om->cols;
331         om->originalColMajor = TRUE;
332         om->aliasedPtr = NULL;
333         om->algebra = NULL;
334         om->objective = NULL;
335         om->currentState = state;
336         om->lastCompute = -1;
337         om->lastRow = -1;
338
339         if(OMX_DEBUG) { Rprintf("Pre-compute call.\n");}
340         omxComputeMatrix(om);
341         if(OMX_DEBUG) { Rprintf("Post-compute call.\n");}
342
343         if(OMX_DEBUG) {
344                 omxPrintMatrix(om, "Finished importing matrix");
345         }
346
347         if(isMxMatrix) {
348                 UNPROTECT(1); // matrix
349         }
350
351         return om;
352 }
353
354 void omxProcessMatrixPopulationList(omxMatrix* matrix, SEXP matStruct) {
355
356         if(OMX_DEBUG) { Rprintf("Processing Population List: %d elements.\n", length(matStruct) - 1); }
357         SEXP subList;
358         SEXP matLoc, xLoc, yLoc;
359
360         if(length(matStruct) > 1) {
361                 int numPopLocs = length(matStruct) - 1;
362                 matrix->numPopulateLocations = numPopLocs;
363                 matrix->populateFrom = (omxMatrix**)R_alloc(numPopLocs, sizeof(omxMatrix*));
364                 matrix->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
365                 matrix->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
366         }
367
368         for(int i = 0; i < length(matStruct)-1; i++) {
369                 PROTECT(subList = AS_INTEGER(VECTOR_ELT(matStruct, i+1)));
370
371                 int* locations = INTEGER(subList);
372                 int loc = locations[2];
373                 if(OMX_DEBUG) { Rprintf("."); } //:::
374                 if(loc < 0) {                   // NOTE: This duplicates some of the functionality of NewMatrixFromMxIndex
375                         matrix->populateFrom[i] = matrix->currentState->matrixList[(~loc)];
376                 } else {
377                         matrix->populateFrom[i] = matrix->currentState->algebraList[(loc)];
378                 }
379
380                 matrix->populateToRow[i] = locations[0];
381                 matrix->populateToCol[i] = locations[1];
382
383                 UNPROTECT(1); // subList
384         }
385 }
386
387 void omxRemoveRowsAndColumns(omxMatrix *om, int numRowsRemoved, int numColsRemoved, int rowsRemoved[], int colsRemoved[])
388 {
389 //      if(OMX_DEBUG) { Rprintf("Removing %d rows and %d columns from 0x%0x.\n", numRowsRemoved, numColsRemoved, om);}
390
391         if(numRowsRemoved < 1 && numColsRemoved < 1) { return; }
392
393         int oldRows, oldCols;
394
395         if(om->aliasedPtr == NULL) {
396                 if(om->originalRows == 0 || om->originalCols == 0) {
397                         om->originalRows = om->rows;
398                         om->originalCols = om->cols;
399                 }
400                 oldRows = om->originalRows;
401                 oldCols = om->originalCols;
402         } else {
403                 oldRows = om->aliasedPtr->rows;
404                 oldCols = om->aliasedPtr->cols;
405         }
406
407         int numCols = 0;
408         int nextCol = 0;
409         int nextRow = 0;
410         int j,k;
411
412         if(om->rows > om->originalRows || om->cols > om->originalCols) {        // sanity check.
413                 error("Aliased Matrix is too small for alias.");
414         }
415
416         om->rows = oldRows - numRowsRemoved;
417         om->cols = oldCols - numColsRemoved;
418
419         // Note:  This really aught to be done using a matrix multiply.  Why isn't it?
420         for(int j = 0; j < oldCols; j++) {
421                 if(OMX_DEBUG) { Rprintf("Handling column %d/%d...", j, oldCols);}
422                 if(colsRemoved[j]) {
423                         if(OMX_DEBUG) { Rprintf("Removed.\n");}
424                         continue;
425                 } else {
426                         nextRow = 0;
427                         if(OMX_DEBUG) { Rprintf("Rows (max %d): ", oldRows); }
428                         for(int k = 0; k < oldRows; k++) {
429                                 if(rowsRemoved[k]) {
430                                         if(OMX_DEBUG) { Rprintf("%d removed....", k);}
431                                         continue;
432                                 } else {
433                                         if(OMX_DEBUG) { Rprintf("%d kept....", k);}
434                                         if(om->aliasedPtr == NULL) {
435                                                 omxSetMatrixElement(om, nextRow, nextCol, omxAliasedMatrixElement(om, k, j));
436                                         } else {
437                                                 omxSetMatrixElement(om, nextRow, nextCol, omxMatrixElement(om->aliasedPtr, k,  j));
438                                         }
439                                         if(OMX_DEBUG) {
440                                                 omxPrint(om, "Now Reads: (:::DEBUG:::)");
441                                         }
442                                         nextRow++;
443                                 }
444                         }
445                         if(OMX_DEBUG) { Rprintf("\n");}
446                         nextCol++;
447                 }
448         }
449
450         omxComputeMatrix(om);
451 }
452
453 /* Function wrappers that switch based on inclusion of algebras */
454 void omxPrint(omxMatrix *source, char* d) {                                     // Pretty-print a (small) matrix
455         if(source->algebra != NULL) omxAlgebraPrint(source->algebra, d);
456         else if(source->objective != NULL) omxObjectivePrint(source->objective, d);
457         else omxPrintMatrix(source, d);
458 }
459
460 unsigned short omxNeedsUpdate(omxMatrix *matrix) {
461         /* Simplest update check: If we're dirty or haven't computed this cycle (iteration or row), we need to. */
462         if(OMX_DEBUG) {Rprintf("MatrixNeedsUpdate?");}
463         if(matrix->isDirty) return TRUE;
464         if(matrix->lastCompute < matrix->currentState->computeCount) return TRUE;       // No need to check args if oa's dirty.
465         if(matrix->lastRow < matrix->currentState->currentRow) return TRUE;                     // Ditto.
466
467         if(matrix->algebra != NULL) return omxAlgebraNeedsUpdate(matrix->algebra);
468         else if(matrix->objective != NULL) return omxObjectiveNeedsUpdate(matrix->objective);
469         else return omxMatrixNeedsUpdate(matrix);
470
471 }
472
473 void inline omxRecompute(omxMatrix *matrix) {
474         if(!omxNeedsUpdate(matrix)) return;
475         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
476         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
477         else omxComputeMatrix(matrix);
478 }
479
480 void inline omxCompute(omxMatrix *matrix) {
481         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
482         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
483         else omxComputeMatrix(matrix);
484 }