Updated copyright to 2013 for R/ demo/ models/passing and src/ folders, and also...
[openmx:openmx.git] / src / omxRFitFunction.c
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
28 void omxDestroyRFitFunction(omxFitFunction *off) {
29         UNPROTECT(5);                   // fitfun, model, flatModel, parameters, and state
30 }
31
32 static void omxCallRFitFunction(omxFitFunction *oo, int want, double *gradient) {
33         omx_omp_set_lock(&GlobalRLock);
34
35         omxState* currentState = oo->matrix->currentState;
36         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
37
38         SEXP theCall, theReturn;
39         PROTECT(theCall = allocVector(LANGSXP, 3));
40         SETCAR(theCall, rFitFunction->fitfun);
41         SETCADR(theCall, rFitFunction->model);
42         SETCADDR(theCall, rFitFunction->state);
43
44         PROTECT(theReturn = eval(theCall, R_GlobalEnv));
45
46         if (LENGTH(theReturn) < 1) {
47                 // seems impossible, but report it if it happens
48                 omxRaiseErrorf(currentState, "FitFunction returned nothing");
49         } else if (LENGTH(theReturn) == 1) {
50                 oo->matrix->data[0] = asReal(theReturn);
51         } else if (LENGTH(theReturn) == 2) {
52                 oo->matrix->data[0] = asReal(VECTOR_ELT(theReturn, 0));
53                 REPROTECT(rFitFunction->state = VECTOR_ELT(theReturn, 1), rFitFunction->stateIndex);
54         } else if (LENGTH(theReturn) > 2) {
55                 omxRaiseErrorf(currentState, "FitFunction returned more than 2 arguments");
56         }
57
58         UNPROTECT(2); // theCall and theReturn
59
60         omx_omp_unset_lock(&GlobalRLock);
61 }
62
63 void omxRepopulateRFitFunction(omxFitFunction* oo, double* x, int n) {
64         omx_omp_set_lock(&GlobalRLock);
65
66         omxRFitFunction* rFitFunction = (omxRFitFunction*)oo->argStruct;
67
68         SEXP theCall, estimate;
69
70         PROTECT(estimate = allocVector(REALSXP, n));
71         double *est = REAL(estimate);
72         for(int i = 0; i < n ; i++) {
73                 est[i] = x[i];
74         }
75
76         PROTECT(theCall = allocVector(LANGSXP, 5));
77
78         SETCAR(theCall, install("imxUpdateModelValues"));
79         SETCADR(theCall, rFitFunction->model);
80         SETCADDR(theCall, rFitFunction->flatModel);
81         SETCADDDR(theCall, rFitFunction->parameters);
82         SETCAD4R(theCall, estimate);
83
84         REPROTECT(rFitFunction->model = eval(theCall, R_GlobalEnv), rFitFunction->modelIndex);
85
86         UNPROTECT(2); // theCall, estimate
87         omx_omp_unset_lock(&GlobalRLock);
88
89         omxMarkDirty(oo->matrix);
90 }
91
92 omxRListElement* omxSetFinalReturnsRFitFunction(omxFitFunction *oo, int *numReturns) {
93         *numReturns = 1;
94         omxRListElement* retVal = (omxRListElement*) R_alloc(1, sizeof(omxRListElement));
95
96         retVal[0].numValues = 1;
97         retVal[0].values = (double*) R_alloc(1, sizeof(double));
98         strncpy(retVal[0].label, "Minus2LogLikelihood", 20);
99         retVal[0].values[0] = oo->matrix->data[0];
100
101         return retVal;
102 }
103
104 void omxInitRFitFunction(omxFitFunction* oo, SEXP rObj) {
105         if(OMX_DEBUG) { Rprintf("Initializing R fit function.\n"); }
106         omxRFitFunction *newObj = (omxRFitFunction*) R_alloc(1, sizeof(omxRFitFunction));
107         
108         /* Set Fit Function Calls to RFitFunction Calls */
109         oo->computeFun = omxCallRFitFunction;
110         oo->setFinalReturns = omxSetFinalReturnsRFitFunction;
111         oo->destructFun = omxDestroyRFitFunction;
112         oo->repopulateFun = omxRepopulateRFitFunction;
113         oo->argStruct = (void*) newObj;
114         
115         PROTECT(newObj->fitfun = GET_SLOT(rObj, install("fitfun")));
116         PROTECT_WITH_INDEX(newObj->model = GET_SLOT(rObj, install("model")), &(newObj->modelIndex));
117         PROTECT(newObj->flatModel = GET_SLOT(rObj, install("flatModel")));
118         PROTECT(newObj->parameters = GET_SLOT(rObj, install("parameters")));
119         PROTECT_WITH_INDEX(newObj->state = GET_SLOT(rObj, install("state")), &(newObj->stateIndex));
120
121 }
122
123