src/omxWLSFitFunction.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 "omxAlgebraFunctions.h"
  24 #include "omxWLSFitFunction.h"
  25
  26 void flattenDataToVector(omxMatrix* cov, omxMatrix* means, omxThresholdColumn* thresholds, int nThresholds, omxMatrix* vector) {
  27     // TODO: vectorize data flattening
  28     int nextLoc = 0;
  29     for(int j = 0; j < cov->rows; j++) {
  30         for(int k = 0; k <= j; k++) {
  31             omxSetVectorElement(vector, nextLoc, omxMatrixElement(cov, k, j)); // Use upper triangle in case of SYMM-style mat.
  32             nextLoc++;
  33         }
  34     }
  35     if (means != NULL) {
  36         for(int j = 0; j < cov->rows; j++) {
  37             omxSetVectorElement(vector, nextLoc, omxVectorElement(means, j));
  38             nextLoc++;
  39         }
  40     }
  41     if (thresholds != NULL) {
  42         for(int j = 0; j < nThresholds; j++) {
  43             omxThresholdColumn* thresh = thresholds + j;
  44             for(int k = 0; k < thresh->numThresholds; k++) {
  45                 omxSetVectorElement(vector, nextLoc, omxMatrixElement(thresh->matrix, k, thresh->column));
  46                 nextLoc++;
  47             }
  48         }
  49     }
  50 }
  51
  52 void omxDestroyWLSFitFunction(omxFitFunction *oo) {
  53
  54         if(OMX_DEBUG) {mxLog("Freeing WLS FitFunction.");}
  55         omxWLSFitFunction* owo = ((omxWLSFitFunction*)oo->argStruct);
  56
  57     if(owo->observedFlattened != NULL) omxFreeMatrixData(owo->observedFlattened);
  58     if(owo->expectedFlattened != NULL) omxFreeMatrixData(owo->expectedFlattened);
  59     if(owo->weights != NULL) omxFreeMatrixData(owo->weights);
  60     if(owo->B != NULL) omxFreeMatrixData(owo->B);
  61     if(owo->P != NULL) omxFreeMatrixData(owo->P);
  62 }
  63
  64 omxRListElement* omxSetFinalReturnsWLSFitFunction(omxFitFunction *oo, int *numReturns) {
  65         *numReturns = 4;
  66         omxRListElement* retVal = (omxRListElement*) R_alloc(*numReturns, sizeof(omxRListElement));
  67
  68         retVal[0].numValues = 1;
  69         retVal[0].values = (double*) R_alloc(1, sizeof(double));
  70         strncpy(retVal[0].label, "Minus2LogLikelihood", 20);
  71         retVal[0].values[0] = omxMatrixElement(oo->matrix, 0, 0);
  72
  73         retVal[1].numValues = 1;
  74         retVal[1].values = (double*) R_alloc(1, sizeof(double));
  75         strncpy(retVal[1].label, "SaturatedLikelihood", 20);
  76     retVal[1].values[0] = 0;
  77
  78         retVal[2].numValues = 1;
  79         retVal[2].values = (double*) R_alloc(1, sizeof(double));
  80         strncpy(retVal[2].label, "IndependenceLikelihood", 23);
  81     retVal[2].values[0] = 0;
  82
  83         retVal[3].numValues = 1;
  84         retVal[3].values = (double*) R_alloc(1, sizeof(double));
  85         strncpy(retVal[3].label, "ADFMisfit", 20);
  86         retVal[3].values[0] = omxMatrixElement(oo->matrix, 0, 0);
  87
  88     // retVal[1].numValues = 1;
  89     // retVal[1].values = (double*) R_alloc(1, sizeof(double));
  90     // strncpy(retVal[1].label, "SaturatedADFMisfit", 20);
  91     // retVal[1].values[0] = (((omxWLSFitFunction*)oo->argStruct)->logDetObserved + ncols) * (((omxWLSFitFunction*)oo->argStruct)->n - 1);
  92     //
  93     // retVal[2].numValues = 1;
  94     // retVal[2].values = (double*) R_alloc(1, sizeof(double));
  95     // strncpy(retVal[2].label, "IndependenceLikelihood", 23);
  96     // // Independence model assumes all-zero manifest covariances.
  97     // // (det(expected) + tr(observed * expected^-1)) * (n - 1);
  98     // // expected is the diagonal of the observed.  Inverse expected is 1/each diagonal value.
  99     // // Therefore the diagonal elements of observed * expected^-1 are each 1.
 100     // // So the trace of the matrix is the same as the number of columns.
 101     // // The determinant of a diagonal matrix is the product of the diagonal elements.
 102     // // Since these are the same in the expected as in the observed, we can get 'em direct.
 103     // for(int i = 0; i < ncols; i++) {
 104     //  // We sum logs instead of logging the product.
 105     //  det += log(omxMatrixElement(cov, i, i));
 106     // }
 107     // if(OMX_DEBUG) { mxLog("det: %f, tr: %f, n= %d, total:%f", det, ncols, ((omxWLSFitFunction*)oo->argStruct)->n, (ncols + det) * (((omxWLSFitFunction*)oo->argStruct)->n - 1)); }
 108     // if(OMX_DEBUG) { omxPrint(cov, "Observed:"); }
 109     // retVal[2].values[0] = (ncols + det) * (((omxWLSFitFunction*)oo->argStruct)->n - 1);
 110
 111         return retVal;
 112 }
 113
 114 static void omxCallWLSFitFunction(omxFitFunction *oo, int want, double *gradient, double *hessian) {    // TODO: Figure out how to give access to other per-iteration structures.
 115
 116         if(OMX_DEBUG) { mxLog("Beginning WLS Evaluation.");}
 117         // Requires: Data, means, covariances.
 118
 119         double sum = 0.0;
 120
 121         omxMatrix *oCov, *oMeans, *eCov, *eMeans, *P, *B, *weights, *oFlat, *eFlat;
 122
 123     omxThresholdColumn *oThresh, *eThresh;
 124
 125         omxWLSFitFunction *owo = ((omxWLSFitFunction*)oo->argStruct);
 126
 127     /* Locals for readability.  Compiler should cut through this. */
 128         oCov            = owo->observedCov;
 129         oMeans          = owo->observedMeans;
 130         oThresh         = owo->observedThresholds;
 131         eCov            = owo->expectedCov;
 132         eMeans          = owo->expectedMeans;
 133         eThresh         = owo->expectedThresholds;
 134         oFlat           = owo->observedFlattened;
 135         eFlat           = owo->expectedFlattened;
 136         weights         = owo->weights;
 137         B                       = owo->B;
 138         P                       = owo->P;
 139     int nThresh = owo->nThresholds;
 140     int onei    = 1;
 141
 142         omxExpectation* expectation = oo->expectation;
 143
 144     /* Recompute and recopy */
 145         omxExpectationCompute(expectation);
 146
 147         flattenDataToVector(oCov, oMeans, oThresh, nThresh, oFlat);
 148         flattenDataToVector(eCov, eMeans, eThresh, nThresh, eFlat);
 149
 150         omxCopyMatrix(B, oFlat);
 151
 152         omxDAXPY(-1.0, eFlat, B);
 153
 154     if(weights != NULL) {
 155         omxDGEMV(TRUE, 1.0, weights, B, 0.0, P);
 156     } else {
 157         // ULS Case: Memcpy faster than dgemv.
 158         // TODO: Better to use an omxMatrix duplicator here.
 159         memcpy(P, B, B->cols*sizeof(double));
 160     }
 161     sum = F77_CALL(ddot)(&(P->cols), P->data, &onei, B->data, &onei);
 162
 163     oo->matrix->data[0] = sum;
 164
 165         if(OMX_DEBUG) { mxLog("WLSFitFunction value comes to: %f.", oo->matrix->data[0]); }
 166
 167 }
 168
 169 void omxPopulateWLSAttributes(omxFitFunction *oo, SEXP algebra) {
 170     if(OMX_DEBUG) { mxLog("Populating WLS Attributes."); }
 171
 172         omxWLSFitFunction *argStruct = ((omxWLSFitFunction*)oo->argStruct);
 173         omxMatrix *expCovInt = argStruct->expectedCov;                  // Expected covariance
 174         omxMatrix *expMeanInt = argStruct->expectedMeans;                       // Expected means
 175         omxMatrix *weightInt = argStruct->weights;                      // Expected means
 176
 177         SEXP expCovExt, expMeanExt, weightExt, gradients;
 178         PROTECT(expCovExt = allocMatrix(REALSXP, expCovInt->rows, expCovInt->cols));
 179         for(int row = 0; row < expCovInt->rows; row++)
 180                 for(int col = 0; col < expCovInt->cols; col++)
 181                         REAL(expCovExt)[col * expCovInt->rows + row] =
 182                                 omxMatrixElement(expCovInt, row, col);
 183
 184         if (expMeanInt != NULL) {
 185                 PROTECT(expMeanExt = allocMatrix(REALSXP, expMeanInt->rows, expMeanInt->cols));
 186                 for(int row = 0; row < expMeanInt->rows; row++)
 187                         for(int col = 0; col < expMeanInt->cols; col++)
 188                                 REAL(expMeanExt)[col * expMeanInt->rows + row] =
 189                                         omxMatrixElement(expMeanInt, row, col);
 190         } else {
 191                 PROTECT(expMeanExt = allocMatrix(REALSXP, 0, 0));
 192         }
 193
 194         PROTECT(weightExt = allocMatrix(REALSXP, weightInt->rows, weightInt->cols));
 195         for(int row = 0; row < weightInt->rows; row++)
 196                 for(int col = 0; col < weightInt->cols; col++)
 197                         REAL(weightExt)[col * weightInt->rows + row] =
 198                                 omxMatrixElement(weightInt, row, col);
 199
 200
 201         if(0) {  /* TODO fix for new internal API
 202                 int nLocs = Global->numFreeParams;
 203                 double gradient[Global->numFreeParams];
 204                 for(int loc = 0; loc < nLocs; loc++) {
 205                         gradient[loc] = NA_REAL;
 206                 }
 207                 //oo->gradientFun(oo, gradient);
 208                 PROTECT(gradients = allocMatrix(REALSXP, 1, nLocs));
 209
 210                 for(int loc = 0; loc < nLocs; loc++)
 211                         REAL(gradients)[loc] = gradient[loc];
 212                  */
 213         } else {
 214                 PROTECT(gradients = allocMatrix(REALSXP, 0, 0));
 215         }
 216
 217         setAttrib(algebra, install("expCov"), expCovExt);
 218         setAttrib(algebra, install("expMean"), expMeanExt);
 219         setAttrib(algebra, install("weights"), weightExt);
 220         setAttrib(algebra, install("gradients"), gradients);
 221
 222         UNPROTECT(4);
 223
 224 }
 225
 226 void omxSetWLSFitFunctionCalls(omxFitFunction* oo) {
 227
 228         /* Set FitFunction Calls to WLS FitFunction Calls */
 229         oo->computeFun = omxCallWLSFitFunction;
 230         oo->destructFun = omxDestroyWLSFitFunction;
 231         oo->setFinalReturns = omxSetFinalReturnsWLSFitFunction;
 232         oo->populateAttrFun = omxPopulateWLSAttributes;
 233 }
 234
 235 void omxInitWLSFitFunction(omxFitFunction* oo) {
 236
 237         omxMatrix *cov, *means, *weights;
 238
 239     if(OMX_DEBUG) { mxLog("Initializing WLS FitFunction function."); }
 240
 241     int vectorSize = 0;
 242
 243         omxSetWLSFitFunctionCalls(oo);
 244
 245         if(OMX_DEBUG) { mxLog("Retrieving expectation.\n"); }
 246         if (!oo->expectation) { error("%s requires an expectation", oo->fitType); }
 247
 248         if(OMX_DEBUG) { mxLog("Retrieving data.\n"); }
 249     omxData* dataMat = oo->expectation->data;
 250
 251         if(strncmp(omxDataType(dataMat), "acov", 4) != 0 && strncmp(omxDataType(dataMat), "cov", 3) != 0) {
 252                 char *errstr = (char*) calloc(250, sizeof(char));
 253                 sprintf(errstr, "WLS FitFunction unable to handle data type %s.  Data must be of type 'acov'.\n", omxDataType(dataMat));
 254                 omxRaiseError(oo->matrix->currentState, -1, errstr);
 255                 free(errstr);
 256                 if(OMX_DEBUG) { mxLog("WLS FitFunction unable to handle data type %s.  Aborting.", omxDataType(dataMat)); }
 257                 return;
 258         }
 259
 260         omxWLSFitFunction *newObj = (omxWLSFitFunction*) R_alloc(1, sizeof(omxWLSFitFunction));
 261
 262     /* Get Expectation Elements */
 263         newObj->expectedCov = omxGetExpectationComponent(oo->expectation, oo, "cov");
 264         newObj->expectedMeans = omxGetExpectationComponent(oo->expectation, oo, "means");
 265     newObj->nThresholds = oo->expectation->numOrdinal;
 266     newObj->expectedThresholds = oo->expectation->thresholds;
 267     // FIXME: threshold structure should be asked for by omxGetExpectationComponent
 268
 269         /* Read and set expected means, variances, and weights */
 270     cov = omxDataMatrix(dataMat, NULL);
 271     means = omxDataMeans(dataMat, NULL, NULL);
 272     weights = omxDataAcov(dataMat, NULL);
 273
 274     newObj->observedCov = cov;
 275     newObj->observedMeans = means;
 276     newObj->weights = weights;
 277     newObj->n = omxDataNumObs(dataMat);
 278     newObj->nThresholds = omxDataNumFactor(dataMat);
 279         UNPROTECT(1);
 280
 281         // Error Checking: Observed/Expected means must agree.
 282         // ^ is XOR: true when one is false and the other is not.
 283         if((newObj->expectedMeans == NULL) ^ (newObj->observedMeans == NULL)) {
 284             if(newObj->expectedMeans != NULL) {
 285                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
 286                             "Observed means not detected, but an expected means matrix was specified.\n  If you  wish to model the means, you must provide observed means.\n");
 287                     return;
 288             } else {
 289                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
 290                             "Observed means were provided, but an expected means matrix was not specified.\n  If you provide observed means, you must specify a model for the means.\n");
 291                     return;
 292             }
 293         }
 294
 295         if((newObj->expectedThresholds == NULL) ^ (newObj->observedThresholds == NULL)) {
 296             if(newObj->expectedMeans != NULL) {
 297                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
 298                             "Observed thresholds not detected, but an expected thresholds matrix was specified.\n   If you wish to model the thresholds, you must provide observed thresholds.\n ");
 299                     return;
 300             } else {
 301                     omxRaiseError(oo->matrix->currentState, OMX_ERROR,
 302                             "Observed thresholds were provided, but an expected thresholds matrix was not specified.\nIf you provide observed thresholds, you must specify a model for the thresholds.\n");
 303                     return;
 304             }
 305         }
 306
 307     /* Error check weight matrix size */
 308     int ncol = newObj->observedCov->cols;
 309     vectorSize = (ncol * (ncol + 1) ) / 2;
 310     if(newObj->expectedMeans != NULL) {
 311         vectorSize = vectorSize + ncol;
 312     }
 313     if(newObj->observedThresholds != NULL) {
 314         for(int i = 0; i < newObj->nThresholds; i++) {
 315             vectorSize = vectorSize + newObj->observedThresholds[i].numThresholds;
 316         }
 317     }
 318
 319     if(weights != NULL && weights->rows != weights->cols && weights->cols != vectorSize) {
 320         omxRaiseError(oo->matrix->currentState, OMX_DEVELOPER_ERROR,
 321          "Developer Error in WLS-based FitFunction object: WLS-based expectation specified an incorrectly-sized weight matrix.\nIf you are not developing a new expectation type, you should probably post this to the OpenMx forums.");
 322      return;
 323     }
 324
 325
 326         // FIXME: More error checking for incoming Fit Functions
 327
 328         /* Temporary storage for calculation */
 329         newObj->observedFlattened = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
 330         newObj->expectedFlattened = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
 331         newObj->P = omxInitMatrix(NULL, 1, vectorSize, TRUE, oo->matrix->currentState);
 332         newObj->B = omxInitMatrix(NULL, vectorSize, 1, TRUE, oo->matrix->currentState);
 333
 334     flattenDataToVector(newObj->observedCov, newObj->observedMeans, newObj->observedThresholds, newObj->nThresholds, newObj->observedFlattened);
 335     flattenDataToVector(newObj->expectedCov, newObj->expectedMeans, newObj->expectedThresholds, newObj->nThresholds, newObj->expectedFlattened);
 336
 337     oo->argStruct = (void*)newObj;
 338
 339 }