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