Remove commented PPML code from LISRELExpectation
[openmx:openmx.git] / src / omxLISRELExpectation.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 "omxExpectation.h"
18 #include "omxBLAS.h"
19 #include "omxFIMLFitFunction.h"
20 #include "omxLISRELExpectation.h"
21 #include "omxRAMExpectation.h"
22
23 extern void omxCreateMLFitFunction(omxFitFunction* oo, SEXP rObj, omxMatrix* cov, omxMatrix* means);
24 // TODO: Merge ML and FIML Fit Functions into one unit.
25
26 void omxCallLISRELExpectation(omxExpectation* oo) {
27     if(OMX_DEBUG) { Rprintf("LISREL Expectation Called.\n"); }
28         omxLISRELExpectation* oro = (omxLISRELExpectation*)(oo->argStruct);
29         
30         omxRecompute(oro->LX);
31         omxRecompute(oro->LY);
32         omxRecompute(oro->BE);
33         omxRecompute(oro->GA);
34         omxRecompute(oro->PH);
35         omxRecompute(oro->PS);
36         omxRecompute(oro->TD);
37         omxRecompute(oro->TE);
38         omxRecompute(oro->TH);
39         if(oro->TX != NULL) {     // Update means?
40                 omxRecompute(oro->TX);
41                 omxRecompute(oro->KA);
42         }
43         if(oro->TY != NULL) {
44                 omxRecompute(oro->TY);
45                 omxRecompute(oro->AL);
46         }
47         
48         omxCalculateLISRELCovarianceAndMeans(oro);
49 }
50
51 void omxDestroyLISRELExpectation(omxExpectation* oo) {
52
53         if(OMX_DEBUG) { Rprintf("Destroying LISREL Expectation.\n"); }
54         
55         omxLISRELExpectation* argStruct = (omxLISRELExpectation*)(oo->argStruct);
56
57         /* We allocated 'em, so we destroy 'em. */
58         if(argStruct->cov != NULL)
59                 omxFreeMatrixData(argStruct->cov);
60
61         if(argStruct->means != NULL)
62                 omxFreeMatrixData(argStruct->means);
63         
64         omxFreeMatrixData(argStruct->A);
65         omxFreeMatrixData(argStruct->B);
66         omxFreeMatrixData(argStruct->C);
67         omxFreeMatrixData(argStruct->D);
68         omxFreeMatrixData(argStruct->E);
69         omxFreeMatrixData(argStruct->F);
70         omxFreeMatrixData(argStruct->G);
71         omxFreeMatrixData(argStruct->H);
72         omxFreeMatrixData(argStruct->I);
73         omxFreeMatrixData(argStruct->J);
74         omxFreeMatrixData(argStruct->K);
75         omxFreeMatrixData(argStruct->L);
76         omxFreeMatrixData(argStruct->TOP);
77         omxFreeMatrixData(argStruct->BOT);
78         omxFreeMatrixData(argStruct->MUX);
79         omxFreeMatrixData(argStruct->MUY);
80 }
81
82 void omxPopulateLISRELAttributes(omxExpectation *oo, SEXP algebra) {
83     if(OMX_DEBUG) { Rprintf("Populating LISREL Attributes.  Currently this does very little!\n"); }
84         
85         /*
86         omxLISRELExpectation* oro = (omxLISRELExpectation*) (oo->argStruct);
87         omxMatrix* LX = oro->LX;
88         omxMatrix* LY = oro->LY;
89         omxMatrix* BE = oro->BE;
90         omxMatrix* GA = oro->GA;
91         omxMatrix* PH = oro->PH;
92         omxMatrix* PS = oro->PS;
93         omxMatrix* TD = oro->TD;
94         omxMatrix* TE = oro->TE;
95         omxMatrix* TH = oro->TH;
96         omxMatrix* LXPH = oro->LXPH;
97         omxMatrix* GAPH = oro->GAPH;
98         omxMatrix* W = oro->W;
99         omxMatrix* U = oro->U;
100         omxMatrix* I = oro->I;
101         int numIters = oro->numIters;
102         double oned = 1.0, zerod = 0.0;
103         
104         omxRecompute(LX);
105         omxRecompute(LY);
106         */ //This block of code works fine but because I do not use any of it later, it throws a huge block of warnings about unused variables.
107         // In general, I do not yet understand what this function needs to do.
108         
109         /*
110         omxFastRAMInverse(numIters, BE, Z, Ax, I ); // Z = (I-A)^-1
111         
112         if(OMX_DEBUG_ALGEBRA) { Rprintf("....DGEMM: %c %c \n %d %d %d \n %f \n %x %d %x %d \n %f %x %d.\n", *(Z->majority), *(S->majority), (Z->rows), (S->cols), (S->rows), oned, Z->data, (Z->leading), S->data, (S->leading), zerod, Ax->data, (Ax->leading));}
113         // F77_CALL(omxunsafedgemm)(Z->majority, S->majority, &(Z->rows), &(S->cols), &(S->rows), &oned, Z->data, &(Z->leading), S->data, &(S->leading), &zerod, Ax->data, &(Ax->leading));     // X = ZS
114         omxDGEMM(FALSE, FALSE, oned, Z, S, zerod, Ax);
115
116         if(OMX_DEBUG_ALGEBRA) { Rprintf("....DGEMM: %c %c %d %d %d %f %x %d %x %d %f %x %d.\n", *(Ax->majority), *(Z->minority), (Ax->rows), (Z->rows), (Z->cols), oned, X->data, (X->leading), Z->data, (Z->lagging), zerod, Ax->data, (Ax->leading));}
117         // F77_CALL(omxunsafedgemm)(Ax->majority, Z->minority, &(Ax->rows), &(Z->rows), &(Z->cols), &oned, Ax->data, &(Ax->leading), Z->data, &(Z->leading), &zerod, Ax->data, &(Ax->leading)); 
118         omxDGEMM(FALSE, TRUE, oned, Ax, Z, zerod, Ax);
119         // Ax = ZSZ' = Covariance matrix including latent variables
120         
121         SEXP expCovExt;
122         PROTECT(expCovExt = allocMatrix(REALSXP, Ax->rows, Ax->cols));
123         for(int row = 0; row < Ax->rows; row++)
124                 for(int col = 0; col < Ax->cols; col++)
125                         REAL(expCovExt)[col * Ax->rows + row] =
126                                 omxMatrixElement(Ax, row, col);
127         setAttrib(algebra, install("UnfilteredExpCov"), expCovExt);
128         UNPROTECT(1);
129         */
130 }
131
132 /* omxFastLISRELInverse would go here */
133
134
135 void omxCalculateLISRELCovarianceAndMeans(omxLISRELExpectation* oro) {
136         omxMatrix* LX = oro->LX;
137         omxMatrix* LY = oro->LY;
138         omxMatrix* BE = oro->BE;
139         omxMatrix* GA = oro->GA;
140         omxMatrix* PH = oro->PH;
141         omxMatrix* PS = oro->PS;
142         omxMatrix* TD = oro->TD;
143         omxMatrix* TE = oro->TE;
144         omxMatrix* TH = oro->TH;
145         omxMatrix* TX = oro->TX;
146         omxMatrix* TY = oro->TY;
147         omxMatrix* KA = oro->KA;
148         omxMatrix* AL = oro->AL;
149         omxMatrix* Cov = oro->cov;
150         omxMatrix* Means = oro->means;
151         int numIters = oro->numIters; //Used for fast RAM/LISREL inverse
152         omxMatrix* A = oro->A;
153         omxMatrix* B = oro->B;
154         omxMatrix* C = oro->C;
155         omxMatrix* D = oro->D;
156         omxMatrix* E = oro->E;
157         omxMatrix* F = oro->F;
158         omxMatrix* G = oro->G;
159         omxMatrix* H = oro->H;
160         omxMatrix* I = oro->I;
161         omxMatrix* J = oro->J;
162         omxMatrix* K = oro->K;
163         omxMatrix* L = oro->L;
164         omxMatrix* TOP = oro->TOP;
165         omxMatrix* BOT = oro->BOT;
166         omxMatrix* MUX = oro->MUX;
167         omxMatrix* MUY = oro->MUY;
168         omxMatrix** args = oro->args;
169         if(OMX_DEBUG) { Rprintf("Running LISREL computation in omxCalculateLISRELCovarianceAndMeans.\n"); }
170         double oned = 1.0, zerod=0.0; //, minusOned = -1.0;
171         //int ipiv[BE->rows], lwork = 4 * BE->rows * BE->cols; //This is copied from omxFastRAMInverse()
172         //double work[lwork];                                                                   // It lets you get the inverse of a matrix via omxDGETRI()
173         
174         
175         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(LX, "....LISREL: LX:");}
176         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(LY, "....LISREL: LY:");}
177         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(BE, "....LISREL: BE:");}
178         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(GA, "....LISREL: GA:");}
179         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(PH, "....LISREL: PH:");}
180         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(PS, "....LISREL: PS:");}
181         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(TD, "....LISREL: TD:");}
182         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(TE, "....LISREL: TE:");}
183         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(TH, "....LISREL: TH:");}
184         
185         /* Calculate the lower right quadrant: the covariance of the Xs */
186         if(LX->cols != 0 && LY->cols != 0) {
187         //if( (LX != NULL) && (LY != NULL) ) {
188                 if(OMX_DEBUG) {Rprintf("Calculating Lower Right Quadrant of Expected Covariance Matrix.\n"); }
189                 omxDGEMM(FALSE, FALSE, oned, LX, PH, zerod, A); // A = LX*PH
190                 omxCopyMatrix(B, TD); // B = TD
191                 omxDGEMM(FALSE, TRUE, oned, A, LX, oned, B);  // B = LX*PH*LX^T + TD
192                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(B, "....LISREL: Lower Right Quadrant of Model-implied Covariance Matrix:");}
193                 
194                 /* Calculate (I-BE)^(-1) and LY*(I-BE)^(-1) */
195                 if(OMX_DEBUG) {Rprintf("Calculating Inverse of I-BE.\n"); }
196                 omxFastRAMInverse(numIters, BE, C, L, I ); // C = (I-BE)^-1
197                 //omxCopyMatrix(C, BE); // C = BE
198                 //omxDGEMM(FALSE, FALSE, oned, I, I, minusOned, C); // C = I - BE
199                 //omxDGETRF(C, ipiv); //LU Decomp
200                 //omxDGETRI(C, ipiv, work, lwork); //Inverse based on LU Decomp ... C = C^(-1) = (I - BE)^(-1)
201                 
202                 
203                 omxDGEMM(FALSE, FALSE, oned, LY, C, zerod, D); // D = LY*C = LY * (I - BE)^(-1)
204                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(D, "....LISREL: LY*(I-BE)^(-1)");}
205                 
206                 /* Calculate the lower left quadrant: the covariance of Xs and Ys, nX by nY */
207                 if(OMX_DEBUG) {Rprintf("Calculating Lower Left Quadrant of Expected Covariance Matrix.\n"); }
208                 omxDGEMM(FALSE, TRUE, oned, A, GA, zerod, E); // E = A*GA^T = LX*PH*GA^T
209                 omxCopyMatrix(F, TH); // F = TH
210                 omxDGEMM(FALSE, TRUE, oned, E, D, oned, F); // F = E*D^T + F = LX*PH*GA^T * (LY * (I - BE)^(-1))^T + TH
211                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(F, "....LISREL: Lower Left Quadrant of Model-implied Covariance Matrix:");}
212         
213                 
214         /* Calculate the upper right quadrant: NOTE THIS IS MERELY THE LOWER LEFT QUADRANT TRANSPOSED. */
215         //DONE as omxTranspose(F)
216                 
217                 
218                 /* Calculate the upper left quadrant: the covariance of the Ys */
219                 if(OMX_DEBUG) {Rprintf("Calculating Upper Left Quadrant of Expected Covariance Matrix.\n"); }
220                 if(OMX_DEBUG_ALGEBRA) {Rprintf("Calculating G = GA*PH.\n");}
221                 omxDGEMM(FALSE, FALSE, oned, GA, PH, zerod, G); // G = GA*PH
222                 if(OMX_DEBUG_ALGEBRA) {Rprintf("Copying C = PS.\n");}
223                 omxCopyMatrix(C, PS); // C = PS
224                 if(OMX_DEBUG_ALGEBRA) {Rprintf("Calculating G = GA*PH.\n");}
225                 omxDGEMM(FALSE, TRUE, oned, G, GA, oned, C); // C = G*GA^T + C = GA*PH*GA^T + PS
226                 omxDGEMM(FALSE, FALSE, oned, D, C, zerod, H); // H = D*C = LY * (I - BE)^(-1) * (GA*PH*GA^T + PS)
227                 omxCopyMatrix(J, TE); // J = TE
228                 omxDGEMM(FALSE, TRUE, oned, H, D, oned, J); // J = H*D^T + J = LY * (I - BE)^(-1) * (GA*PH*GA^T + PS) * (LY * (I - BE)^(-1))^T + TE
229                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(J, "....LISREL: Upper Left Quadrant of Model-implied Covariance Matrix:");}
230         
231         
232         /* Construct the full model-implied covariance matrix from the blocks previously calculated */
233         // SigmaHat = ( J  t(F) )
234         //            ( F    B  )
235                 args[0] = F;
236                 args[1] = B;
237                 omxMatrixHorizCat(args, 2, BOT);
238                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(BOT, "....LISREL: BOT = cbind(F, B):");}
239                 args[0] = J;
240                 omxTransposeMatrix(F);
241                 args[1] = F;
242                 omxMatrixHorizCat(args, 2, TOP);
243                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(args[0], "....LISREL: TOP Debugging args[0] = J:");}
244                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(args[1], "....LISREL: TOP Debugging args[1] = F:");}
245                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(F, "....LISREL: TOP Debugging F (should be 2 rows):");}
246                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(TOP, "....LISREL: TOP = cbind(J, t(F)):");}
247                 omxTransposeMatrix(F); // So that it's back where it was.
248                 args[0] = TOP;
249                 args[1] = BOT;
250                 omxMatrixVertCat(args, 2, Cov);
251         
252                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(Cov, "....LISREL: Model-implied Covariance Matrix:");}
253         
254         
255                 /* Now Calculate the Expected Means */
256                 if(Means != NULL) {
257                         /* Mean of the Xs */
258                         //if(TX != NULL) {
259                                 omxCopyMatrix(MUX, TX);
260                                 omxDGEMV(FALSE, oned, LX, KA, oned, MUX);
261                         //}
262                         
263                         /* Mean of Ys */
264                         //if(TY != NULL) {
265                                 omxCopyMatrix(K, AL);
266                                 omxDGEMV(FALSE, oned, GA, KA, oned, K);
267                                 omxCopyMatrix(MUY, TY);
268                                 omxDGEMV(FALSE, oned, D, K, oned, MUY);
269                         //}
270                 
271                         /* Build means from blocks */
272                         args[0] = MUY;
273                         args[1] = MUX;
274                         omxMatrixVertCat(args, 2, Means);
275                         
276                         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(Means, "....LISREL: Model-implied Means Vector:");}
277                 }
278         }
279         else if(LX->cols != 0) { /* IF THE MODEL ONLY HAS Xs */
280         //else if(LX != NULL) { /* IF THE MODEL ONLY HAS Xs */
281                 if(OMX_DEBUG) {Rprintf("Calculating Lower Right Quadrant of Expected Covariance Matrix.\n"); }
282                 omxDGEMM(FALSE, FALSE, oned, LX, PH, zerod, A); // A = LX*PH
283                 omxCopyMatrix(Cov, TD); // Cov = TD
284                 omxDGEMM(FALSE, TRUE, oned, A, LX, oned, Cov);  // Cov = LX*PH*LX^T + Cov
285                 if(Means != NULL) {
286                                 /* Mean of the Xs */
287                                 omxCopyMatrix(Means, TX);
288                                 omxDGEMV(FALSE, oned, LX, KA, oned, Means);
289                 }
290         }
291         
292         /* IF THE MODEL ONLY HAS Ys */
293         else if(LY->cols != 0) {
294         //else if(LY != NULL) {
295                 /* Calculate (I-BE)^(-1) and LY*(I-BE)^(-1) */
296                 if(OMX_DEBUG) {Rprintf("Calculating Inverse of I-BE.\n"); }
297                 omxFastRAMInverse(numIters, BE, C, L, I ); // C = (I-BE)^-1
298                 //omxCopyMatrix(C, BE); // C = BE
299                 //omxDGEMM(FALSE, FALSE, oned, I, I, minusOned, C); // C = I - BE
300                 //omxDGETRF(C, ipiv); //LU Decomp
301                 //omxDGETRI(C, ipiv, work, lwork); //Inverse based on LU Decomp ... C = C^(-1) = (I - BE)^(-1)
302                 omxDGEMM(FALSE, FALSE, oned, LY, C, zerod, D); // D = LY*C = LY * (I - BE)^(-1)
303                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(D, "....LISREL: LY*(I-BE)^(-1)");}
304                 /* Calculate the upper left quadrant: the covariance of the Ys */
305                 if(OMX_DEBUG) {Rprintf("Calculating Upper Left Quadrant of Expected Covariance Matrix.\n"); }
306                 if(OMX_DEBUG_ALGEBRA) {Rprintf("Copying C = PS.\n");}
307                 omxDGEMM(FALSE, FALSE, oned, D, PS, zerod, H); // H = D*PS = LY * (I - BE)^(-1) * PS
308                 omxCopyMatrix(Cov, TE); // Cov = TE
309                 omxDGEMM(FALSE, TRUE, oned, H, D, oned, Cov); // Cov = H*D^T + Cov = LY * (I - BE)^(-1) * PS * (LY * (I - BE)^(-1))^T + TE
310                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(J, "....LISREL: Upper Left Quadrant of Model-implied Covariance Matrix:");}
311                 if(Means != NULL) {
312                                 omxCopyMatrix(Means, TY);
313                                 omxDGEMV(FALSE, oned, D, AL, oned, Means);
314                 }
315         }
316 /*      
317         if(OMX_DEBUG) { Rprintf("Running RAM computation."); }
318                 
319         double oned = 1.0, zerod=0.0;
320         
321         if(Ax == NULL || I == NULL || Z == NULL || Y == NULL || X == NULL) {
322                 error("Internal Error: RAM Metadata improperly populated.  Please report this to the OpenMx development team.");
323         }
324                 
325         if(Cov == NULL && Means == NULL) {
326                 return; // We're not populating anything, so why bother running the calculation?
327         }
328         
329         // if(   (Cov->rows != Cov->cols)  || (A->rows  != A->cols)  // Conformance check
330         //      || (X->rows  != Cov->cols)  || (X->cols  != A->rows)
331         //      || (Y->rows  != Cov->cols)  || (Y->cols  != A->rows)
332         //      || (Ax->rows != Cov->cols)  || (Ax->cols != A->rows)
333         //      || (I->rows  != Cov->cols)  || (I->cols  != Cov->rows)
334         //      || (Y->rows  != Cov->cols)  || (Y->cols  != A->rows)
335         //      || (M->cols  != Cov->cols)  || (M->rows  != 1)
336         //      || (Means->rows != 1)       || (Means->cols != Cov->cols) ) {
337         //              error("INTERNAL ERROR: Incorrectly sized matrices being passed to omxRAMObjective Calculation.\n Please report this to the OpenMx development team.");
338         // }
339         
340         omxFastRAMInverse(numIters, A, Z, Ax, I );
341                 
342         // IMPORTANT: Cov = FZSZ'F'
343         if(OMX_DEBUG_ALGEBRA) { Rprintf("....DGEMM: %c %c %d %d %d %f %x %d %x %d %f %x %d.\n", *(F->majority), *(Z->majority), (F->rows), (Z->cols), (Z->rows), oned, F->data, (F->leading), Z->data, (Z->leading), zerod, Y->data, (Y->leading));}
344         // F77_CALL(omxunsafedgemm)(F->majority, Z->majority, &(F->rows), &(Z->cols), &(Z->rows), &oned, F->data, &(F->leading), Z->data, &(Z->leading), &zerod, Y->data, &(Y->leading));       // Y = FZ
345         omxDGEMM(FALSE, FALSE, 1.0, F, Z, 0.0, Y);
346
347         if(OMX_DEBUG_ALGEBRA) { Rprintf("....DGEMM: %c %c %d %d %d %f %x %d %x %d %f %x %d.\n", *(Y->majority), *(S->majority), (Y->rows), (S->cols), (S->rows), oned, Y->data, (Y->leading), S->data, (S->leading), zerod, X->data, (X->leading));}
348         // F77_CALL(omxunsafedgemm)(Y->majority, S->majority, &(Y->rows), &(S->cols), &(S->rows), &oned, Y->data, &(Y->leading), S->data, &(S->leading), &zerod, X->data, &(X->leading));       // X = FZS
349         omxDGEMM(FALSE, FALSE, 1.0, Y, S, 0.0, X);
350
351         if(OMX_DEBUG_ALGEBRA) { Rprintf("....DGEMM: %c %c %d %d %d %f %x %d %x %d %f %x %d.\n", *(X->majority), *(Y->minority), (X->rows), (Y->rows), (Y->cols), oned, X->data, (X->leading), Y->data, (Y->lagging), zerod, Cov->data, (Cov->leading));}
352         // F77_CALL(omxunsafedgemm)(X->majority, Y->minority, &(X->rows), &(Y->rows), &(Y->cols), &oned, X->data, &(X->leading), Y->data, &(Y->leading), &zerod, Cov->data, &(Cov->leading));
353         omxDGEMM(FALSE, TRUE, 1.0, X, Y, 0.0, Cov);
354          // Cov = FZSZ'F' (Because (FZ)' = Z'F')
355         
356         if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(Cov, "....RAM: Model-implied Covariance Matrix:");}
357         
358         if(M != NULL && Means != NULL) {
359                 // F77_CALL(omxunsafedgemv)(Y->majority, &(Y->rows), &(Y->cols), &oned, Y->data, &(Y->leading), M->data, &onei, &zerod, Means->data, &onei);
360                 omxDGEMV(FALSE, 1.0, Y, M, 0.0, Means);
361                 if(OMX_DEBUG_ALGEBRA) {omxPrintMatrix(Means, "....RAM: Model-implied Means Vector:");}
362         }
363 */
364 }
365
366 void omxInitLISRELExpectation(omxExpectation* oo, SEXP rObj) {
367         
368         if(OMX_DEBUG) { Rprintf("Initializing LISREL Expectation.\n"); }
369                 
370         int nx, nxi, ny, neta, ntotal;
371         
372         SEXP slotValue;
373         
374         /* Create and fill expectation */
375         oo->expType = "omxLISRELExpectation";
376         omxLISRELExpectation *LISobj = (omxLISRELExpectation*) R_alloc(1, sizeof(omxLISRELExpectation));
377         omxState* currentState = oo->currentState;
378         
379         /* Set Expectation Calls and Structures */
380         oo->computeFun = omxCallLISRELExpectation;
381         oo->destructFun = omxDestroyLISRELExpectation;
382         oo->componentFun = omxGetLISRELExpectationComponent;
383         oo->populateAttrFun = omxPopulateLISRELAttributes;
384         oo->argStruct = (void*) LISobj;
385         
386         /* Set up expectation structures */
387         if(OMX_DEBUG) { Rprintf("Initializing LISREL Meta Data for expectation.\n"); }
388         
389         if(OMX_DEBUG) { Rprintf("Processing LX.\n"); }
390         LISobj->LX = omxNewMatrixFromIndexSlot(rObj, currentState, "LX");
391         
392         if(OMX_DEBUG) { Rprintf("Processing LY.\n"); }
393         LISobj->LY = omxNewMatrixFromIndexSlot(rObj, currentState, "LY");
394         
395         if(OMX_DEBUG) { Rprintf("Processing BE.\n"); }
396         LISobj->BE = omxNewMatrixFromIndexSlot(rObj, currentState, "BE");
397         
398         if(OMX_DEBUG) { Rprintf("Processing GA.\n"); }
399         LISobj->GA = omxNewMatrixFromIndexSlot(rObj, currentState, "GA");
400         
401         if(OMX_DEBUG) { Rprintf("Processing PH.\n"); }
402         LISobj->PH = omxNewMatrixFromIndexSlot(rObj, currentState, "PH");
403         
404         if(OMX_DEBUG) { Rprintf("Processing PS.\n"); }
405         LISobj->PS = omxNewMatrixFromIndexSlot(rObj, currentState, "PS");
406         
407         if(OMX_DEBUG) { Rprintf("Processing TD.\n"); }
408         LISobj->TD = omxNewMatrixFromIndexSlot(rObj, currentState, "TD");
409         
410         if(OMX_DEBUG) { Rprintf("Processing TE.\n"); }
411         LISobj->TE = omxNewMatrixFromIndexSlot(rObj, currentState, "TE");
412         
413         if(OMX_DEBUG) { Rprintf("Processing TH.\n"); }
414         LISobj->TH = omxNewMatrixFromIndexSlot(rObj, currentState, "TH");
415
416         if(OMX_DEBUG) { Rprintf("Processing TX.\n"); }
417         LISobj->TX = omxNewMatrixFromIndexSlot(rObj, currentState, "TX");
418
419         if(OMX_DEBUG) { Rprintf("Processing TY.\n"); }
420         LISobj->TY = omxNewMatrixFromIndexSlot(rObj, currentState, "TY");
421
422         if(OMX_DEBUG) { Rprintf("Processing KA.\n"); }
423         LISobj->KA = omxNewMatrixFromIndexSlot(rObj, currentState, "KA");
424
425         if(OMX_DEBUG) { Rprintf("Processing AL.\n"); }
426         LISobj->AL = omxNewMatrixFromIndexSlot(rObj, currentState, "AL");
427         
428         
429         if(LISobj->LY == NULL) {
430                 LISobj->LY = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
431                 LISobj->PS = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
432                 LISobj->BE = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
433                 LISobj->TE = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
434         }
435         
436         if(LISobj->LX == NULL) {
437                 LISobj->LX = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
438                 LISobj->PH = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
439                 LISobj->TD = omxInitMatrix(NULL, 0, 0, TRUE, currentState);
440         }
441         
442         if(LISobj->LY->cols == 0 || LISobj->LX->cols == 0) {
443                 LISobj->GA = omxInitMatrix(NULL, LISobj->LY->cols, LISobj->LX->cols, TRUE, currentState);
444                 LISobj->TH = omxInitMatrix(NULL, LISobj->LX->rows, LISobj->LY->rows, TRUE, currentState);
445         }
446         
447         
448         /* Identity Matrix, Size Of BE */
449         if(OMX_DEBUG) { Rprintf("Generating I.\n"); }
450         LISobj->I = omxNewIdentityMatrix(LISobj->BE->rows, currentState);
451         omxRecompute(LISobj->I);
452         
453         
454         /* Get the nilpotency index of the BE matrix for I-BE inverse speedup */
455         if(OMX_DEBUG) { Rprintf("Processing expansion iteration depth.\n"); }
456         PROTECT(slotValue = GET_SLOT(rObj, install("depth")));
457         LISobj->numIters = INTEGER(slotValue)[0];
458         if(OMX_DEBUG) { Rprintf("Using %d iterations.", LISobj->numIters); }
459         UNPROTECT(1);
460         
461         
462         /* Initialize the place holder matrices used in calculations */
463         nx = LISobj->LX->rows;
464         nxi = LISobj->LX->cols;
465         ny = LISobj->LY->rows;
466         neta = LISobj->LY->cols;
467         ntotal = nx + ny;
468         
469         
470         if(OMX_DEBUG) { Rprintf("Generating internals for computation.\n"); }
471         
472         LISobj->A =     omxInitMatrix(NULL, nx, nxi, TRUE, currentState);
473         LISobj->B =     omxInitMatrix(NULL, nx, nx, TRUE, currentState);
474         LISobj->C =     omxInitMatrix(NULL, neta, neta, TRUE, currentState);
475         LISobj->D =     omxInitMatrix(NULL, ny, neta, TRUE, currentState);
476         LISobj->E =     omxInitMatrix(NULL, nx, neta, TRUE, currentState);
477         LISobj->F =     omxInitMatrix(NULL, nx, ny, TRUE, currentState);
478         LISobj->G =     omxInitMatrix(NULL, neta, nxi, TRUE, currentState);
479         LISobj->H =     omxInitMatrix(NULL, ny, neta, TRUE, currentState);
480         LISobj->J =     omxInitMatrix(NULL, ny, ny, TRUE, currentState);
481         LISobj->K =     omxInitMatrix(NULL, neta, 1, TRUE, currentState);
482         LISobj->L =     omxInitMatrix(NULL, neta, neta, TRUE, currentState);
483         LISobj->TOP =   omxInitMatrix(NULL, ny, ntotal, TRUE, currentState);
484         LISobj->BOT =   omxInitMatrix(NULL, nx, ntotal, TRUE, currentState);
485         LISobj->MUX =   omxInitMatrix(NULL, nx, 1, TRUE, currentState);
486         LISobj->MUY =   omxInitMatrix(NULL, ny, 1, TRUE, currentState);
487         
488         
489         LISobj->cov =   omxInitMatrix(NULL, ntotal, ntotal, TRUE, currentState);
490
491         LISobj->args = (omxMatrix**) R_alloc(2, sizeof(omxMatrix*));
492         
493         /* Means */
494         if(LISobj->TX != NULL || LISobj->TY != NULL || LISobj->KA != NULL || LISobj->AL != NULL) {
495                 LISobj->means =         omxInitMatrix(NULL, 1, ntotal, TRUE, currentState);
496         } else LISobj->means  =         NULL;
497         //TODO: Adjust means processing to allow only Xs or only Ys
498
499 }
500
501 omxMatrix* omxGetLISRELExpectationComponent(omxExpectation* ox, omxFitFunction* off, const char* component) {
502         omxLISRELExpectation* ore = (omxLISRELExpectation*)(ox->argStruct);
503         omxMatrix* retval = NULL;
504
505         if(!strncmp("cov", component, 3)) {
506                 retval = ore->cov;
507         } else if(!strncmp("mean", component, 4)) {
508                 retval = ore->means;
509         } else if(!strncmp("pvec", component, 4)) {
510                 // Once implemented, change compute function and return pvec
511         }
512         
513         if(OMX_DEBUG) { Rprintf("Returning 0x%x.\n", retval); }
514
515         return retval;
516 }