ifa: Fixup
[openmx:openmx.git] / src / omxRFitFunction.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 "omxRFitFunction.h"
25 #include "omxOpenmpWrap.h"
26 #include "npsolWrap.h"
27 #include "Compute.h"
28
29 void omxDestroyRFitFunction(omxFitFunction *off) {
30         UNPROTECT(4);                   // fitfun, model, flatModel, and state
31 }
32
33 static void omxCallRFitFunction(omxFitFunction *oo, int want, double *gradient) {
34         omxState* currentState = oo->matrix->currentState;
35         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
36
37         SEXP theCall, theReturn;
38         PROTECT(theCall = allocVector(LANGSXP, 3));
39         SETCAR(theCall, rFitFunction->fitfun);
40         SETCADR(theCall, rFitFunction->model);
41         SETCADDR(theCall, rFitFunction->state);
42
43         PROTECT(theReturn = eval(theCall, R_GlobalEnv));
44
45         if (LENGTH(theReturn) < 1) {
46                 // seems impossible, but report it if it happens
47                 omxRaiseErrorf(currentState, "FitFunction returned nothing");
48         } else if (LENGTH(theReturn) == 1) {
49                 oo->matrix->data[0] = asReal(theReturn);
50         } else if (LENGTH(theReturn) == 2) {
51                 oo->matrix->data[0] = asReal(VECTOR_ELT(theReturn, 0));
52                 REPROTECT(rFitFunction->state = VECTOR_ELT(theReturn, 1), rFitFunction->stateIndex);
53         } else if (LENGTH(theReturn) > 2) {
54                 omxRaiseErrorf(currentState, "FitFunction returned more than 2 arguments");
55         }
56
57         UNPROTECT(2); // theCall and theReturn
58 }
59
60 omxRListElement* omxSetFinalReturnsRFitFunction(omxFitFunction *oo, int *numReturns) {
61         *numReturns = 1;
62         omxRListElement* retVal = (omxRListElement*) R_alloc(1, sizeof(omxRListElement));
63
64         retVal[0].numValues = 1;
65         retVal[0].values = (double*) R_alloc(1, sizeof(double));
66         strncpy(retVal[0].label, "Minus2LogLikelihood", 20);
67         retVal[0].values[0] = oo->matrix->data[0];
68
69         return retVal;
70 }
71
72 void omxInitRFitFunction(omxFitFunction* oo) {
73         FitContext::setRFitFunction(oo);
74
75         if(OMX_DEBUG) { mxLog("Initializing R fit function."); }
76         omxRFitFunction *newObj = (omxRFitFunction*) R_alloc(1, sizeof(omxRFitFunction));
77         
78         SEXP rObj = oo->rObj;
79
80         /* Set Fit Function Calls to RFitFunction Calls */
81         oo->computeFun = omxCallRFitFunction;
82         oo->setFinalReturns = omxSetFinalReturnsRFitFunction;
83         oo->destructFun = omxDestroyRFitFunction;
84         oo->argStruct = (void*) newObj;
85         
86         PROTECT(newObj->fitfun = GET_SLOT(rObj, install("fitfun")));
87         PROTECT_WITH_INDEX(newObj->model = GET_SLOT(rObj, install("model")), &(newObj->modelIndex));
88         PROTECT(newObj->flatModel = GET_SLOT(rObj, install("flatModel")));
89         PROTECT_WITH_INDEX(newObj->state = GET_SLOT(rObj, install("state")), &(newObj->stateIndex));
90
91 }
92
93