PageRenderTime 146ms CodeModel.GetById 10ms app.highlight 126ms RepoModel.GetById 1ms app.codeStats 0ms

/src/C/factor_model_util.c

http://github.com/beechung/Latent-Factor-Models
C | 1148 lines | 827 code | 170 blank | 151 comment | 429 complexity | fdd3e70f2ce58d51fe86ff8ea227cf07 MD5 | raw file
   1/*
   2	Copyright (c) 2011, Yahoo! Inc.  All rights reserved.
   3	Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
   4
   5    Author: Bee-Chung Chen
   6*/
   7
   8#include <R.h>
   9#include <Rinternals.h>
  10#include <Rmath.h>
  11#include <R_ext/Lapack.h>
  12#include <stdio.h>
  13#include <time.h>
  14#include "util.h"
  15
  16
  17void gaussianPosterior_mainEffect(
  18    // OUTPUT
  19    double* outSample, double* outMean, double* outVar,
  20    //INPUT
  21    const int* option /*1:Sample, 2:Mean&Var, 3:Sample&Mean&Var*/,
  22    const int* thisEffIndex /*user or item*/,
  23    const double* rest /*o in the paper*/,
  24    const double* fittedEff /*g0*x_user or d0*x_item*/,
  25    const double* multiplier /*NULL*/,
  26    const double* var_y, const double* var_eff /*var_alpha or var_beta*/,
  27    const int* nObs, const int* nThisEff, const int* nVar_y, const int* nVar_eff,
  28    // OTHER
  29    const int* debug
  30){
  31    int outputSample=0, outputMeanVar=0;
  32
  33    if(*option == 1){
  34        outputSample = 1; outputMeanVar = 0;
  35    }else if(*option == 2){
  36        outputSample = 0; outputMeanVar = 1;
  37    }else if(*option == 3){
  38        outputSample = 1; outputMeanVar = 1;
  39    }else error("Unknown option: %d", *option);
  40
  41    double *sum_ivar = (double*)Calloc(*nThisEff, double); // sum_ivar[effect_i] = sum_{k having effect_i} 1/var_y[k]
  42    double *sum_o    = (double*)Calloc(*nThisEff, double); // sum_o[effect_i] = sum_{k having effect_i} o[k]/var_y[k]
  43
  44    for(int k=0; k<*nObs; k++){
  45        const int thisIndex  = thisEffIndex[k]-1;
  46        if(*debug > 0) CHK_C_INDEX(thisIndex,  *nThisEff);
  47
  48        double o = rest[k]; // for alpha and beta: o.  for gamma: o * x_dyad * b
  49        double square = 1;  // for alpha and beta: 1.  for gamma: (x_dyad * b)^2
  50        if(multiplier != NULL){
  51            o *= multiplier[k];
  52            square = multiplier[k] * multiplier[k];
  53        }
  54
  55        if((*nVar_y) == 1){
  56            sum_ivar[thisIndex] += square / var_y[0];
  57            sum_o[   thisIndex] +=      o / var_y[0];
  58        }else if((*nVar_y) == (*nObs)){
  59            sum_ivar[thisIndex] += square / var_y[k];
  60            sum_o[   thisIndex] +=      o / var_y[k];
  61        }else error("nVar_y = %d, nObs = %d", *nVar_y, *nObs);
  62    }
  63
  64    if(outputSample) GetRNGstate();
  65    for(int i=0; i<*nThisEff; i++){
  66    	double mean = 0, var=0;
  67        if((*nVar_eff) == 1){
  68            var  = 1.0 / (sum_ivar[i] + (1.0 / var_eff[0]));
  69            mean = var * (sum_o[i] + (fittedEff[i] / var_eff[0]));
  70        }else if((*nVar_eff) == (*nThisEff)){
  71            var  = 1.0 / (sum_ivar[i] + (1.0 / var_eff[i]));
  72            mean = var * (sum_o[i] + (fittedEff[i] / var_eff[i]));
  73        }else error("nVar_eff = %d, nThisEff = %d", *nVar_eff, *nThisEff);
  74
  75        if(outputMeanVar){
  76            outMean[i] = mean;
  77            outVar[i]  = var;
  78        }
  79        if(outputSample){
  80            outSample[i] = rnorm(mean, sqrt(var));
  81        }
  82    }
  83    if(outputSample) PutRNGstate();
  84
  85    Free(sum_ivar);
  86    Free(sum_o);
  87}
  88
  89SEXP gaussianPosterior_mainEffect_Call(
  90    // OUTPUT
  91    SEXP outSample, SEXP outMean, SEXP outVar,
  92    //INPUT
  93    SEXP option /*1:Sample, 2:Mean&Var, 3:Sample&Mean&Var*/,
  94    SEXP thisEffIndex /*user or item*/,
  95    SEXP rest /*o in the paper*/,
  96    SEXP fittedEff /*g0*x_user or d0*x_item*/,
  97    SEXP multiplier /*NULL*/,
  98    SEXP var_y, SEXP var_eff /*var_alpha or var_beta*/,
  99    SEXP nObs, SEXP nThisEff, SEXP nVar_y, SEXP nVar_eff,
 100    // OTHER
 101    SEXP debug
 102){
 103  gaussianPosterior_mainEffect(
 104      // OUTPUT
 105      MY_REAL(outSample), MY_REAL(outMean), MY_REAL(outVar),
 106      //INPUT
 107      MY_INTEGER(option) /*1:Sample, 2:Mean&Var, 3:Sample&Mean&Var*/,
 108      MY_INTEGER(thisEffIndex) /*user or item*/,
 109      MY_REAL(rest) /*o in the paper*/,
 110      MY_REAL(fittedEff) /*g0*x_user or d0*x_item*/,
 111      MY_REAL(multiplier) /*NULL*/,
 112      MY_REAL(var_y), MY_REAL(var_eff) /*var_alpha or var_beta*/,
 113      MY_INTEGER(nObs), MY_INTEGER(nThisEff), MY_INTEGER(nVar_y), MY_INTEGER(nVar_eff),
 114      // OTHER
 115      MY_INTEGER(debug)
 116  );
 117  return R_NilValue;
 118}
 119
 120void gaussianPosterior_2WayInteraction(
 121    // OUTPUT
 122    double* sample, double* posteriorMean, double* posteriorVar,
 123    // INPUT
 124    const int* option /*1:Sample, 2:Mean&Var, 3:Sample&Mean&Var*/,
 125    const int* thisEffIndex /*user or item*/, const int* otherEffIndex /*item or user*/, const double* obs /*o in the paper*/,
 126    const double* priorMean /*Gw or Dz*/, const double* otherEff /*v or u*/,
 127    const double* obsVar, const double* priorVar /*var_u or var_v*/,
 128    const int* nObs, const int* nLevelsThisEff, const int* nLevelsOtherEff,
 129    const int* nFactors, const int* nObsVar, const int *nPriorVar,
 130    const int* obsIndex, const int* oiStart, const int* oiNum,
 131    // OTHER
 132    const int* debug
 133){
 134    double *sum_vv /*nFactors x nFactors*/, *sum_ov /*nFactors x 1*/, o, *vj /*nFactors x 1*/, *Gwi /*nFactors x 1*/,
 135           *work, size, *mean, *var, *temp, *rnd;
 136    int i,j,k,m, thisIndex, otherIndex, outputSample, outputMeanVar, oIndex, lwork=-1, info, symCheck;
 137    char jobz = 'V', uplo = 'L';
 138
 139    symCheck = (*debug) - 2;
 140
 141    if(*option == 1){
 142        outputSample = 1; outputMeanVar = 0;
 143    }else if(*option == 2){
 144        outputSample = 0; outputMeanVar = 1;
 145    }else if(*option == 3){
 146        outputSample = 1; outputMeanVar = 1;
 147    }else error("Unknown option: %d", *option);
 148
 149    vj     = (double*)Calloc(*nFactors, double);
 150    Gwi    = (double*)Calloc(*nFactors, double);
 151    sum_ov = (double*)Calloc(*nFactors, double);
 152    sum_vv = (double*)Calloc((*nFactors)*(*nFactors), double);
 153    mean   = (double*)Calloc(*nFactors, double);
 154    var    = (double*)Calloc((*nFactors)*(*nFactors), double);
 155    temp   = (double*)Calloc((*nFactors)*(*nFactors), double);
 156
 157    if(outputSample) GetRNGstate();
 158
 159    for(i=0; i<*nLevelsThisEff; i++){
 160
 161        thisIndex = i+1;
 162        for(j=0; j<*nFactors; j++) sum_ov[j] = 0;
 163        for(j=0; j<(*nFactors)*(*nFactors); j++) sum_vv[j] = 0;
 164
 165        for(j=0; j<oiNum[i]; j++){
 166
 167            oIndex = obsIndex[R_VEC(oiStart[i]+j)];
 168            otherIndex = otherEffIndex[R_VEC(oIndex)];
 169
 170            if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
 171            if(*debug > 0) CHK_R_INDEX(otherIndex, *nLevelsOtherEff);
 172            if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
 173
 174            o = obs[R_VEC(oIndex)];
 175            for(k=1; k<=*nFactors; k++) vj[R_VEC(k)] = otherEff[R_MAT(otherIndex,k,*nLevelsOtherEff)];
 176
 177            double var_y_thisObs = 0;
 178            if((*nObsVar) == 1)            var_y_thisObs = obsVar[0];
 179            else if((*nObsVar) == (*nObs)) var_y_thisObs = obsVar[R_VEC(oIndex)];
 180            else error("nVar_y = %d, nObs = %d", *nObsVar, *nObs);
 181
 182            for(k=0; k<*nFactors; k++) sum_ov[k] += (o * vj[k]) / var_y_thisObs;
 183            for(k=0; k<*nFactors; k++)
 184                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += (vj[k] * vj[m]) / var_y_thisObs;
 185        }
 186
 187        for(k=1; k<=*nFactors; k++) Gwi[R_VEC(k)] = priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)];
 188
 189        if((*nPriorVar) == 1){
 190            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[0]);
 191            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[0];
 192        }else if((*nPriorVar) == (*nLevelsThisEff)*(*nFactors)*(*nFactors)){
 193
 194            for(k=0; k<*nFactors; k++)
 195                for(m=0; m<*nFactors; m++) temp[C_MAT(k,m,*nFactors)] = priorVar[C_3DA(i,k,m,*nLevelsThisEff,*nFactors)];
 196            sym_inv_byCholesky(temp, nFactors, &symCheck);
 197            // Now, temp is var_u[i]^-1
 198
 199            for(k=0; k<*nFactors; k++)
 200                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += temp[C_MAT(k,m,*nFactors)];
 201
 202            // reuse the space of vj
 203            for(k=0; k<*nFactors; k++) vj[k] = Gwi[k];
 204            for(k=0; k<*nFactors; k++){
 205                Gwi[k] = 0;
 206                for(m=0; m<*nFactors; m++) Gwi[k] += temp[C_MAT(k,m,*nFactors)] * vj[m];
 207            }
 208
 209        }else error("nVar_eff = %d, nThisEff = %d", *nPriorVar, *nLevelsThisEff);
 210
 211        // Now, sum_vv = var^-1
 212        //      Gwi = var_u[i]^-1 G wi
 213
 214        if((*nFactors) == 1){
 215            var[0]  = 1/sum_vv[0];
 216            mean[0] = var[0] * (sum_ov[0] + Gwi[0]);
 217            if(outputSample){
 218                sample[i] = rnorm(mean[0], sqrt(var[0]));
 219            }
 220            if(outputMeanVar){
 221                posteriorMean[i] = mean[0];
 222                posteriorVar[i]  = var[0];
 223            }
 224        }else{
 225            double *eigen_val, *eigen_vec;
 226            eigen_val = vj;
 227            eigen_vec = sum_vv;
 228
 229            if(*debug > 2) CHK_SYMMETRIC(eigen_vec, *nFactors);
 230
 231            //
 232            // Compute the variance-covariance matrix
 233            //
 234            // Allocate workspace for eigen value decomposition (only allocate once)
 235            if(lwork == -1){
 236                F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, &size, &lwork, &info);
 237                if(info != 0) error("error in dsyev(...)");
 238                lwork = (int)size;
 239                work  = (double*)Calloc(lwork, double);
 240            }
 241
 242            F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, work, &lwork, &info);
 243            if(info != 0) error("error in dsyev(...)");
 244            for(j=0; j<*nFactors; j++) eigen_val[j] = 1/eigen_val[j];
 245            // Now, eigen_val, eigen_vec are the eigen values and vectors of var
 246
 247            for(j=0; j<*nFactors; j++)
 248                for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 249            for(j=0; j<*nFactors; j++)
 250                for(k=0; k<*nFactors; k++){
 251                    var[C_MAT(j,k,*nFactors)] = 0;
 252                    for(m=0; m<*nFactors; m++) var[C_MAT(j,k,*nFactors)] += temp[C_MAT(j,m,*nFactors)] * eigen_vec[C_MAT(k,m,*nFactors)];
 253                }
 254            // Now, var is the variance-covariance matrix
 255
 256            //
 257            // Compute the mean vector
 258            //
 259
 260            // print_vector("sum_ov: ", sum_ov, *nFactors);
 261            // print_vector("Gwi: ", Gwi, *nFactors);
 262
 263            for(j=0; j<*nFactors; j++) sum_ov[j] += Gwi[j];
 264            for(j=0; j<*nFactors; j++){
 265                mean[j] = 0;
 266                for(k=0; k<*nFactors; k++) mean[j] += var[C_MAT(j,k,*nFactors)] * sum_ov[k];
 267            }
 268
 269            if(outputMeanVar){
 270                for(j=0; j<*nFactors; j++) posteriorMean[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 271                for(j=0; j<*nFactors; j++)
 272                    for(k=0; k<*nFactors; k++) posteriorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)] = var[C_MAT(j,k,*nFactors)];
 273            }
 274
 275            // DEBUG CODE
 276            // if(i==38){
 277        	//     printf("i = %d\n",i);
 278            //     print_vector("mean: ", mean, *nFactors);
 279            //     print_matrix(" var: ", var, *nFactors, *nFactors);
 280            // }
 281
 282            if(*debug >= 2 && oiNum[i] == 0){
 283                for(k=1; k<=*nFactors; k++){
 284                    CHK_SAME_NUMBER("mean[k] != fittedEff[k]", mean[R_VEC(k)], priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)]);
 285                }
 286                if((*nPriorVar) == 1){
 287                    for(j=0; j<*nFactors; j++)
 288                        for(k=0; k<*nFactors; k++){
 289                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[0]);}
 290                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 291                        }
 292                }else{
 293                    for(j=0; j<*nFactors; j++)
 294                        for(k=0; k<*nFactors; k++){
 295                            CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)]);
 296                        }
 297                }
 298            }
 299            //
 300            //  Generate the random vector
 301            //
 302            if(outputSample){
 303                // reuse the space allocated for Gwi
 304                rnd = Gwi;
 305
 306                // DEBUG CODE
 307                // if(i==38){
 308                //     print_vector("eigenvalue: ", eigen_val, *nFactors);
 309                //     Rprintf("eigenvector:\n");  print_matrix("    ", eigen_vec, *nFactors, *nFactors);
 310                // }
 311
 312                for(j=0; j<*nFactors; j++) eigen_val[j] = sqrt(eigen_val[j]);
 313                for(j=0; j<*nFactors; j++)
 314                    for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 315
 316                for(j=0; j<*nFactors; j++) rnd[j] = norm_rand();
 317
 318                // DEBUG CODE
 319                // if(i==38){
 320                //     Rprintf("temp:\n");
 321                //     print_matrix("    ", temp, *nFactors, *nFactors);
 322                //     print_vector("rnd: ", rnd, *nFactors);
 323                // }
 324
 325                for(j=0; j<*nFactors; j++){
 326                    sample[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 327                    for(k=0; k<*nFactors; k++) sample[C_MAT(i,j,*nLevelsThisEff)] += temp[C_MAT(j,k,*nFactors)] * rnd[k];
 328                }
 329            }
 330        }
 331    }
 332    if(outputSample) PutRNGstate();
 333
 334    Free(sum_ov);
 335    Free(sum_vv);
 336    Free(vj);
 337    Free(Gwi);
 338    Free(mean);
 339    Free(var);
 340    Free(temp);
 341    if(lwork > 0) Free(work);
 342}
 343
 344
 345void gaussianPosterior_SelfInteraction(
 346	// IN/OUT
 347    double* sample /* v */,
 348    // OUTPUT
 349    double* posteriorMean, double* posteriorVar,
 350    // INPUT
 351    const int* option /*1:Sample only, 2:Sample&Mean&Var, 3:Sample&Mean&Var*/,
 352    const int* fromIndex, const int* toIndex, const double* obs /*o in the paper*/,
 353    const double* priorMean /*Gx_user*/,
 354    const double* obsVar, const double* priorVar /* var_v */,
 355    const int* nObs, const int* nLevelsThisEff, const int* nFactors,
 356    const int* nObsVar, const int *nPriorVar,
 357    const int* author_obsIndex, const int* author_oiStart, const int* author_oiNum,
 358    const int*  voter_obsIndex, const int*  voter_oiStart, const int*  voter_oiNum,
 359    // OTHER
 360    const int* debug
 361){
 362    double *sum_vv /*nFactors x nFactors*/, *sum_ov /*nFactors x 1*/, o, *vj /*nFactors x 1*/, *Gwi /*nFactors x 1*/,
 363           *work, size, *mean, *var, *temp, *rnd;
 364    int i,j,k,m, thisIndex, otherIndex, outputMeanVar, oIndex, lwork=-1, info, symCheck, role;
 365    char jobz = 'V', uplo = 'L';
 366
 367    symCheck = (*debug) - 2;
 368
 369    if(*option == 1){
 370        outputMeanVar = 0;
 371    }else if(*option == 2){
 372        outputMeanVar = 1;
 373    }else if(*option == 3){
 374        outputMeanVar = 1;
 375    }else error("Unknown option: %d", *option);
 376
 377    vj     = (double*)Calloc(*nFactors, double);
 378    Gwi    = (double*)Calloc(*nFactors, double);
 379    sum_ov = (double*)Calloc(*nFactors, double);
 380    sum_vv = (double*)Calloc((*nFactors)*(*nFactors), double);
 381    mean   = (double*)Calloc(*nFactors, double);
 382    var    = (double*)Calloc((*nFactors)*(*nFactors), double);
 383    temp   = (double*)Calloc((*nFactors)*(*nFactors), double);
 384
 385    GetRNGstate();
 386
 387    for(i=0; i<*nLevelsThisEff; i++){
 388
 389        thisIndex = i+1;
 390        for(j=0; j<*nFactors; j++) sum_ov[j] = 0;
 391        for(j=0; j<(*nFactors)*(*nFactors); j++) sum_vv[j] = 0;
 392        int nObs_for_i = 0;
 393
 394        // begin: compute sum_ov and sum_vv
 395        for(role=0; role<2; role++){
 396        	int *obsIndex = NULL, *oiNum = NULL, *oiStart=NULL, *thisEffIndex=NULL, *otherEffIndex=NULL;
 397        	if(role == 0){
 398        		// user i is an author
 399        		obsIndex = (int*)author_obsIndex;  oiNum = (int*)author_oiNum;  oiStart = (int*)author_oiStart;
 400        		thisEffIndex = (int*)toIndex;  otherEffIndex = (int*)fromIndex;
 401        	}else if(role == 1){
 402        		// user i is an voter
 403        		obsIndex = (int*)voter_obsIndex;   oiNum = (int*)voter_oiNum;   oiStart = (int*)voter_oiStart;
 404        		thisEffIndex = (int*)fromIndex; otherEffIndex = (int*)toIndex;
 405        	}else DIE_HERE;
 406
 407			for(j=0; j<oiNum[i]; j++){
 408
 409				nObs_for_i++;
 410
 411				oIndex = obsIndex[R_VEC(oiStart[i]+j)];
 412				otherIndex = otherEffIndex[R_VEC(oIndex)];
 413				if(otherIndex == i+1) error("self votes are not allowed");
 414
 415				if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
 416				if(*debug > 0) CHK_R_INDEX(otherIndex, *nLevelsThisEff);
 417				if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
 418
 419				o = obs[R_VEC(oIndex)];
 420				for(k=1; k<=*nFactors; k++) vj[R_VEC(k)] = sample[R_MAT(otherIndex,k,*nLevelsThisEff)];
 421
 422				double var_y_thisObs = 0;
 423				if((*nObsVar) == 1)            var_y_thisObs = obsVar[0];
 424				else if((*nObsVar) == (*nObs)) var_y_thisObs = obsVar[R_VEC(oIndex)];
 425				else error("nVar_y = %d, nObs = %d", *nObsVar, *nObs);
 426
 427				for(k=0; k<*nFactors; k++) sum_ov[k] += (o * vj[k]) / var_y_thisObs;
 428				for(k=0; k<*nFactors; k++)
 429					for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += (vj[k] * vj[m]) / var_y_thisObs;
 430			}
 431        }
 432        // end: compute sum_ov and sum_vv
 433
 434        for(k=1; k<=*nFactors; k++) Gwi[R_VEC(k)] = priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)];
 435
 436        if((*nPriorVar) == 1){
 437            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[0]);
 438            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[0];
 439        }else if((*nPriorVar) == (*nLevelsThisEff)*(*nFactors)*(*nFactors)){
 440
 441            for(k=0; k<*nFactors; k++)
 442                for(m=0; m<*nFactors; m++) temp[C_MAT(k,m,*nFactors)] = priorVar[C_3DA(i,k,m,*nLevelsThisEff,*nFactors)];
 443            sym_inv_byCholesky(temp, nFactors, &symCheck);
 444            // Now, temp is var_u[i]^-1
 445
 446            for(k=0; k<*nFactors; k++)
 447                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += temp[C_MAT(k,m,*nFactors)];
 448
 449            // reuse the space of vj
 450            for(k=0; k<*nFactors; k++) vj[k] = Gwi[k];
 451            for(k=0; k<*nFactors; k++){
 452                Gwi[k] = 0;
 453                for(m=0; m<*nFactors; m++) Gwi[k] += temp[C_MAT(k,m,*nFactors)] * vj[m];
 454            }
 455
 456        }else error("nVar_eff = %d, nThisEff = %d", *nPriorVar, *nLevelsThisEff);
 457
 458        // Now, sum_vv = var^-1
 459        //      Gwi = var_u[i]^-1 G wi
 460
 461        if((*nFactors) == 1){
 462            var[0]  = 1/sum_vv[0];
 463            mean[0] = var[0] * (sum_ov[0] + Gwi[0]);
 464            sample[i] = rnorm(mean[0], sqrt(var[0]));
 465            if(outputMeanVar){
 466                posteriorMean[i] = mean[0];
 467                posteriorVar[i]  = var[0];
 468            }
 469        }else{
 470            double *eigen_val, *eigen_vec;
 471            eigen_val = vj;
 472            eigen_vec = sum_vv;
 473
 474            if(*debug > 2) CHK_SYMMETRIC(eigen_vec, *nFactors);
 475
 476            //
 477            // Compute the variance-covariance matrix
 478            //
 479            // Allocate workspace for eigen value decomposition (only allocate once)
 480            if(lwork == -1){
 481                F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, &size, &lwork, &info);
 482                if(info != 0) error("error in dsyev(...)");
 483                lwork = (int)size;
 484                work  = (double*)Calloc(lwork, double);
 485            }
 486
 487            F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, work, &lwork, &info);
 488            if(info != 0) error("error in dsyev(...)");
 489            for(j=0; j<*nFactors; j++) eigen_val[j] = 1/eigen_val[j];
 490            // Now, eigen_val, eigen_vec are the eigen values and vectors of var
 491
 492            for(j=0; j<*nFactors; j++)
 493                for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 494            for(j=0; j<*nFactors; j++)
 495                for(k=0; k<*nFactors; k++){
 496                    var[C_MAT(j,k,*nFactors)] = 0;
 497                    for(m=0; m<*nFactors; m++) var[C_MAT(j,k,*nFactors)] += temp[C_MAT(j,m,*nFactors)] * eigen_vec[C_MAT(k,m,*nFactors)];
 498                }
 499            // Now, var is the variance-covariance matrix
 500
 501            //
 502            // Compute the mean vector
 503            //
 504
 505            // print_vector("sum_ov: ", sum_ov, *nFactors);
 506            // print_vector("Gwi: ", Gwi, *nFactors);
 507
 508            for(j=0; j<*nFactors; j++) sum_ov[j] += Gwi[j];
 509            for(j=0; j<*nFactors; j++){
 510                mean[j] = 0;
 511                for(k=0; k<*nFactors; k++) mean[j] += var[C_MAT(j,k,*nFactors)] * sum_ov[k];
 512            }
 513
 514            if(outputMeanVar){
 515                for(j=0; j<*nFactors; j++) posteriorMean[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 516                for(j=0; j<*nFactors; j++)
 517                    for(k=0; k<*nFactors; k++) posteriorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)] = var[C_MAT(j,k,*nFactors)];
 518            }
 519
 520            if(*debug >= 2 && nObs_for_i == 0){
 521                for(k=1; k<=*nFactors; k++){
 522                    CHK_SAME_NUMBER("mean[k] != fittedEff[k]", mean[R_VEC(k)], priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)]);
 523                }
 524                if((*nPriorVar) == 1){
 525                    for(j=0; j<*nFactors; j++)
 526                        for(k=0; k<*nFactors; k++){
 527                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[0]);}
 528                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 529                        }
 530                }else{
 531                    for(j=0; j<*nFactors; j++)
 532                        for(k=0; k<*nFactors; k++){
 533                            CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)]);
 534                        }
 535                }
 536            }
 537            //
 538            //  Generate the random vector
 539            //
 540			// reuse the space allocated for Gwi
 541			rnd = Gwi;
 542
 543			for(j=0; j<*nFactors; j++) eigen_val[j] = sqrt(eigen_val[j]);
 544			for(j=0; j<*nFactors; j++)
 545				for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 546
 547			for(j=0; j<*nFactors; j++) rnd[j] = norm_rand();
 548
 549			for(j=0; j<*nFactors; j++){
 550				sample[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 551				for(k=0; k<*nFactors; k++) sample[C_MAT(i,j,*nLevelsThisEff)] += temp[C_MAT(j,k,*nFactors)] * rnd[k];
 552			}
 553        }
 554    }
 555    PutRNGstate();
 556
 557    Free(sum_ov);  Free(sum_vv);  Free(vj);    Free(Gwi);
 558    Free(mean);    Free(var);     Free(temp);
 559    if(lwork > 0) Free(work);
 560}
 561
 562void gaussianPosterior_3WayInteraction(
 563    // OUTPUT
 564    double* sample,        /* nLevelsThisEff x nFactors */
 565    double* posteriorMean, /* nLevelsThisEff x nFactors */
 566    double* posteriorVar,  /* nLevelsThisEff x nFactors x nFactors */
 567    // INPUT
 568    const int* option /*1: output sample, 2: output mean & var, 3: output sample & Mean & Var*/,
 569    const int* thisEffIndex  /* nObs x 1 */,
 570    const int* otherEffIndex /* nObs x 1 */,
 571    const int* thirdEffIndex /* nObs x 1 */,
 572    const double* obs        /* nObs x 1 */,
 573    const double* priorMean   /* nLevelsThisEff  x nFactors */,
 574    const double* otherEff /* nLevelsOtherEff x nFactors */,
 575    const double* thirdEff /* nLevelsThirdEff x nFactors */,
 576    const double* obsVar   /* nObsVar x 1 */,
 577    const double* priorVar /* nPriorVar x 1 */,
 578    const int* nObs, const int* nLevelsThisEff, const int* nLevelsOtherEff, const int* nLevelsThirdEff,
 579    const int* nFactors, const int* nObsVar   /* 1 or nObs */,
 580    const int *nPriorVar /* 1 or nLevelsThisEff*nFactors*nFactors */,
 581    const int* obsIndex, const int* oiStart, const int* oiNum,
 582    // OTHER
 583    const int* debug
 584){
 585    double *sum_vv /*nFactors x nFactors*/, *sum_ov /*nFactors x 1*/, o, *vj /*nFactors x 1*/, *Gwi /*nFactors x 1*/,
 586           *work, size, *mean, *var, *temp, *rnd;
 587    int i,j,k,m, thisIndex, otherIndex, outputSample, outputMeanVar, oIndex, lwork=-1, info, symCheck;
 588    char jobz = 'V', uplo = 'L';
 589
 590    symCheck = (*debug) - 2;
 591
 592    if(*option == 1){
 593        outputSample = 1; outputMeanVar = 0;
 594    }else if(*option == 2){
 595        outputSample = 0; outputMeanVar = 1;
 596    }else if(*option == 3){
 597        outputSample = 1; outputMeanVar = 1;
 598    }else error("Unknown option: %d", *option);
 599
 600    vj     = (double*)Calloc(*nFactors, double);
 601    Gwi    = (double*)Calloc(*nFactors, double);
 602    sum_ov = (double*)Calloc(*nFactors, double);
 603    sum_vv = (double*)Calloc((*nFactors)*(*nFactors), double);
 604    mean   = (double*)Calloc(*nFactors, double);
 605    var    = (double*)Calloc((*nFactors)*(*nFactors), double);
 606    temp   = (double*)Calloc((*nFactors)*(*nFactors), double);
 607
 608    if(outputSample) GetRNGstate();
 609
 610    for(i=0; i<*nLevelsThisEff; i++){
 611
 612        thisIndex = i+1;
 613        for(j=0; j<*nFactors; j++) sum_ov[j] = 0;
 614        for(j=0; j<(*nFactors)*(*nFactors); j++) sum_vv[j] = 0;
 615
 616        for(j=0; j<oiNum[i]; j++){
 617
 618            oIndex = obsIndex[R_VEC(oiStart[i]+j)];
 619            otherIndex = otherEffIndex[R_VEC(oIndex)];
 620
 621            if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
 622            if(*debug > 0) CHK_R_INDEX(otherIndex, *nLevelsOtherEff);
 623            if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
 624
 625            o = obs[R_VEC(oIndex)];
 626
 627            if((*nLevelsThirdEff) > 0){
 628                int thirdIndex = thirdEffIndex[R_VEC(oIndex)];
 629                if(*debug > 0) CHK_R_INDEX(thirdIndex, *nLevelsThirdEff);
 630            	for(k=1; k<=*nFactors; k++)
 631            		vj[R_VEC(k)] = otherEff[R_MAT(otherIndex,k,*nLevelsOtherEff)] *
 632            		               thirdEff[R_MAT(thirdIndex,k,*nLevelsThirdEff)];
 633            }else{
 634            	for(k=1; k<=*nFactors; k++) vj[R_VEC(k)] = otherEff[R_MAT(otherIndex,k,*nLevelsOtherEff)];
 635            }
 636
 637            double var_y_thisObs = 0;
 638            if((*nObsVar) == 1)            var_y_thisObs = obsVar[0];
 639            else if((*nObsVar) == (*nObs)) var_y_thisObs = obsVar[R_VEC(oIndex)];
 640            else error("nVar_y = %d, nObs = %d", *nObsVar, *nObs);
 641
 642            for(k=0; k<*nFactors; k++) sum_ov[k] += (o * vj[k]) / var_y_thisObs;
 643            for(k=0; k<*nFactors; k++)
 644                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += (vj[k] * vj[m]) / var_y_thisObs;
 645        }
 646
 647        for(k=1; k<=*nFactors; k++) Gwi[R_VEC(k)] = priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)];
 648
 649        if((*nPriorVar) == 1){
 650            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[0]);
 651            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[0];
 652        }else if((*nPriorVar) == (*nFactors)){
 653            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[k]);
 654            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[k];
 655        }else if((*nPriorVar) == (*nLevelsThisEff)*(*nFactors)*(*nFactors)){
 656
 657            for(k=0; k<*nFactors; k++)
 658                for(m=0; m<*nFactors; m++) temp[C_MAT(k,m,*nFactors)] = priorVar[C_3DA(i,k,m,*nLevelsThisEff,*nFactors)];
 659            sym_inv_byCholesky(temp, nFactors, &symCheck);
 660            // Now, temp is var_u[i]^-1
 661
 662            for(k=0; k<*nFactors; k++)
 663                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += temp[C_MAT(k,m,*nFactors)];
 664
 665            // reuse the space of vj
 666            for(k=0; k<*nFactors; k++) vj[k] = Gwi[k];
 667            for(k=0; k<*nFactors; k++){
 668                Gwi[k] = 0;
 669                for(m=0; m<*nFactors; m++) Gwi[k] += temp[C_MAT(k,m,*nFactors)] * vj[m];
 670            }
 671
 672        }else error("nVar_eff = %d, nThisEff = %d", *nPriorVar, *nLevelsThisEff);
 673
 674        // Now, sum_vv = var^-1
 675        //      Gwi = var_u[i]^-1 G wi
 676
 677        if((*nFactors) == 1){
 678            var[0]  = 1/sum_vv[0];
 679            mean[0] = var[0] * (sum_ov[0] + Gwi[0]);
 680            if(outputSample){
 681                sample[i] = rnorm(mean[0], sqrt(var[0]));
 682            }
 683            if(outputMeanVar){
 684                posteriorMean[i] = mean[0];
 685                posteriorVar[i]  = var[0];
 686            }
 687        }else{
 688            double *eigen_val, *eigen_vec;
 689            eigen_val = vj;
 690            eigen_vec = sum_vv;
 691
 692            if(*debug > 2) CHK_SYMMETRIC(eigen_vec, *nFactors);
 693
 694            //
 695            // Compute the variance-covariance matrix
 696            //
 697            // Allocate workspace for eigen value decomposition (only allocate once)
 698            if(lwork == -1){
 699                F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, &size, &lwork, &info);
 700                if(info != 0) error("error in dsyev(...)");
 701                lwork = (int)size;
 702                work  = (double*)Calloc(lwork, double);
 703            }
 704
 705            F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, work, &lwork, &info);
 706            if(info != 0) error("error in dsyev(...)");
 707            for(j=0; j<*nFactors; j++) eigen_val[j] = 1/eigen_val[j];
 708            // Now, eigen_val, eigen_vec are the eigen values and vectors of var
 709
 710            for(j=0; j<*nFactors; j++)
 711                for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 712            for(j=0; j<*nFactors; j++)
 713                for(k=0; k<*nFactors; k++){
 714                    var[C_MAT(j,k,*nFactors)] = 0;
 715                    for(m=0; m<*nFactors; m++) var[C_MAT(j,k,*nFactors)] += temp[C_MAT(j,m,*nFactors)] * eigen_vec[C_MAT(k,m,*nFactors)];
 716                }
 717            // Now, var is the variance-covariance matrix
 718
 719            //
 720            // Compute the mean vector
 721            //
 722            for(j=0; j<*nFactors; j++) sum_ov[j] += Gwi[j];
 723            for(j=0; j<*nFactors; j++){
 724                mean[j] = 0;
 725                for(k=0; k<*nFactors; k++) mean[j] += var[C_MAT(j,k,*nFactors)] * sum_ov[k];
 726            }
 727
 728            if(outputMeanVar){
 729                for(j=0; j<*nFactors; j++) posteriorMean[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 730                for(j=0; j<*nFactors; j++)
 731                    for(k=0; k<*nFactors; k++) posteriorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)] = var[C_MAT(j,k,*nFactors)];
 732            }
 733
 734            if(*debug >= 2 && oiNum[i] == 0){
 735				// printf("i = %d\n",i);
 736				// print_vector("     mean: ", mean, *nFactors);
 737				// print_matrix("      var: ", var, *nFactors, *nFactors);
 738				// print_vector("prior_var: ", priorVar, *nPriorVar);
 739                for(k=1; k<=*nFactors; k++){
 740                    CHK_SAME_NUMBER("mean[k] != fittedEff[k]", mean[R_VEC(k)], priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)]);
 741                }
 742                if((*nPriorVar) == 1){
 743                    for(j=0; j<*nFactors; j++)
 744                        for(k=0; k<*nFactors; k++){
 745                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[0]);}
 746                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 747                        }
 748                }else if((*nPriorVar) == (*nFactors)){
 749                    for(j=0; j<*nFactors; j++)
 750                        for(k=0; k<*nFactors; k++){
 751                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[k]);}
 752                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 753                        }
 754                }else{
 755                    for(j=0; j<*nFactors; j++)
 756                        for(k=0; k<*nFactors; k++){
 757                            CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)]);
 758                        }
 759                }
 760            }
 761            //
 762            //  Generate the random vector
 763            //
 764            if(outputSample){
 765                // reuse the space allocated for Gwi
 766                rnd = Gwi;
 767
 768                for(j=0; j<*nFactors; j++) eigen_val[j] = sqrt(eigen_val[j]);
 769                for(j=0; j<*nFactors; j++)
 770                    for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 771
 772                for(j=0; j<*nFactors; j++) rnd[j] = norm_rand();
 773
 774                for(j=0; j<*nFactors; j++){
 775                    sample[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 776                    for(k=0; k<*nFactors; k++) sample[C_MAT(i,j,*nLevelsThisEff)] += temp[C_MAT(j,k,*nFactors)] * rnd[k];
 777                }
 778            }
 779        }
 780    }
 781    if(outputSample) PutRNGstate();
 782
 783    Free(sum_ov);
 784    Free(sum_vv);
 785    Free(vj);
 786    Free(Gwi);
 787    Free(mean);
 788    Free(var);
 789    Free(temp);
 790    if(lwork > 0) Free(work);
 791}
 792
 793void gaussianPosterior_SelfPlusOneInteraction(
 794	// IN/OUT
 795	double* sample /* the current sample: nLevelsThisEff x nFactors */,
 796	// OUTPUT
 797	double* posteriorMean, /* nLevelsThisEff x nFactors */
 798	double* posteriorVar,  /* nLevelsThisEff x nFactors x nFactors */
 799	// INPUT
 800	const int* option /*1: output sample, 2: output mean & var, 3: output sample & Mean & Var*/,
 801	const int* fromIndex /* nObs x 1 */,
 802	const int*   toIndex /* nObs x 1 */,
 803	const int* thirdEffIndex /* nObs x 1 */,
 804	const double* obs    /* nObs x 1 */,
 805	const double* priorMean /* nLevelsThisEff x nFactors */,
 806	const double* obsVar   /* nObsVar x 1 */,
 807	const double* priorVar /* nPriorVar x 1 */,
 808	const double* thirdEff /* nLevelsThirdEff x nFactors */,
 809	const int* nObs, const int* nLevelsThisEff, const int* nLevelsThirdEff, const int* nFactors,
 810	const int* nObsVar   /* 1 or nObs */,
 811	const int *nPriorVar /* 1 or nLevelsThisEff*nFactors*nFactors */,
 812	const int* from_obsIndex, const int* from_oiStart, const int* from_oiNum,
 813	const int*   to_obsIndex, const int*   to_oiStart, const int*   to_oiNum,
 814	// OTHER
 815	const int* debug
 816){
 817    double *sum_vv /*nFactors x nFactors*/, *sum_ov /*nFactors x 1*/, o, *vj /*nFactors x 1*/, *Gwi /*nFactors x 1*/,
 818           *work, size, *mean, *var, *temp, *rnd;
 819    int i,j,k,m, thisIndex, otherIndex, outputMeanVar, oIndex, lwork=-1, info, symCheck, role;
 820    char jobz = 'V', uplo = 'L';
 821
 822    symCheck = (*debug) - 2;
 823
 824    if(*option == 1){
 825        outputMeanVar = 0;
 826    }else if(*option == 2){
 827        outputMeanVar = 1;
 828    }else if(*option == 3){
 829        outputMeanVar = 1;
 830    }else error("Unknown option: %d", *option);
 831
 832    vj     = (double*)Calloc(*nFactors, double);
 833    Gwi    = (double*)Calloc(*nFactors, double);
 834    sum_ov = (double*)Calloc(*nFactors, double);
 835    sum_vv = (double*)Calloc((*nFactors)*(*nFactors), double);
 836    mean   = (double*)Calloc(*nFactors, double);
 837    var    = (double*)Calloc((*nFactors)*(*nFactors), double);
 838    temp   = (double*)Calloc((*nFactors)*(*nFactors), double);
 839
 840    GetRNGstate();
 841
 842    for(i=0; i<*nLevelsThisEff; i++){
 843
 844        thisIndex = i+1;
 845        for(j=0; j<*nFactors; j++) sum_ov[j] = 0;
 846        for(j=0; j<(*nFactors)*(*nFactors); j++) sum_vv[j] = 0;
 847        int nObs_for_i = 0;
 848
 849        // begin: compute sum_ov and sum_vv
 850        for(role=0; role<2; role++){
 851        	int *obsIndex = NULL, *oiNum = NULL, *oiStart=NULL, *thisEffIndex=NULL, *otherEffIndex=NULL;
 852        	if(role == 0){
 853        		// user i is an author
 854        		obsIndex = (int*)to_obsIndex;  oiNum = (int*)to_oiNum;  oiStart = (int*)to_oiStart;
 855        		thisEffIndex = (int*)toIndex;  otherEffIndex = (int*)fromIndex;
 856        	}else if(role == 1){
 857        		// user i is an voter
 858        		obsIndex = (int*)from_obsIndex; oiNum = (int*)from_oiNum;   oiStart = (int*)from_oiStart;
 859        		thisEffIndex = (int*)fromIndex; otherEffIndex = (int*)toIndex;
 860        	}else DIE_HERE;
 861
 862			for(j=0; j<oiNum[i]; j++){
 863
 864				nObs_for_i++;
 865
 866				oIndex = obsIndex[R_VEC(oiStart[i]+j)];
 867				otherIndex = otherEffIndex[R_VEC(oIndex)];
 868				if(otherIndex == i+1) error("self votes are not allowed");
 869
 870				if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
 871				if(*debug > 0) CHK_R_INDEX(otherIndex, *nLevelsThisEff);
 872				if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
 873
 874				o = obs[R_VEC(oIndex)];
 875
 876	            if((*nLevelsThirdEff) > 0){
 877	                int thirdIndex = thirdEffIndex[R_VEC(oIndex)];
 878	                if(*debug > 0) CHK_R_INDEX(thirdIndex, *nLevelsThirdEff);
 879	            	for(k=1; k<=*nFactors; k++)
 880	            		vj[R_VEC(k)] = sample[R_MAT(otherIndex,k,*nLevelsThisEff)] *
 881	            		               thirdEff[R_MAT(thirdIndex,k,*nLevelsThirdEff)];
 882	            }else{
 883	            	for(k=1; k<=*nFactors; k++) vj[R_VEC(k)] = sample[R_MAT(otherIndex,k,*nLevelsThisEff)];
 884	            }
 885
 886				double var_y_thisObs = 0;
 887				if((*nObsVar) == 1)            var_y_thisObs = obsVar[0];
 888				else if((*nObsVar) == (*nObs)) var_y_thisObs = obsVar[R_VEC(oIndex)];
 889				else error("nVar_y = %d, nObs = %d", *nObsVar, *nObs);
 890
 891				for(k=0; k<*nFactors; k++) sum_ov[k] += (o * vj[k]) / var_y_thisObs;
 892				for(k=0; k<*nFactors; k++)
 893					for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += (vj[k] * vj[m]) / var_y_thisObs;
 894			}
 895        }
 896        // end: compute sum_ov and sum_vv
 897
 898        for(k=1; k<=*nFactors; k++) Gwi[R_VEC(k)] = priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)];
 899
 900        if((*nPriorVar) == 1){
 901            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[0]);
 902            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[0];
 903        }else if((*nPriorVar) == (*nFactors)){
 904            for(k=0; k<*nFactors; k++) sum_vv[C_MAT(k,k,*nFactors)] += (1/priorVar[k]);
 905            for(k=0; k<*nFactors; k++) Gwi[k] /= priorVar[k];
 906        }else if((*nPriorVar) == (*nLevelsThisEff)*(*nFactors)*(*nFactors)){
 907
 908            for(k=0; k<*nFactors; k++)
 909                for(m=0; m<*nFactors; m++) temp[C_MAT(k,m,*nFactors)] = priorVar[C_3DA(i,k,m,*nLevelsThisEff,*nFactors)];
 910            sym_inv_byCholesky(temp, nFactors, &symCheck);
 911            // Now, temp is var_u[i]^-1
 912
 913            for(k=0; k<*nFactors; k++)
 914                for(m=0; m<*nFactors; m++) sum_vv[C_MAT(k,m,*nFactors)] += temp[C_MAT(k,m,*nFactors)];
 915
 916            // reuse the space of vj
 917            for(k=0; k<*nFactors; k++) vj[k] = Gwi[k];
 918            for(k=0; k<*nFactors; k++){
 919                Gwi[k] = 0;
 920                for(m=0; m<*nFactors; m++) Gwi[k] += temp[C_MAT(k,m,*nFactors)] * vj[m];
 921            }
 922
 923        }else error("nVar_eff = %d, nThisEff = %d", *nPriorVar, *nLevelsThisEff);
 924
 925        // Now, sum_vv = var^-1
 926        //      Gwi = var_u[i]^-1 G wi
 927
 928        if((*nFactors) == 1){
 929            var[0]  = 1/sum_vv[0];
 930            mean[0] = var[0] * (sum_ov[0] + Gwi[0]);
 931            sample[i] = rnorm(mean[0], sqrt(var[0]));
 932            if(outputMeanVar){
 933                posteriorMean[i] = mean[0];
 934                posteriorVar[i]  = var[0];
 935            }
 936        }else{
 937            double *eigen_val, *eigen_vec;
 938            eigen_val = vj;
 939            eigen_vec = sum_vv;
 940
 941            if(*debug > 2) CHK_SYMMETRIC(eigen_vec, *nFactors);
 942
 943            //
 944            // Compute the variance-covariance matrix
 945            //
 946            // Allocate workspace for eigen value decomposition (only allocate once)
 947            if(lwork == -1){
 948                F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, &size, &lwork, &info);
 949                if(info != 0) error("error in dsyev(...)");
 950                lwork = (int)size;
 951                work  = (double*)Calloc(lwork, double);
 952            }
 953
 954            F77_NAME(dsyev)(&jobz, &uplo, nFactors, eigen_vec, nFactors, eigen_val, work, &lwork, &info);
 955            if(info != 0) error("error in dsyev(...)");
 956            for(j=0; j<*nFactors; j++) eigen_val[j] = 1/eigen_val[j];
 957            // Now, eigen_val, eigen_vec are the eigen values and vectors of var
 958
 959            for(j=0; j<*nFactors; j++)
 960                for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
 961            for(j=0; j<*nFactors; j++)
 962                for(k=0; k<*nFactors; k++){
 963                    var[C_MAT(j,k,*nFactors)] = 0;
 964                    for(m=0; m<*nFactors; m++) var[C_MAT(j,k,*nFactors)] += temp[C_MAT(j,m,*nFactors)] * eigen_vec[C_MAT(k,m,*nFactors)];
 965                }
 966            // Now, var is the variance-covariance matrix
 967
 968            //
 969            // Compute the mean vector
 970            //
 971            for(j=0; j<*nFactors; j++) sum_ov[j] += Gwi[j];
 972            for(j=0; j<*nFactors; j++){
 973                mean[j] = 0;
 974                for(k=0; k<*nFactors; k++) mean[j] += var[C_MAT(j,k,*nFactors)] * sum_ov[k];
 975            }
 976
 977            if(outputMeanVar){
 978                for(j=0; j<*nFactors; j++) posteriorMean[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
 979                for(j=0; j<*nFactors; j++)
 980                    for(k=0; k<*nFactors; k++) posteriorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)] = var[C_MAT(j,k,*nFactors)];
 981            }
 982
 983            if(*debug >= 2 && nObs_for_i == 0){
 984                for(k=1; k<=*nFactors; k++){
 985                    CHK_SAME_NUMBER("mean[k] != fittedEff[k]", mean[R_VEC(k)], priorMean[R_MAT(thisIndex,k,*nLevelsThisEff)]);
 986                }
 987                if((*nPriorVar) == 1){
 988                    for(j=0; j<*nFactors; j++)
 989                        for(k=0; k<*nFactors; k++){
 990                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[0]);}
 991                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 992                        }
 993                }else if((*nPriorVar) == (*nFactors)){
 994                    for(j=0; j<*nFactors; j++)
 995                        for(k=0; k<*nFactors; k++){
 996                            if(j==k){ CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[k]);}
 997                            else{     CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], 0);}
 998                        }
 999                }else{
1000                    for(j=0; j<*nFactors; j++)
1001                        for(k=0; k<*nFactors; k++){
1002                            CHK_SAME_NUMBER("var != var_eff", var[C_MAT(j,k,*nFactors)], priorVar[C_3DA(i,j,k,*nLevelsThisEff,*nFactors)]);
1003                        }
1004                }
1005            }
1006            //
1007            //  Generate the random vector
1008            //
1009			// reuse the space allocated for Gwi
1010			rnd = Gwi;
1011
1012			for(j=0; j<*nFactors; j++) eigen_val[j] = sqrt(eigen_val[j]);
1013			for(j=0; j<*nFactors; j++)
1014				for(k=0; k<*nFactors; k++) temp[C_MAT(j,k,*nFactors)] = eigen_vec[C_MAT(j,k,*nFactors)] * eigen_val[k];
1015
1016			for(j=0; j<*nFactors; j++) rnd[j] = norm_rand();
1017
1018			for(j=0; j<*nFactors; j++){
1019				sample[C_MAT(i,j,*nLevelsThisEff)] = mean[j];
1020				for(k=0; k<*nFactors; k++) sample[C_MAT(i,j,*nLevelsThisEff)] += temp[C_MAT(j,k,*nFactors)] * rnd[k];
1021			}
1022        }
1023    }
1024    PutRNGstate();
1025
1026    Free(sum_ov);  Free(sum_vv);  Free(vj);    Free(Gwi);
1027    Free(mean);    Free(var);     Free(temp);
1028    if(lwork > 0) Free(work);
1029}
1030
1031/**
1032 * contextSample can be set the same as contextPosteriorMean
1033 * globalSample  can be set the same as globalPosteriorMean
1034 * If so, the output will be the sample only.
1035 */
1036void gaussianPosterior_mainEffect_2Levels(
1037	// OUTPUT
1038	double* contextSample /* nLevelsThisEff x nContexts */,
1039	double* globalSample  /* nLevelsThisEff x 1 */,
1040	double* contextPosteriorMean /* nLevelsThisEff x nContexts: This is the posterior given the globalSample */,
1041	double* contextPosteriorVar  /* nLevelsThisEff x nContexts: This is the posterior given the globalSample */,
1042	double* globalPosteriorMean /* nLevelsThisEff x 1 */,
1043	double* globalPosteriorVar  /* nLevelsThisEff x 1 */,
1044	//INPUT
1045	const int* thisEffIndex /* nObs x 1 */,
1046	const int* context /* nObs x 1 */,
1047	const double* obs /* nObs x 1 */,
1048	const double* q /* nContext x 1 */,
1049	const double* contextOffset, /* nLevelsThisEff x nContexts */
1050	const double* obsVar   /* nObsVar x 1 */,
1051	const double* contextPriorVar /* see nContextPriorVar */,
1052	const double* globalPriorVar /*  see nGlobalPriorVar */,
1053	const int* numObs, const int* numLevelsThisEff, const int* numContexts,
1054	const int* numObsVar   /* 1 or nObs */,
1055	const int* numContextPriorVar /* nContexts or nLevelsThisEff*nContexts */,
1056	const int* numGlobalPriorVar  /* 1 or nLevelsThisEff */,
1057	const int* numContextOffset /* 0 or nLevelsThisEff or nLevelsThisEff*nContexts */,
1058	// OTHER
1059	const int* debug, const int* verbose
1060){
1061
1062	const int nObs = (*numObs), nItems = (*numLevelsThisEff), nCategories = (*numContexts),
1063		      nObsVar = (*numObsVar), nContextPriorVar = (*numContextPriorVar), nGlobalPriorVar = (*numGlobalPriorVar),
1064		      nCatOffset = (*numContextOffset);
1065	const int *itemIndex = thisEffIndex, *categoryIndex=context;
1066	double *b_mean=contextPosteriorMean, *b_var=contextPosteriorVar,
1067		   *a_mean=globalPosteriorMean,  *a_var=globalPosteriorVar;
1068	const double *b_offset=contextOffset, *var_obs=obsVar, *var_b=contextPriorVar, *var_a=globalPriorVar;
1069	// Model:
1070	// 		obs[n] ~ N(mean=b[i,k],                    var=var_obs[n])
1071	//      b[i,k] ~ N(mean=b_offset[i,k] + q[k]*a[i], var=var_b[i,k])
1072	//      a[i]   ~ N(mean=0,                         var=var_a[i])
1073
1074	if(nObsVar != 1 && nObsVar != nObs) error("nObsVar != 1 && nObsVar != nObs");
1075	if(nContextPriorVar != nCategories && nContextPriorVar != nItems*nCategories) error("nContextPriorVar != nCategories && nContextPriorVar != nItem*nCategories");
1076	if(nGlobalPriorVar != 1 && nGlobalPriorVar != nItems) error("nContextPriorVar != nCategories && nContextPriorVar != nItem*nCategories");
1077
1078	// Compute sufficient statistics
1079	// 		b_mean[i,k] = sum_j { o_{ijk}/var_obs_{ijk} } = C_ik
1080	//       b_var[i,k] = sum_j {       1/var_obs_{ijk} } = F_ik
1081	for(int n=0; n < nItems*nCategories; n++){ b_mean[n] = 0; b_var[n] = 0; }
1082	for(int j=0; j<nObs; j++){
1083		int i = itemIndex[j] - 1;     CHK_C_INDEX(i,nItems);
1084		int k = categoryIndex[j] - 1; CHK_C_INDEX(k,nCategories);
1085		int ik = C_MAT(i,k,nItems);
1086		double var = (nObsVar==1 ? var_obs[0] : var_obs[j]);
1087
1088		double o = obs[j];
1089		if(nCatOffset == nItems)                  o -= b_offset[i];
1090		else if(nCatOffset == nItems*nCategories) o -= b_offset[ik];
1091		else if(nCatOffset != 0) error("nCatOffset = %d", nCatOffset);
1092
1093		b_mean[ik] += o / var;
1094		b_var[ik]  += 1 / var;
1095	}
1096
1097	GetRNGstate();
1098
1099	// Compute E[a[i] | obs] and Var[a[i] | obs]
1100	for(int i=0; i<nItems; i++){
1101		double sum_for_mean = 0; // sum_k E[a_i | obs_ik]/Var[a_i | obs_ik]
1102		double sum_for_var  = 0; // sum_k (1/Var[a_i | obs_ik] - 1/var_a)
1103		for(int k=0; k<nCategories; k++){
1104			double var_b_this = (nContextPriorVar==nCategories ? var_b[k] : var_b[C_MAT(i,k,nItems)]);
1105			double F = b_var[C_MAT(i,k,nItems)];
1106			if(F == 0) continue;
1107			double C = b_mean[C_MAT(i,k,nItems)];
1108			if(var_b_this > 0){
1109				double A = (1/var_b_this) + F;
1110				// E[a_i | obs_ik]/Var[a_i | obs_ik] = (C * q[k]) / (A * var_b[k])
1111				// (1/Var[a_i | obs_ik] - 1/var_a)   = (q[k]^2 / var_b[k]) * (1 - 1/(A * var_b[k]))
1112				sum_for_mean += (C * q[k]) / (A * var_b_this);
1113				sum_for_var  += (q[k]*q[k] / var_b_this) * (1 - 1/(A * var_b_this));
1114			}else if(var_b_this == 0){
1115				STOP("not yet support prior var = 0");
1116			}else STOP("var_b_this < 0");
1117		}
1118		if((*verbose) >= 100) printf("   i=%d:  sum.for.a.mean=%f, sum.for.a.var=%f\n", i+1, sum_for_mean, sum_for_var);
1119
1120		// Compute E[a[i] | obs] and Var[a[i] | obs]
1121		double var_a_this = (nGlobalPriorVar==1 ? var_a[0] : var_a[i]);
1122		a_var[i]  = 1/( (1/var_a_this) + sum_for_var );
1123		a_mean[i] = a_var[i] * sum_for_mean;
1124
1125		// draw for the globalSample
1126		globalSample[i] = rnorm(a_mean[i], sqrt(a_var[i]));
1127	}
1128
1129	// Compute E[b[i,k] | a, obs] and Var[b[i,k] | a, obs]
1130	for(int i=0; i<nItems; i++){ for(int k=0; k<nCategories; k++){
1131		int ik = C_MAT(i,k,nItems);
1132		double var_b_this = (nContextPriorVar==nCategories ? var_b[k] : var_b[C_MAT(i,k,nItems)]);
1133		double A = 1 / ( (1/var_b_this) + b_var[ik] );
1134		double D = (A * q[k]) / var_b_this;
1135
1136		b_mean[ik] = D * globalSample[i] + A * b_mean[ik];
1137		b_var[ik]  = A;
1138
1139		if(nCatOffset == nItems)                  b_mean[ik] += b_offset[i];
1140		else if(nCatOffset == nItems*nCategories) b_mean[ik] += b_offset[ik];
1141		else if(nCatOffset != 0) error("nCatOffset = %d", nCatOffset);
1142
1143		contextSample[ik] = rnorm(b_mean[ik], sqrt(b_var[ik]));
1144	}}
1145
1146    PutRNGstate();
1147}
1148