Revert "Store algebraList in std::vector"
[openmx:openmx.git] / src / omxRowFitFunction.cpp
1 /*
2  *  Copyright 2007-2013 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 #include <R.h>
18 #include <Rinternals.h>
19 #include <Rdefines.h>
20 #include <R_ext/Rdynload.h>
21 #include <R_ext/BLAS.h>
22 #include <R_ext/Lapack.h>
23 #include "omxDefines.h"
24 #include "omxAlgebraFunctions.h"
25 #include "omxSymbolTable.h"
26 #include "omxData.h"
27 #include "omxRowFitFunction.h"
28 #include "omxFIMLFitFunction.h"
29
30 void omxDestroyRowFitFunction(omxFitFunction *oo) {
31
32         omxRowFitFunction* argStruct = (omxRowFitFunction*)(oo->argStruct);
33
34         omxFreeMatrixData(argStruct->dataRow);
35 }
36
37 omxRListElement* omxSetFinalReturnsRowFitFunction(omxFitFunction *oo, int *numReturns) {
38         *numReturns = 0;
39         omxRListElement* retVal = (omxRListElement*) R_alloc(1, sizeof(omxRListElement));
40
41         retVal[0].numValues = 0;
42
43         return retVal;
44 }
45
46
47 void omxCopyMatrixToRow(omxMatrix* source, int row, omxMatrix* target) {
48         
49         int i;
50         for(i = 0; i < source->cols; i++) {
51                 omxSetMatrixElement(target, row, i, omxMatrixElement(source, 0, i));
52         }
53
54 }
55
56 void markDataRowDependencies(omxState* os, omxRowFitFunction* orff) {
57
58         int numDeps = orff->numDataRowDeps;
59         int *deps = orff->dataRowDeps;
60
61         omxMatrix** algebraList = os->algebraList;
62
63         for (int i = 0; i < numDeps; i++) {
64                 int value = deps[i];
65
66                 if(value < 0) {
67                         omxMarkDirty(os->matrixList[~value]);
68                 } else {
69                         omxMarkDirty(algebraList[value]);
70                 }
71         }
72
73 }
74
75 void omxRowFitFunctionSingleIteration(omxFitFunction *localobj, omxFitFunction *sharedobj, int rowbegin, int rowcount) {
76
77     omxRowFitFunction* oro = ((omxRowFitFunction*) localobj->argStruct);
78     omxRowFitFunction* shared_oro = ((omxRowFitFunction*) sharedobj->argStruct);
79
80         int numDefs;
81
82     omxMatrix *rowAlgebra, *rowResults;
83     omxMatrix *filteredDataRow, *dataRow, *existenceVector;
84     omxMatrix *dataColumns;
85         omxDefinitionVar* defVars;
86         omxData *data;
87         int isContiguous, contiguousStart, contiguousLength;
88     double* oldDefs;
89     int numCols, numRemoves;
90
91         rowAlgebra          = oro->rowAlgebra;
92         rowResults          = shared_oro->rowResults;
93         data                = oro->data;
94         defVars             = oro->defVars;
95         numDefs             = oro->numDefs;
96     oldDefs         = oro->oldDefs;
97     dataColumns     = oro->dataColumns;
98     dataRow         = oro->dataRow;
99     filteredDataRow = oro->filteredDataRow;
100     existenceVector = oro->existenceVector;
101     
102     isContiguous    = oro->contiguous.isContiguous;
103         contiguousStart = oro->contiguous.start;
104         contiguousLength = oro->contiguous.length;
105
106         numCols = dataColumns->cols;
107         int *toRemove = (int*) malloc(sizeof(int) * dataColumns->cols);
108         int *zeros = (int*) calloc(dataColumns->cols, sizeof(int));
109
110     resetDefinitionVariables(oldDefs, numDefs);
111
112         for(int row = rowbegin; row < data->rows && (row - rowbegin) < rowcount; row++) {
113
114                 // Handle Definition Variables.
115         if(OMX_DEBUG_ROWS(row)) { Rprintf("numDefs is %d", numDefs);}
116                 if(numDefs != 0) {              // With defs, just copy repeatedly to the rowResults matrix.
117                         handleDefinitionVarList(data, localobj->matrix->currentState, row, defVars, oldDefs, numDefs);
118                 }
119
120                 omxStateNextRow(localobj->matrix->currentState);                                                // Advance row
121                 
122         // Populate data row
123                 numRemoves = 0;
124         
125                 if (isContiguous) {
126                         omxContiguousDataRow(data, row, contiguousStart, contiguousLength, dataRow);
127                 } else {
128                         omxDataRow(data, row, dataColumns, dataRow);    // Populate data row
129                 }
130
131                 markDataRowDependencies(localobj->matrix->currentState, oro);
132                 
133                 for(int j = 0; j < dataColumns->cols; j++) {
134                         double dataValue = omxVectorElement(dataRow, j);
135                         if(isnan(dataValue)) {
136                                 numRemoves++;
137                                 toRemove[j] = 1;
138                 omxSetVectorElement(existenceVector, j, 0);
139                         } else {
140                             toRemove[j] = 0;
141                 omxSetVectorElement(existenceVector, j, 1);
142                         }
143                 }               
144                 // TODO: Determine if this is the correct response.
145                 
146                 if(numRemoves == numCols) {
147                         char *errstr = (char*) calloc(250, sizeof(char));
148                         sprintf(errstr, "Row %d completely missing.  omxRowFitFunction cannot have completely missing rows.", omxDataIndex(data, row));
149                         omxRaiseError(localobj->matrix->currentState, -1, errstr);
150                         free(errstr);
151                         continue;
152                 }
153
154                 omxResetAliasedMatrix(filteredDataRow);                         // Reset the row
155                 omxRemoveRowsAndColumns(filteredDataRow, 0, numRemoves, zeros, toRemove);
156
157                 omxRecompute(rowAlgebra);                                                       // Compute this row
158
159                 omxCopyMatrixToRow(rowAlgebra, omxDataIndex(data, row), rowResults);
160         }
161         free(toRemove);
162         free(zeros);
163 }
164
165 static void omxCallRowFitFunction(omxFitFunction *oo, int want, double *gradient) {     // TODO: Figure out how to give access to other per-iteration structures.
166     if(OMX_DEBUG) { Rprintf("Beginning Row Evaluation.\n");}
167         // Requires: Data, means, covariances.
168
169         omxMatrix* objMatrix  = oo->matrix;
170         omxState* parentState = objMatrix->currentState;
171         int numChildren = parentState->numChildren;
172
173     omxMatrix *reduceAlgebra;
174         omxData *data;
175
176     omxRowFitFunction* oro = ((omxRowFitFunction*) oo->argStruct);
177
178         reduceAlgebra   = oro->reduceAlgebra;
179         data                = oro->data;
180
181         /* Michael Spiegel, 7/31/12
182         * The demo "RowFitFunctionSimpleExamples" will fail in the parallel 
183         * Hessian calculation if the resizing operation is performed.
184         *
185         omxMatrix *rowAlgebra, *rowResults
186         rowAlgebra          = oro->rowAlgebra;
187         rowResults          = oro->rowResults;
188
189         if(rowResults->cols != rowAlgebra->cols || rowResults->rows != data->rows) {
190                 if(OMX_DEBUG_ROWS(1)) { 
191                         Rprintf("Resizing rowResults from %dx%d to %dx%d.\n", 
192                                 rowResults->rows, rowResults->cols, 
193                                 data->rows, rowAlgebra->cols); 
194                 }
195                 omxResizeMatrix(rowResults, data->rows, rowAlgebra->cols, FALSE);
196         }
197         */
198                 
199     int parallelism = (numChildren == 0) ? 1 : numChildren;
200
201         if (parallelism > data->rows) {
202                 parallelism = data->rows;
203         }
204
205         if (parallelism > 1) {
206                 int stride = (data->rows / parallelism);
207
208                 #pragma omp parallel for num_threads(parallelism) 
209                 for(int i = 0; i < parallelism; i++) {
210                         omxMatrix *childMatrix = omxLookupDuplicateElement(parentState->childList[i], objMatrix);
211                         omxFitFunction *childFit = childMatrix->fitFunction;
212                         if (i == parallelism - 1) {
213                                 omxRowFitFunctionSingleIteration(childFit, oo, stride * i, data->rows - stride * i);
214                         } else {
215                                 omxRowFitFunctionSingleIteration(childFit, oo, stride * i, stride);
216                         }
217                 }
218
219                 for(int i = 0; i < parallelism; i++) {
220                         if (parentState->childList[i]->statusCode < 0) {
221                                 parentState->statusCode = parentState->childList[i]->statusCode;
222                                 strncpy(parentState->statusMsg, parentState->childList[i]->statusMsg, 249);
223                                 parentState->statusMsg[249] = '\0';
224                         }
225                 }
226
227         } else {
228                 omxRowFitFunctionSingleIteration(oo, oo, 0, data->rows);
229         }
230
231         omxRecompute(reduceAlgebra);
232
233         omxCopyMatrix(oo->matrix, reduceAlgebra);
234
235 }
236
237 void omxInitRowFitFunction(omxFitFunction* oo, SEXP rObj) {
238
239         if(OMX_DEBUG) { Rprintf("Initializing Row/Reduce fit function.\n"); }
240
241         SEXP nextMatrix, itemList, nextItem;
242         int nextDef, index, numDeps;
243
244         omxRowFitFunction *newObj = (omxRowFitFunction*) R_alloc(1, sizeof(omxRowFitFunction));
245
246         if(OMX_DEBUG) {Rprintf("Accessing data source.\n"); }
247         PROTECT(nextMatrix = GET_SLOT(rObj, install("data")));
248         newObj->data = omxDataLookupFromState(nextMatrix, oo->matrix->currentState);
249         if(newObj->data == NULL) {
250                 char *errstr = (char*) calloc(250, sizeof(char));
251                 sprintf(errstr, "No data provided to omxRowFitFunction.");
252                 omxRaiseError(oo->matrix->currentState, -1, errstr);
253                 free(errstr);
254         }
255         UNPROTECT(1); // nextMatrix
256
257         PROTECT(nextMatrix = GET_SLOT(rObj, install("rowAlgebra")));
258         newObj->rowAlgebra = omxMatrixLookupFromState1(nextMatrix, oo->matrix->currentState);
259         if(newObj->rowAlgebra == NULL) {
260                 char *errstr = (char*) calloc(250, sizeof(char));
261                 sprintf(errstr, "No row-wise algebra in omxRowFitFunction.");
262                 omxRaiseError(oo->matrix->currentState, -1, errstr);
263                 free(errstr);
264         }
265         UNPROTECT(1);// nextMatrix
266
267         PROTECT(nextMatrix = GET_SLOT(rObj, install("filteredDataRow")));
268         newObj->filteredDataRow = omxMatrixLookupFromState1(nextMatrix, oo->matrix->currentState);
269         if(newObj->filteredDataRow == NULL) {
270                 char *errstr = (char*) calloc(250, sizeof(char));
271                 sprintf(errstr, "No row results matrix in omxRowFitFunction.");
272                 omxRaiseError(oo->matrix->currentState, -1, errstr);
273                 free(errstr);
274         }
275         // Create the original data row from which to filter.
276     newObj->dataRow = omxInitMatrix(NULL, newObj->filteredDataRow->rows, newObj->filteredDataRow->cols, TRUE, oo->matrix->currentState);
277     omxAliasMatrix(newObj->filteredDataRow, newObj->dataRow);
278         UNPROTECT(1);// nextMatrix
279
280         PROTECT(nextMatrix = GET_SLOT(rObj, install("existenceVector")));
281         newObj->existenceVector = omxMatrixLookupFromState1(nextMatrix, oo->matrix->currentState);
282     // Do we allow NULL existence?  (Whoa, man. That's, like, deep, or something.)
283         if(newObj->existenceVector == NULL) {
284                 char *errstr = (char*) calloc(250, sizeof(char));
285                 sprintf(errstr, "No existance matrix in omxRowFitFunction.");
286                 omxRaiseError(oo->matrix->currentState, -1, errstr);
287                 free(errstr);
288         }
289         UNPROTECT(1);// nextMatrix
290
291
292         PROTECT(nextMatrix = GET_SLOT(rObj, install("rowResults")));
293         newObj->rowResults = omxMatrixLookupFromState1(nextMatrix, oo->matrix->currentState);
294         if(newObj->rowResults == NULL) {
295                 char *errstr = (char*) calloc(250, sizeof(char));
296                 sprintf(errstr, "No row results matrix in omxRowFitFunction.");
297                 omxRaiseError(oo->matrix->currentState, -1, errstr);
298                 free(errstr);
299         }
300         UNPROTECT(1);// nextMatrix
301
302         PROTECT(nextMatrix = GET_SLOT(rObj, install("reduceAlgebra")));
303         newObj->reduceAlgebra = omxMatrixLookupFromState1(nextMatrix, oo->matrix->currentState);
304         if(newObj->reduceAlgebra == NULL) {
305                 char *errstr = (char*) calloc(250, sizeof(char));
306                 sprintf(errstr, "No row reduction algebra in omxRowFitFunction.");
307                 omxRaiseError(oo->matrix->currentState, -1, errstr);
308                 free(errstr);
309         }
310         UNPROTECT(1);// nextMatrix
311         
312         if(OMX_DEBUG) {Rprintf("Accessing variable mapping structure.\n"); }
313         PROTECT(nextMatrix = GET_SLOT(rObj, install("dataColumns")));
314         newObj->dataColumns = omxNewMatrixFromRPrimitive(nextMatrix, oo->matrix->currentState, 0, 0);
315         if(OMX_DEBUG) { omxPrint(newObj->dataColumns, "Variable mapping"); }
316         UNPROTECT(1);
317
318         if(OMX_DEBUG) {Rprintf("Accessing data row dependencies.\n"); }
319         PROTECT(nextItem = GET_SLOT(rObj, install("dataRowDeps")));
320         numDeps = LENGTH(nextItem);
321         newObj->numDataRowDeps = numDeps;
322         newObj->dataRowDeps = (int*) R_alloc(numDeps, sizeof(int));
323         for(int i = 0; i < numDeps; i++) {
324                 newObj->dataRowDeps[i] = INTEGER(nextItem)[i];
325         }
326         UNPROTECT(1);
327
328         if(OMX_DEBUG) {Rprintf("Accessing definition variables structure.\n"); }
329         PROTECT(nextMatrix = GET_SLOT(rObj, install("definitionVars")));
330         newObj->numDefs = length(nextMatrix);
331         newObj->oldDefs = (double *) R_alloc(newObj->numDefs, sizeof(double));          // Storage for Def Vars
332         if(OMX_DEBUG) {Rprintf("Number of definition variables is %d.\n", newObj->numDefs); }
333         newObj->defVars = (omxDefinitionVar *) R_alloc(newObj->numDefs, sizeof(omxDefinitionVar));
334         for(nextDef = 0; nextDef < newObj->numDefs; nextDef++) {
335                 SEXP dataSource, columnSource, depsSource; 
336
337                 PROTECT(itemList = VECTOR_ELT(nextMatrix, nextDef));
338                 PROTECT(dataSource = VECTOR_ELT(itemList, 0));
339                 if(OMX_DEBUG) {Rprintf("Data source number is %d.\n", INTEGER(dataSource)[0]); }
340                 newObj->defVars[nextDef].data = INTEGER(dataSource)[0];
341                 newObj->defVars[nextDef].source = oo->matrix->currentState->dataList[INTEGER(dataSource)[0]];
342                 PROTECT(columnSource = VECTOR_ELT(itemList, 1));
343                 if(OMX_DEBUG) {Rprintf("Data column number is %d.\n", INTEGER(columnSource)[0]); }
344                 newObj->defVars[nextDef].column = INTEGER(columnSource)[0];
345                 PROTECT(depsSource = VECTOR_ELT(itemList, 2));
346                 numDeps = LENGTH(depsSource);
347                 newObj->defVars[nextDef].numDeps = numDeps;
348                 newObj->defVars[nextDef].deps = (int*) R_alloc(numDeps, sizeof(int));
349                 for(int i = 0; i < numDeps; i++) {
350                         newObj->defVars[nextDef].deps[i] = INTEGER(depsSource)[i];
351                 }
352                 UNPROTECT(3); // unprotect dataSource, columnSource, and depsSource
353
354                 newObj->defVars[nextDef].numLocations = length(itemList) - 3;
355                 newObj->defVars[nextDef].matrices = (int *) R_alloc(length(itemList) - 3, sizeof(int));
356                 newObj->defVars[nextDef].rows = (int *) R_alloc(length(itemList) - 3, sizeof(int));
357                 newObj->defVars[nextDef].cols = (int *) R_alloc(length(itemList) - 3, sizeof(int));
358
359                 for(index = 3; index < length(itemList); index++) {
360                         PROTECT(nextItem = VECTOR_ELT(itemList, index));
361                         newObj->defVars[nextDef].matrices[index-3] = INTEGER(nextItem)[0];
362                         newObj->defVars[nextDef].rows[index-3]     = INTEGER(nextItem)[1];
363                         newObj->defVars[nextDef].cols[index-3]     = INTEGER(nextItem)[2];
364                         UNPROTECT(1); // unprotect nextItem
365                 }
366                 UNPROTECT(1); // unprotect itemList
367         }
368         UNPROTECT(1); // unprotect nextMatrix
369         
370         /* Set up data columns */
371         omxSetContiguousDataColumns(&(newObj->contiguous), newObj->data, newObj->dataColumns);
372
373         oo->computeFun = omxCallRowFitFunction;
374         oo->setFinalReturns = omxSetFinalReturnsRowFitFunction;
375         oo->destructFun = omxDestroyRowFitFunction;
376         oo->repopulateFun = NULL;
377         oo->usesChildModels = TRUE;
378
379         oo->argStruct = (void*) newObj;
380 }
381
382