Added 'name=' argument to omxSetParameters()
[openmx:openmx.git] / src / omxMatrix.c
1 /*
2  *  Copyright 2007-2012 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 // forward declarations
31 omxMatrix* fillMatrixHelperFunction(omxMatrix* om, SEXP matrix, omxState* state,
32         unsigned short hasMatrixNumber, int matrixNumber);
33
34 const char omxMatrixMajorityList[3] = "Tn";             // BLAS Column Majority.
35
36 void omxPrintMatrix(omxMatrix *source, char* header) {
37         int j, k;
38
39         Rprintf("%s: (%d x %d) [%s-major]\n", header, source->rows, source->cols, (source->colMajor?"col":"row"));
40         if(OMX_DEBUG_MATRIX) {Rprintf("Matrix Printing is at %0x\n", source);}
41
42         if(source->colMajor) {
43                 for(j = 0; j < source->rows; j++) {
44                         for(k = 0; k < source->cols; k++) {
45                                 Rprintf("\t%3.6f", source->data[k*source->rows+j]);
46                         }
47                         Rprintf("\n");
48                 }
49         } else {
50                 for(j = 0; j < source->cols; j++) {
51                         for(k = 0; k < source->rows; k++) {
52                                 Rprintf("\t%3.6f", source->data[k*source->cols+j]);
53                         }
54                         Rprintf("\n");
55                 }
56         }
57 }
58
59 omxMatrix* omxInitMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
60
61         if(om == NULL) om = (omxMatrix*) R_alloc(1, sizeof(omxMatrix));
62         if(OMX_DEBUG_MATRIX) { Rprintf("Initializing matrix 0x%0x to (%d, %d) with state at 0x%x.\n", om, nrows, ncols, os); }
63
64         om->hasMatrixNumber = 0;
65         om->rows = nrows;
66         om->cols = ncols;
67         om->colMajor = (isColMajor?1:0);
68
69         om->originalRows = om->rows;
70         om->originalCols = om->cols;
71         om->originalColMajor=om->colMajor;
72
73         if(om->rows == 0 || om->cols == 0) {
74                 om->data = NULL;
75                 om->localData = FALSE;
76         } else {
77                 om->data = (double*) Calloc(nrows * ncols, double);
78                 om->localData = TRUE;
79         }
80
81         om->populateFrom = NULL;
82         om->populateFromCol = NULL;
83         om->populateFromRow = NULL;
84         om->populateToCol = NULL;
85         om->populateToRow = NULL;
86
87         om->numPopulateLocations = 0;
88
89         om->aliasedPtr = NULL;
90         om->algebra = NULL;
91         om->objective = NULL;
92
93         om->currentState = os;
94         om->lastCompute = -2;
95         om->lastRow = -2;
96         om->isTemporary = FALSE;
97
98         omxMatrixCompute(om);
99
100         return om;
101
102 }
103
104 omxMatrix* omxInitTemporaryMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
105
106         if(om == NULL) {
107                 om = (omxMatrix*) Calloc(1, omxMatrix);
108         }
109
110         om = omxInitMatrix(om, nrows, ncols, isColMajor, os);
111         om->isTemporary = TRUE;
112         
113         return(om);
114
115 }
116
117 void omxUpdateMatrixHelper(omxMatrix *dest, omxMatrix *orig) {
118
119         dest->rows = orig->rows;
120         dest->cols = orig->cols;
121         dest->colMajor = orig->colMajor;
122         dest->originalRows = dest->rows;
123         dest->originalCols = dest->cols;
124         dest->originalColMajor = dest->colMajor;
125         dest->lastCompute = orig->lastCompute;
126         dest->lastRow = orig->lastRow;
127         dest->localData = orig->localData;
128
129         if ((dest->data == NULL) && (dest->rows > 0) && (dest->cols > 0)) {
130                 dest->data = (double*) Calloc(dest->rows * dest->cols, double);
131         }
132
133         memcpy(dest->data, orig->data, dest->rows * dest->cols * sizeof(double));
134
135         if (dest->aliasedPtr != NULL && orig->aliasedPtr != NULL) {
136                 omxUpdateMatrix(dest->aliasedPtr, orig->aliasedPtr);
137         }
138
139         omxMatrixCompute(dest);
140 }
141
142 void omxUpdateMatrix(omxMatrix *dest, omxMatrix *orig) {
143         omxUpdateMatrixHelper(dest, orig);
144
145         if (orig->algebra) {
146                 omxUpdateAlgebra(dest->algebra, orig->algebra);
147         } else if (orig->objective) {
148                 omxUpdateObjectiveFunction(dest->objective, orig->objective);
149         }
150
151 }
152
153 void omxCopyMatrix(omxMatrix *dest, omxMatrix *orig) {
154         /* Copy a matrix.  NOTE: Matrix maintains its algebra bindings. */
155
156         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("omxCopyMatrix"); }
157
158         int regenerateMemory = TRUE;
159         int numPopLocs = orig->numPopulateLocations;
160
161         if(dest->localData && (dest->originalRows == orig->rows && dest->originalCols == orig->cols)) {
162                 regenerateMemory = FALSE;                               // If it's local data and the right size, we can keep memory.
163         }
164
165         dest->rows = orig->rows;
166         dest->cols = orig->cols;
167         dest->colMajor = orig->colMajor;
168         dest->originalRows = dest->rows;
169         dest->originalCols = dest->cols;
170         dest->originalColMajor = dest->colMajor;
171         dest->lastCompute = orig->lastCompute;
172         dest->lastRow = orig->lastRow;
173
174         dest->numPopulateLocations = numPopLocs;
175         if (numPopLocs > 0) {
176                 dest->populateFrom = (int*)R_alloc(numPopLocs, sizeof(int));
177                 dest->populateFromRow = (int*)R_alloc(numPopLocs, sizeof(int));
178                 dest->populateFromCol = (int*)R_alloc(numPopLocs, sizeof(int));
179                 dest->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
180                 dest->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
181                 
182                 memcpy(dest->populateFrom, orig->populateFrom, numPopLocs * sizeof(int));
183                 memcpy(dest->populateFromRow, orig->populateFromRow, numPopLocs * sizeof(int));
184                 memcpy(dest->populateFromCol, orig->populateFromCol, numPopLocs * sizeof(int));
185                 memcpy(dest->populateToRow, orig->populateToRow, numPopLocs * sizeof(int));
186                 memcpy(dest->populateToCol, orig->populateToCol, numPopLocs * sizeof(int));
187         }
188
189         if(dest->rows == 0 || dest->cols == 0) {
190                 omxFreeMatrixData(dest);
191                 dest->data = NULL;
192                 dest->localData=FALSE;
193         } else {
194                 if(regenerateMemory) {
195                         omxFreeMatrixData(dest);                                                                                        // Free and regenerate memory
196                         dest->data = (double*) Calloc(dest->rows * dest->cols, double);
197                 }
198                 memcpy(dest->data, orig->data, dest->rows * dest->cols * sizeof(double));
199                 dest->localData = TRUE;
200         }
201
202         dest->aliasedPtr = NULL;
203
204         omxMatrixCompute(dest);
205 }
206
207 void omxAliasMatrix(omxMatrix *dest, omxMatrix *src) {
208         omxCopyMatrix(dest, src);
209         dest->aliasedPtr = src;                                 // Alias now follows back matrix precisely.
210 }
211
212 void omxFreeMatrixData(omxMatrix * om) {
213
214         if(om->localData && om->data != NULL) {
215                 if(OMX_DEBUG_MATRIX) { Rprintf("Freeing matrix at 0x%0x. Localdata = %d.\n", om->data, om->localData); }
216                 Free(om->data);
217                 om->data = NULL;
218                 om->localData = FALSE;
219         }
220
221 }
222
223 void omxFreeAllMatrixData(omxMatrix *om) {
224     
225     if(om == NULL) return;
226
227         if(OMX_DEBUG) { 
228             Rprintf("Freeing matrix at 0x%0x with data = 0x%x, algebra 0x%x, and objective 0x%x.\n", 
229                 om, om->data, om->algebra); 
230         }
231
232         if(om->localData && om->data != NULL) {
233                 Free(om->data);
234                 om->data = NULL;
235                 om->localData = FALSE;
236         }
237
238         if(om->algebra != NULL) {
239                 omxFreeAlgebraArgs(om->algebra);
240                 om->algebra = NULL;
241         }
242
243         if(om->objective != NULL) {
244                 omxFreeObjectiveArgs(om->objective);
245                 om->objective = NULL;
246         }
247         
248         if(om->isTemporary) {
249                 Free(om);
250                 om = NULL;
251         }
252
253 }
254
255 void omxZeroByZeroMatrix(omxMatrix *om) {
256         if (om->rows > 0 || om->cols > 0) {
257                 omxResizeMatrix(om, 0, 0, FALSE);
258         }
259 }
260
261 omxMatrix* omxNewIdentityMatrix(int nrows, omxState* state) {
262         omxMatrix* newMat = NULL;
263         int l,k;
264
265         newMat = omxInitMatrix(newMat, nrows, nrows, FALSE, state);
266         for(k =0; k < newMat->rows; k++) {
267                 for(l = 0; l < newMat->cols; l++) {
268                         if(l == k) {
269                                 omxSetMatrixElement(newMat, k, l, 1);
270                         } else {
271                                 omxSetMatrixElement(newMat, k, l, 0);
272                         }
273                 }
274         }
275         return newMat;
276 }
277
278 omxMatrix* omxDuplicateMatrix(omxMatrix* src, omxState* newState) {
279         omxMatrix* newMat;
280     
281         if(src == NULL) return NULL;
282         newMat = omxInitMatrix(NULL, src->rows, src->cols, FALSE, newState);
283         omxCopyMatrix(newMat, src);
284         newMat->hasMatrixNumber = src->hasMatrixNumber;
285         newMat->matrixNumber    = src->matrixNumber;
286     
287     return newMat;    
288 }
289
290 void omxResizeMatrix(omxMatrix *om, int nrows, int ncols, unsigned short keepMemory) {
291         // Always Recompute() before you Resize().
292         if(OMX_DEBUG_MATRIX) { 
293                 Rprintf("Resizing matrix from (%d, %d) to (%d, %d) (keepMemory: %d)", 
294                         om->rows, om->cols, 
295                         nrows, ncols, keepMemory);
296         }
297         if((keepMemory == FALSE) && (om->rows != nrows || om->cols != ncols)) {
298                 if(OMX_DEBUG_MATRIX) { Rprintf(" and regenerating memory to do it"); }
299                 omxFreeMatrixData(om);
300                 om->data = (double*) Calloc(nrows * ncols, double);
301                 om->localData = TRUE;
302         } else if(om->originalRows * om->originalCols < nrows * ncols) {
303                 warning("Upsizing an existing matrix may cause undefined behavior.\n"); // TODO: Define this behavior?
304         }
305
306         if(OMX_DEBUG_MATRIX) { Rprintf(".\n"); }
307         om->rows = nrows;
308         om->cols = ncols;
309         if(keepMemory == FALSE) {
310                 om->originalRows = om->rows;
311                 om->originalCols = om->cols;
312         }
313
314         omxMatrixCompute(om);
315 }
316
317 void omxResetAliasedMatrix(omxMatrix *om) {
318         om->rows = om->originalRows;
319         om->cols = om->originalCols;
320         if(om->aliasedPtr != NULL) {
321                 memcpy(om->data, om->aliasedPtr->data, om->rows*om->cols*sizeof(double));
322                 om->colMajor = om->aliasedPtr->colMajor;
323         }
324         omxMatrixCompute(om);
325 }
326
327 void omxMatrixCompute(omxMatrix *om) {
328
329         if(OMX_DEBUG_MATRIX) { Rprintf("Matrix compute: 0x%0x, 0x%0x, algebra: 0x%x.\n", om, om->currentState, om->algebra); }
330         om->majority = &(omxMatrixMajorityList[(om->colMajor?1:0)]);
331         om->minority = &(omxMatrixMajorityList[(om->colMajor?0:1)]);
332         om->leading = (om->colMajor?om->rows:om->cols);
333         om->lagging = (om->colMajor?om->cols:om->rows);
334
335         for(int i = 0; i < om->numPopulateLocations; i++) {
336                 int index = om->populateFrom[i];
337                 omxMatrix* sourceMatrix;
338                 if (index < 0) {
339                         sourceMatrix = om->currentState->matrixList[~index];
340                 } else {
341                         sourceMatrix = om->currentState->algebraList[index];
342                 }
343                 if (sourceMatrix != NULL) {
344                         omxRecompute(sourceMatrix);                             // Make sure it's up to date
345                         double value = omxMatrixElement(sourceMatrix, om->populateFromRow[i], om->populateFromCol[i]);
346                         omxSetMatrixElement(om, om->populateToRow[i], om->populateToCol[i], value);
347                         // And then fill in the details.  Use the Helper here in case of transposition/downsampling.
348                 }
349         }
350
351         om->isDirty = FALSE;
352         om->lastCompute = om->currentState->computeCount;
353         om->lastRow = om->currentState->currentRow;
354
355 }
356
357 double* omxLocationOfMatrixElement(omxMatrix *om, int row, int col) {
358         int index = 0;
359         if(om->colMajor) {
360                 index = col * om->rows + row;
361         } else {
362                 index = row * om->cols + col;
363         }
364         return om->data + index;
365 }
366
367 void vectorElementError(int index, int numrow, int numcol) {
368         char *errstr = calloc(250, sizeof(char));
369         sprintf(errstr, "Requested improper index (%d) from (%d, %d) vector.", 
370                 index, numrow, numcol);
371         error(errstr);
372         free(errstr);
373 }
374
375 void setMatrixError(omxMatrix *om, int row, int col, int numrow, int numcol) {
376         char *errstr = calloc(250, sizeof(char));
377         static const char *matrixString = "matrix";
378         static const char *algebraString = "algebra";
379         static const char *objectiveString = "objective";
380         const char *typeString;
381         if (om->algebra != NULL) {
382                 typeString = algebraString;
383         } else if (om->objective != NULL) {
384                 typeString = objectiveString;
385         } else {
386                 typeString = matrixString;
387         }
388         if (om->name == NULL) {
389                 sprintf(errstr, "Attempted to set row and column (%d, %d) in %s with dimensions %d x %d.", 
390                         row, col, typeString, numrow, numcol);
391         } else {
392                 sprintf(errstr, "Attempted to set row and column (%d, %d) in %s \"%s\" with dimensions %d x %d.", 
393                         row, col, typeString, om->name, numrow, numcol);
394         }
395         error(errstr);
396         free(errstr);
397 }
398
399 void matrixElementError(int row, int col, int numrow, int numcol) {
400         char *errstr = calloc(250, sizeof(char));
401         sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.",
402                 row, col, numrow, numcol);
403         error(errstr);
404         free(errstr);
405 }
406
407 void setVectorError(int index, int numrow, int numcol) {
408         char *errstr = calloc(250, sizeof(char));
409         sprintf(errstr, "Setting improper index (%d) from (%d, %d) vector.", 
410                 index, numrow, numcol);
411         error(errstr);
412         free(errstr);
413 }
414
415 double omxAliasedMatrixElement(omxMatrix *om, int row, int col) {
416         int index = 0;
417         if(row >= om->originalRows || col >= om->originalCols) {
418                 char *errstr = calloc(250, sizeof(char));
419                 sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.", row+1, col+1, om->originalRows, om->originalCols);
420                 error(errstr);
421                 free(errstr);
422         return (NA_REAL);
423         }
424         if(om->colMajor) {
425                 index = col * om->originalRows + row;
426         } else {
427                 index = row * om->originalCols + col;
428         }
429         return om->data[index];
430 }
431
432 void omxMarkDirty(omxMatrix *om) { om->isDirty = TRUE; }
433
434 unsigned short omxMatrixNeedsUpdate(omxMatrix *om) {
435         for(int i = 0; i < om->numPopulateLocations; i++) {
436                 int index = om->populateFrom[i];
437                 omxMatrix* sourceMatrix;
438                 if (index < 0) {
439                         sourceMatrix = om->currentState->matrixList[~index];
440                 } else {
441                         sourceMatrix = om->currentState->algebraList[index];
442                 }
443                 if(omxNeedsUpdate(sourceMatrix)) return TRUE;   // Make sure it's up to date
444         }
445     return FALSE;
446 };
447
448 omxMatrix* omxNewMatrixFromMxMatrix(SEXP mxMatrix, omxState* state) {
449 /* Creates and populates an omxMatrix with details from an R MxMatrix object. */
450         omxMatrix *om = NULL;
451         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
452         return omxFillMatrixFromMxMatrix(om, mxMatrix, state);
453 }
454
455 omxMatrix* omxNewMatrixFromRPrimitive(SEXP rObject, omxState* state, 
456         unsigned short hasMatrixNumber, int matrixNumber) {
457 /* Creates and populates an omxMatrix with details from an R matrix object. */
458         omxMatrix *om = NULL;
459         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
460         return omxFillMatrixFromRPrimitive(om, rObject, state, hasMatrixNumber, matrixNumber);
461 }
462
463 omxMatrix* omxFillMatrixFromMxMatrix(omxMatrix* om, SEXP mxMatrix, omxState* state) {
464 /* Populates the fields of a omxMatrix with details from an R Matrix. */
465         if(inherits(mxMatrix, "MxMatrix")) {
466                 if(OMX_DEBUG) { Rprintf("R matrix is Mx Matrix.  Processing.\n"); }
467                 SEXP matrix;
468                 PROTECT(matrix = GET_SLOT(mxMatrix,  install("values")));
469                 om = fillMatrixHelperFunction(om, matrix, state, 0, 0);
470                 UNPROTECT(1);
471         } else {
472                 error("Recieved unknown matrix type in omxFillMatrixFromMxMatrix.");
473         }
474         return(om);
475 }
476
477 omxMatrix* omxFillMatrixFromRPrimitive(omxMatrix* om, SEXP rObject, omxState* state,
478         unsigned short hasMatrixNumber, int matrixNumber) {
479 /* Populates the fields of a omxMatrix with details from an R object. */
480         if(!isMatrix(rObject) && !isVector(rObject)) { // Sanity Check
481                 error("Recieved unknown matrix type in omxFillMatrixFromRPrimitive.");
482         }
483         return(fillMatrixHelperFunction(om, rObject, state, hasMatrixNumber, matrixNumber));
484 }
485
486
487
488 omxMatrix* fillMatrixHelperFunction(omxMatrix* om, SEXP matrix, omxState* state,
489         unsigned short hasMatrixNumber, int matrixNumber) {
490
491         SEXP matrixDims;
492         int* dimList;
493
494         if(OMX_DEBUG) { Rprintf("Filling omxMatrix from R matrix.\n"); }
495
496         if(om == NULL) {
497                 om = omxInitMatrix(NULL, 0, 0, FALSE, state);
498         }
499
500         om->data = REAL(AS_NUMERIC(matrix));    // TODO: Class-check first?
501
502         if(isMatrix(matrix)) {
503                 PROTECT(matrixDims = getAttrib(matrix, R_DimSymbol));
504                 dimList = INTEGER(matrixDims);
505                 om->rows = dimList[0];
506                 om->cols = dimList[1];
507                 UNPROTECT(1);   // MatrixDims
508         } else if (isVector(matrix)) {          // If it's a vector, assume it's a row vector. BLAS doesn't care.
509                 if(OMX_DEBUG) { Rprintf("Vector discovered.  Assuming rowity.\n"); }
510                 om->rows = 1;
511                 om->cols = length(matrix);
512         }
513         if(OMX_DEBUG) { Rprintf("Matrix connected to (%d, %d) matrix or MxMatrix.\n", om->rows, om->cols); }
514
515         om->localData = FALSE;
516         om->colMajor = TRUE;
517         om->originalRows = om->rows;
518         om->originalCols = om->cols;
519         om->originalColMajor = TRUE;
520         om->aliasedPtr = NULL;
521         om->algebra = NULL;
522         om->objective = NULL;
523         om->currentState = state;
524         om->lastCompute = -1;
525         om->lastRow = -1;
526         om->hasMatrixNumber = hasMatrixNumber;
527         om->matrixNumber = matrixNumber;
528
529
530         if(OMX_DEBUG) { Rprintf("Pre-compute call.\n");}
531         omxMatrixCompute(om);
532         if(OMX_DEBUG) { Rprintf("Post-compute call.\n");}
533
534         if(OMX_DEBUG) {
535                 omxPrintMatrix(om, "Finished importing matrix");
536         }
537
538         return om;
539 }
540
541 void omxProcessMatrixPopulationList(omxMatrix* matrix, SEXP matStruct) {
542
543         if(OMX_DEBUG) { Rprintf("Processing Population List: %d elements.\n", length(matStruct) - 1); }
544         SEXP subList;
545
546         if(length(matStruct) > 1) {
547                 int numPopLocs = length(matStruct) - 1;
548                 matrix->numPopulateLocations = numPopLocs;
549                 matrix->populateFrom = (int*)R_alloc(numPopLocs, sizeof(int));
550                 matrix->populateFromRow = (int*)R_alloc(numPopLocs, sizeof(int));
551                 matrix->populateFromCol = (int*)R_alloc(numPopLocs, sizeof(int));
552                 matrix->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
553                 matrix->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
554         }
555
556         for(int i = 0; i < length(matStruct)-1; i++) {
557                 PROTECT(subList = AS_INTEGER(VECTOR_ELT(matStruct, i+1)));
558
559                 int* locations = INTEGER(subList);
560                 if(OMX_DEBUG) { Rprintf("."); } //:::
561                 matrix->populateFrom[i] = locations[0];
562                 matrix->populateFromRow[i] = locations[1];
563                 matrix->populateFromCol[i] = locations[2];
564                 matrix->populateToRow[i] = locations[3];
565                 matrix->populateToCol[i] = locations[4];
566
567                 UNPROTECT(1); // subList
568         }
569 }
570
571 void omxToggleRowColumnMajor(omxMatrix *mat) {
572
573         int i, j;
574         int nrows = mat->rows;
575         int ncols = mat->cols;
576         
577         double *newdata = (double*) Calloc(nrows * ncols, double);
578         double *olddata = mat->data;
579
580         if (mat->colMajor) {
581                 for(i = 0; i < ncols; i++) {
582                         for(j = 0; j < nrows; j++) {
583                                 newdata[i + ncols * j] = olddata[i * nrows + j];
584                         }
585                 }
586         } else {
587                 for(i = 0; i < nrows; i++) {
588                         for(j = 0; j < ncols; j++) {
589                                 newdata[i + nrows * j] = olddata[i * ncols + j];
590                         }
591                 }
592         }
593
594         if (mat->localData) {
595                 Free(mat->data);
596         }
597
598         mat->localData = TRUE;
599         mat->data = newdata;
600         mat->colMajor = !mat->colMajor;
601 }
602
603 void omxTransposeMatrix(omxMatrix *mat) {
604         mat->colMajor = !mat->colMajor;
605         
606         if(mat->rows != mat->cols){
607         int mid = mat->rows;
608         mat->rows = mat->cols;
609         mat->cols = mid;
610         }
611         
612         omxMatrixCompute(mat);
613 }
614
615 void omxRemoveRowsAndColumns(omxMatrix *om, int numRowsRemoved, int numColsRemoved, int rowsRemoved[], int colsRemoved[])
616 {
617     // TODO: Create short-circuit form of omxRemoveRowsAndCols to remove just rows or just columns.
618 //      if(OMX_DEBUG_MATRIX) { Rprintf("Removing %d rows and %d columns from 0x%0x.\n", numRowsRemoved, numColsRemoved, om);}
619
620         if(numRowsRemoved < 1 && numColsRemoved < 1) { return; }
621
622         int oldRows, oldCols;
623
624         if(om->aliasedPtr == NULL) {
625                 if(om->originalRows == 0 || om->originalCols == 0) {
626                         om->originalRows = om->rows;
627                         om->originalCols = om->cols;
628                 }
629                 oldRows = om->originalRows;
630                 oldCols = om->originalCols;
631         } else {
632                 oldRows = om->aliasedPtr->rows;
633                 oldCols = om->aliasedPtr->cols;
634         }
635
636         int nextCol = 0;
637         int nextRow = 0;
638
639         if(om->rows > om->originalRows || om->cols > om->originalCols) {        // sanity check.
640                 error("Aliased Matrix is too small for alias.");
641         }
642
643         om->rows = oldRows - numRowsRemoved;
644         om->cols = oldCols - numColsRemoved;
645
646         // Note:  This really aught to be done using a matrix multiply.  Why isn't it?
647         for(int j = 0; j < oldCols; j++) {
648                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Handling column %d/%d...", j, oldCols);}
649                 if(colsRemoved[j]) {
650                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Removed.\n");}
651                         continue;
652                 } else {
653                         nextRow = 0;
654                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Rows (max %d): ", oldRows); }
655                         for(int k = 0; k < oldRows; k++) {
656                                 if(rowsRemoved[k]) {
657                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d removed....", k);}
658                                         continue;
659                                 } else {
660                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d kept....", k);}
661                                         if(om->aliasedPtr == NULL) {
662                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Self-aliased matrix access.\n");}
663                                                 omxSetMatrixElement(om, nextRow, nextCol, omxAliasedMatrixElement(om, k, j));
664                                         } else {
665                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Matrix 0x%x re-aliasing to 0x%x.\n", om, om->aliasedPtr);}
666                                                 omxSetMatrixElement(om, nextRow, nextCol, omxMatrixElement(om->aliasedPtr, k,  j));
667                                         }
668                                         nextRow++;
669                                 }
670                         }
671                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("\n");}
672                         nextCol++;
673                 }
674         }
675
676         omxMatrixCompute(om);
677 }
678
679 /* Function wrappers that switch based on inclusion of algebras */
680 void omxPrint(omxMatrix *source, char* d) {                                     // Pretty-print a (small) matrix
681         if(source->algebra != NULL) omxAlgebraPrint(source->algebra, d);
682         else if(source->objective != NULL) omxObjectivePrint(source->objective, d);
683         else omxPrintMatrix(source, d);
684 }
685
686 unsigned short omxNeedsUpdate(omxMatrix *matrix) {
687         unsigned short retval;
688         /* Simplest update check: If we're dirty or haven't computed this cycle (iteration or row), we need to. */
689         // TODO : Implement a dependency-tree-based dirtiness propagation system
690         if(OMX_DEBUG_MATRIX) {Rprintf("Matrix 0x%x NeedsUpdate?", matrix);}
691
692         if(matrix == NULL) {
693                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix argument is NULL. ");}
694                 retval = FALSE;         // Not existing means never having to say you need to recompute.
695         } else if(matrix->isDirty) {
696                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix is dirty. ");}
697                 retval = TRUE;
698         } else if(matrix->lastCompute < matrix->currentState->computeCount) {
699                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix last compute is less than current compute count. ");}
700                 retval = TRUE;          // No need to check args if oa's dirty.
701         } else if(matrix->lastRow != matrix->currentState->currentRow) {
702                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix last row is less than current row. ");}
703                 retval = TRUE;                  // Ditto.
704         } else if(matrix->algebra != NULL) {
705                 if(OMX_DEBUG_MATRIX) {Rprintf("checking algebra needs update. ");}
706                 retval = omxAlgebraNeedsUpdate(matrix->algebra);
707         } else if(matrix->objective != NULL) {
708                 if(OMX_DEBUG_MATRIX) {Rprintf("checking objective function needs update. ");}
709                 retval = omxObjectiveNeedsUpdate(matrix->objective);
710         } else {
711                 if(OMX_DEBUG_MATRIX) {Rprintf("checking matrix needs update. ");}
712                 retval = omxMatrixNeedsUpdate(matrix);
713         }
714         if(OMX_DEBUG_MATRIX && retval) {Rprintf("Yes.\n");}
715         if(OMX_DEBUG_MATRIX && !retval) {Rprintf("No.\n");}
716         return(retval);
717 }
718
719 void omxRecompute(omxMatrix *matrix) {
720         if(!omxNeedsUpdate(matrix)) return;
721         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
722         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
723         else omxMatrixCompute(matrix);
724 }
725
726 void omxCompute(omxMatrix *matrix) {
727         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
728         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
729         else omxMatrixCompute(matrix);
730 }