Skip to content
Browse files

RLFM logistic

  • Loading branch information...
1 parent 167c55f commit 622fa896796f0cad8025fda08d015739ae8b8f28 Liang Zhang committed Feb 29, 2012
View
406 src/RLFM-ars-logistic/C/MCEM_EStep_logistic.c
@@ -0,0 +1,406 @@
+/*
+ Copyright (c) 2011, Yahoo! Inc. All rights reserved.
+ Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
+
+ Author: Liang Zhang
+*/
+
+/*
+To Compile:
+R CMD SHLIB util.c MCEM_EStep.c -o MCEM_EStep.so
+
+*/
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include "util.h"
+//#include "logistic.h"
+#include "arsspline.h"
+
+///// centered logistic ars code starts from here
+// FUNCTION: mainEffect_condMeanVar
+// thisEffIndex: nObs x 1 otherEffIndex: nObs x 1 y_minus_xb_minus_uv: nObs x 1
+// fittedEff: nThisEff x 1 otherEff: nOtherEff x 1
+// For alpha, call
+// mainEffect_condMeanVar(1, user, item, y-xb-uv, g0w, beta, var_y, var_alpha, nObs, nUsers, nItems, sample, NULL, NULL, ...);
+// For beta, call
+// mainEffect_condMeanVar(1, item, user, y-xb-uv, d0z, alpha, var_y, var_beta, nObs, nItems, nUsers, sample, NULL, NULL, ...);
+// Observation index (consider, say, user i) - R indices (starting from 1, NOT 0)
+// obsIndex[ oiStart[i]+0 ], ..., obsIndex[ oiStart[i]+oiNum[i]-1 ] are the indices of user i's observations
+// in y, x, user, item
+void mainEffect_condMeanVarSample_arscid(
+ // OUTPUT
+ double* outSample, //int* naccepts,
+ double* ars_XI, /* ars_ninit X ncov X nThisEff*/
+ //INPUT
+ const int* thisEffIndex /*user or item*/, const int* otherEffIndex /*item or user*/,
+ const double* y, const double* offset, /*xb+uv*/
+ const double* fittedEff /*g0w or d0z*/,
+ const double* otherEff /*beta or alpha*/,
+ const double* var_eff /*var_alpha or var_beta*/,
+ const int* nObs, const int* nThisEff, const int* nOtherEff,
+ const int* obsIndex, const int* oiStart, const int* oiNum,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug
+
+ ){
+ int i,j,k,m, thisIndex, otherIndex, oIndex;
+ //*naccepts = 0;
+ GetRNGstate();
+ int neval = 0;
+ for(i=0; i<*nThisEff; i++){
+ thisIndex = i+1;
+ if(oiNum[i]>0){
+ double* y_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the observations for this user
+ double* offset_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the offsets for this user
+ double* X = (double*)calloc(oiNum[i], sizeof(double)); // the design matrix nobs*ncov, ncov=1 for main effects
+ for (j=0;j<oiNum[i];j++)
+ X[j] = 1;
+ for(j=0; j<oiNum[i]; j++){
+ oIndex = obsIndex[R_VEC(oiStart[i]+j)];
+ otherIndex = otherEffIndex[R_VEC(oIndex)];
+ if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
+ if(*debug > 0) CHK_R_INDEX(otherIndex, *nOtherEff);
+ if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
+ y_thisEff[j] = y[R_VEC(oIndex)];
+ offset_thisEff[j] = offset[R_VEC(oIndex)] + otherEff[R_VEC(otherIndex)];
+ }
+ //double new_thisEff = 0;
+ //int accept = 0;
+ int ncov = 1;
+ //MHlogistic(&outSample[i], X, y_thisEff, offset_thisEff, &fittedEff[i], var_eff, &oiNum[i], &ncov, &new_thisEff, &accept);
+ ARSLOGISTICSPLINE(offset_thisEff, &fittedEff[i], var_eff, X, y_thisEff, &ncov, &oiNum[i], &outSample[i], ars_qcent, ars_ninit, ars_ninit, ars_xl, ars_xu, &ars_XI[i*(*ars_ninit)], ars_alpha, &neval);
+ //printf("%f\n",outSample[i]);
+ //outSample[i] = new_thisEff;
+ //if (accept==1) *naccepts = *naccepts + 1;
+ free(y_thisEff);
+ free(offset_thisEff);
+ free(X);
+ }
+ }
+ PutRNGstate();
+}
+void factor_condMeanVarSample_arscid(
+ // OUTPUT
+ double* outSample, //int* naccepts,
+ double* ars_XI, /* ars_ninit X ncov X nThisEff*/
+ // INPUT
+ const int* thisEffIndex /*user or item*/, const int* otherEffIndex /*item or user*/,
+ const double* y, const double* offset, /*xb+alpha+beta*/
+ const double* fittedEff /*Gw or Dz*/,
+ const double* otherEff /*v or u*/,
+ const double* var_eff /*var_u or var_v*/,
+ const int* nObs, const int* nThisEff, const int* nOtherEff, const int* nFactors,
+ const int* obsIndex, const int* oiStart, const int* oiNum,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug
+ ){
+ int i,j,k,m, thisIndex, otherIndex, oIndex;
+ int* neval = (int*)calloc(*nFactors,sizeof(double));
+ //*naccepts = 0;
+ GetRNGstate();
+ for(i=0; i<*nThisEff; i++){
+ thisIndex = i+1;
+ if(oiNum[i]>0)
+ {
+ double* y_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the observations for this user
+ double* offset_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the offsets for this user
+ double* vj = (double*)calloc(oiNum[i]*(*nFactors), sizeof(double)); // the design matrix nobs*ncov, ncov=nfactors
+ double* thisEff_i = (double*)calloc((*nFactors),sizeof(double));
+ //double* new_thisEff_i = (double*)calloc((*nFactors),sizeof(double));
+ double* fittedEff_i = (double*)calloc((*nFactors),sizeof(double));
+ double* var_this_eff = (double*)calloc((*nFactors),sizeof(double));
+
+ for (j=0;j<*nFactors;j++)
+ {
+ thisEff_i[j] = outSample[C_MAT(i,j,*nThisEff)];
+ fittedEff_i[j] = fittedEff[C_MAT(i,j,*nThisEff)];
+ var_this_eff[j] = var_eff[j];
+ }
+ for(j=0; j<oiNum[i]; j++){
+ oIndex = obsIndex[R_VEC(oiStart[i]+j)];
+ otherIndex = otherEffIndex[R_VEC(oIndex)];
+ if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
+ if(*debug > 0) CHK_R_INDEX(otherIndex, *nOtherEff);
+ if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
+ y_thisEff[j] = y[R_VEC(oIndex)];
+ offset_thisEff[j] = offset[R_VEC(oIndex)];
+ for(k=1; k<=*nFactors; k++) {
+ vj[R_MAT((j+1),k,oiNum[i])] = otherEff[R_MAT(otherIndex,k,*nOtherEff)];
+ }
+ }
+ //double new_thisEff = 0;
+ // int accept = 0;
+ // MHlogistic(thisEff_i, vj, y_thisEff, offset_thisEff, fittedEff_i, var_eff, &oiNum[i], nFactors, new_thisEff_i, &accept);
+ //for (j=0;j<*nFactors;j++)
+ // printf("%f ",thisEff_i[j]);
+ //printf("\n");
+ ARSLOGISTICSPLINE(offset_thisEff, fittedEff_i, var_this_eff, vj, y_thisEff, nFactors, &oiNum[i], thisEff_i, ars_qcent, ars_ninit, ars_ninit, ars_xl, ars_xu, &ars_XI[i*(*ars_ninit)*(*nFactors)], ars_alpha, neval);
+
+ for (j=0;j<*nFactors;j++)
+ {
+ //printf("%f ",thisEff_i[j]);
+ outSample[C_MAT(i,j,*nThisEff)] = thisEff_i[j];
+ //if (accept==1) *naccepts = *naccepts + 1;
+ }
+ //printf("\n");
+ free(y_thisEff);
+ free(offset_thisEff);
+ free(vj);
+ free(thisEff_i);
+ //free(new_thisEff_i);
+ free(fittedEff_i);
+ free(var_this_eff);
+ }
+ }
+ PutRNGstate();
+ free(neval);
+}
+
+
+void computeMeanSumvar_arsc(
+ // OUTPUT
+ double *mean, double *sumvar,
+ // INPUT
+ const double *sum, const double *sos /* sum of squares */, const int length, const int nSamples
+ ){
+ int k;
+ double ratio = ((double)nSamples) / (nSamples - 1.0);
+ sumvar[0] = 0;
+ for(k=0; k<length; k++){
+ mean[k] = sum[k] / nSamples;
+ sumvar[0] += (sos[k] / (nSamples - 1.0)) - (ratio * mean[k] * mean[k]);
+ }
+}
+
+void computeMeanSumvarFactor_arsc(
+ // OUTPUT
+ double *mean, double *sumvar /* nFactors x 1*/,
+ // INPUT
+ const double *sum, const double *sos /* sum of squares */, const int length, const int nFactors, const int nSamples
+ ){
+ int k,l;
+ double ratio = ((double)nSamples) / (nSamples - 1.0);
+ for (k=0;k<nFactors;k++)
+ sumvar[k] = 0;
+ for (l=0; l<nFactors;l++)
+ for(k=0; k<length; k++){
+ mean[C_MAT(k,l,length)] = sum[C_MAT(k,l,length)] / nSamples;
+ sumvar[l] += (sos[C_MAT(k,l,length)] / (nSamples - 1.0)) - (ratio * mean[C_MAT(k,l,length)] * mean[C_MAT(k,l,length)]);
+ }
+}
+
+void InitializeXI_arsc(double* xi, const int ninit, const int ncov, const int nthisEff, const double xu, const double xl)
+{
+ int i,j,k;
+ for (i=0;i<nthisEff;i++)
+ for (j=0;j<ncov;j++)
+ for (k=0;k<ninit;k++)
+ {
+ double tmp = xl + (k + 2.0)*(xu - xl)/(ninit + 1.0);
+ if (tmp>=xu) tmp = tmp - .1;
+ if (tmp<=xl) tmp = tmp + .1;
+ xi[i*(ncov*ninit)+j*ninit+k] = tmp;
+ }
+}
+void MCEM_EStep_arscid(
+ // OUTPUT
+ double* o_mean /*nObs x 1*/,
+ double* alpha_mean/*nUsers x 1*/, double* alpha_sumvar/*1x1*/,
+ double* beta_mean/*nItems x 1*/, double* beta_sumvar/*1x1*/,
+ double* u_mean/*nUsers x nFactors*/, double* u_sumvar/*nFactors x 1*/,
+ double* v_mean/*nItems x nFactors*/, double* v_sumvar/*nFactors x 1*/,
+ double* ars_XI_alpha, double* ars_XI_beta,
+ double* ars_XI_u, double* ars_XI_v,
+ //double* acceptrate_maineff, // the overall acceptance rate for main effects
+ //double* acceptrate_fact, // the overall acceptance rate for factors
+ // INPUT
+ const int* nSamples, const int* nBurnin,
+ const int* user/*nObs x 1*/, const int* item/*nObs x 1*/,
+ const double* y/*nObs x 1*/, const double* xb/*nObs x 1*/,
+ const double* g0w/*nUsers x 1*/, const double* d0z/*nItems x 1*/,
+ const double* Gw/*nUsers x nFactors*/, const double* Dz/*nItems x nFactors*/,
+ const double* alpha_in/*nUsers x 1*/, const double* beta_in/*nItems x 1*/,
+ const double* u_in/*nUsers x nFactors*/, const double* v_in/*nItems x nFactors*/,
+ const double* var_alpha, const double* var_beta,
+ const double* var_u, const double* var_v, /*nFactors x 1*/
+ const int* nObs, const int* nUsers, const int* nItems, const int* nFactors,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug, const int* main_effects, const int* beta_int, const int* center
+
+ ){
+ int *obsIndex_user, *oiStart_user, *oiNum_user,
+ *obsIndex_item, *oiStart_item, *oiNum_item,
+ s, k, f, user_i, item_j, naccepts;
+ double *alpha, *beta, *u, *v,
+ *alpha_sum, *alpha_sos, *beta_sum, *beta_sos, *u_sum, *u_sos, *v_sum, *v_sos,
+ *temp, *xb_plus_uv, *xb_plus_alpha_beta, uv;
+ double *ars_xl_v = (double*)calloc(1, sizeof(double));
+ ars_xl_v[0] = 0;
+
+ // double *ars_XI_alpha, *ars_XI_beta, *ars_XI_u, *ars_XI_v;
+ //double accept_maineff_denom = 0, accept_maineff_numeri = 0,
+ // accept_fact_denom = 0, accept_fact_numeri = 0;
+
+
+ // Allocate space for sum and sum-of-squares
+ alpha_sum = (double*)calloc(*nUsers, sizeof(double));
+ beta_sum = (double*)calloc(*nItems, sizeof(double));
+ u_sum = (double*)calloc((*nUsers)*(*nFactors), sizeof(double));
+ v_sum = (double*)calloc((*nItems)*(*nFactors), sizeof(double));
+ alpha_sos = (double*)calloc(*nUsers, sizeof(double));
+ beta_sos = (double*)calloc(*nItems, sizeof(double));
+ u_sos = (double*)calloc((*nUsers)*(*nFactors), sizeof(double));
+ v_sos = (double*)calloc((*nItems)*(*nFactors), sizeof(double));
+ // Allocate space for the observation indices
+ obsIndex_user = (int*)calloc(*nObs, sizeof(int));
+ oiStart_user = (int*)calloc(*nUsers, sizeof(int));
+ oiNum_user = (int*)calloc(*nUsers, sizeof(int));
+ obsIndex_item = (int*)calloc(*nObs, sizeof(int));
+ oiStart_item = (int*)calloc(*nItems, sizeof(int));
+ oiNum_item = (int*)calloc(*nItems, sizeof(int));
+ // Allocate space for XI ars
+ //ars_XI_alpha = (double*)calloc((*ars_ninit)*(*nUsers), sizeof(double));
+ //ars_XI_beta = (double*)calloc((*ars_ninit)*(*nItems), sizeof(double));
+ //ars_XI_u = (double*)calloc((*ars_ninit)*(*nFactors)*(*nUsers), sizeof(double));
+ //ars_XI_v = (double*)calloc((*ars_ninit)*(*nFactors)*(*nItems), sizeof(double));
+
+ // Allocate temp space
+ temp = (double*)calloc(*nObs, sizeof(double));
+
+ memset(o_mean,0,(*nObs)*sizeof(double));
+ // Use the memory space of the output to store the current alpha, beta, u and v
+ alpha = alpha_mean; beta = beta_mean; u = u_mean; v = v_mean;
+ // Use the temp space for both xb_plus_uv and xb_plus_alpha_beta
+ xb_plus_uv = temp;
+ xb_plus_alpha_beta = temp;
+
+
+ // Create Observation indices for users and items
+ generateObsIndex(obsIndex_user, oiStart_user, oiNum_user, user, nObs, nUsers, debug);
+ generateObsIndex(obsIndex_item, oiStart_item, oiNum_item, item, nObs, nItems, debug);
+
+ // Initialize alpha, beta, u, v
+ for(k=0; k<*nUsers; k++) alpha[k] = alpha_in[k];
+ for(k=0; k<*nItems; k++) beta[k] = beta_in[k];
+ for(k=0; k<(*nUsers)*(*nFactors); k++) u[k] = u_in[k];
+ for(k=0; k<(*nItems)*(*nFactors); k++) v[k] = v_in[k];
+
+ // Initialize ars_XI_alpha, etc.
+ InitializeXI_arsc(ars_XI_alpha, *ars_ninit, 1, *nUsers, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_beta, *ars_ninit, 1, *nItems, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_u, *ars_ninit, *nFactors, *nUsers, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_v, *ars_ninit, *nFactors, *nItems, *ars_xu, *ars_xl_v);
+
+ for(s=0; s<(*nSamples+*nBurnin); s++){
+
+ // Compute xb+uv
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ uv = 0;
+ for(f=1; f<=*nFactors; f++) uv += u[R_MAT(user_i,f,*nUsers)] * v[R_MAT(item_j,f,*nItems)];
+ xb_plus_uv[k] = xb[k] + uv;
+ }
+ // Sample alpha
+ //mainEffect_condMeanVarSample_arsc(alpha, &naccepts, user, item, y, xb_plus_uv, g0w, beta, var_alpha, nObs, nUsers, nItems, obsIndex_user, oiStart_user, oiNum_user, debug);
+ //accept_maineff_denom += *nUsers; accept_maineff_numeri += naccepts;
+ mainEffect_condMeanVarSample_arscid(alpha, ars_XI_alpha, user, item, y, xb_plus_uv, g0w, beta, var_alpha, nObs, nUsers, nItems, obsIndex_user, oiStart_user, oiNum_user, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ //center
+ if(*center==1){ center_array(alpha, *nUsers); }
+ // Sample beta
+ //mainEffect_condMeanVarSample_arsc(beta, &naccepts, item, user, y, xb_plus_uv, d0z, alpha, var_beta, nObs, nItems, nUsers, obsIndex_item, oiStart_item, oiNum_item, debug);
+ //accept_maineff_denom += *nItems; accept_maineff_numeri += naccepts;
+ mainEffect_condMeanVarSample_arscid(beta, ars_XI_beta, item, user, y, xb_plus_uv, d0z, alpha, var_beta, nObs, nItems, nUsers, obsIndex_item, oiStart_item, oiNum_item, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ //subtract mean from the betas ...
+ if(*center==1 && *beta_int==0){ center_array(beta, *nItems); }
+
+ // Compute y - (xb + alpha + beta)
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ xb_plus_alpha_beta[k] = xb[k] + alpha[R_VEC(user_i)] + beta[R_VEC(item_j)];
+ }
+
+ if(*main_effects==0){
+ // Sample u
+ //factor_condMeanVarSample_arsc(u, &naccepts, user, item, y, xb_plus_alpha_beta, Gw, v, var_u, nObs, nUsers, nItems, nFactors, obsIndex_user, oiStart_user, oiNum_user, debug);
+ //accept_fact_denom += *nUsers; accept_fact_numeri += naccepts;
+ factor_condMeanVarSample_arscid(u, ars_XI_u, user, item, y, xb_plus_alpha_beta, Gw, v, var_u, nObs, nUsers, nItems, nFactors, obsIndex_user, oiStart_user, oiNum_user, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ if(*center==1){ center_array_2d(u, *nUsers, *nFactors, 2); }
+ // Sample v
+ //factor_condMeanVarSample_arsc(v, &naccepts, item, user, y, xb_plus_alpha_beta, Dz, u, var_v, nObs, nItems, nUsers, nFactors, obsIndex_item, oiStart_item, oiNum_item, debug);
+ //accept_fact_denom += *nItems; accept_fact_numeri += naccepts;
+ factor_condMeanVarSample_arscid(v, ars_XI_v, item, user, y, xb_plus_alpha_beta, Dz, u, var_v, nObs, nItems, nUsers, nFactors, obsIndex_item, oiStart_item, oiNum_item, ars_ninit, ars_qcent, ars_xl_v, ars_xu,ars_alpha,debug);
+ //if(*center==1){ center_array_2d(v, *nItems, *nFactors, 2); }
+ }
+ // Ignore the first several samples
+ if(s >= *nBurnin){
+ // update o
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ uv = 0;
+ for(f=1; f<=*nFactors; f++) uv += u[R_MAT(user_i,f,*nUsers)] * v[R_MAT(item_j,f,*nItems)];
+ double o = alpha[R_VEC(user_i)] + beta[R_VEC(item_j)] + uv;
+ o_mean[k] += o/(*nSamples);
+ }
+ // update alpha
+ for(k=0; k<*nUsers; k++){
+ alpha_sum[k] += alpha[k];
+ alpha_sos[k] += alpha[k]*alpha[k];
+ }
+ // update beta
+ for(k=0; k<*nItems; k++){
+ beta_sum[k] += beta[k];
+ beta_sos[k] += beta[k]*beta[k];
+ }
+ // update u
+ for(k=0; k<(*nUsers)*(*nFactors); k++){
+ u_sum[k] += u[k];
+ u_sos[k] += u[k]*u[k];
+ }
+ // update v
+ for(k=0; k<(*nItems)*(*nFactors); k++){
+ v_sum[k] += v[k];
+ v_sos[k] += v[k]*v[k];
+ }
+ }
+ }
+
+ computeMeanSumvar_arsc(alpha_mean, alpha_sumvar, alpha_sum, alpha_sos, *nUsers, *nSamples);
+ computeMeanSumvar_arsc(beta_mean, beta_sumvar, beta_sum, beta_sos, *nItems, *nSamples);
+ computeMeanSumvarFactor_arsc(u_mean, u_sumvar, u_sum, u_sos, *nUsers, *nFactors, *nSamples);
+ computeMeanSumvarFactor_arsc(v_mean, v_sumvar, v_sum, v_sos, *nItems, *nFactors, *nSamples);
+ //*acceptrate_maineff = accept_maineff_numeri/accept_maineff_denom;
+ //*acceptrate_fact = accept_fact_numeri/accept_fact_denom;
+
+ Free(alpha_sos);
+ Free(beta_sos);
+ Free(u_sos);
+ Free(v_sos);
+ Free(alpha_sum);
+ Free(beta_sum);
+ Free(u_sum);
+ Free(v_sum);
+
+ Free(obsIndex_user);
+ Free(oiStart_user);
+ Free(oiNum_user);
+ Free(obsIndex_item);
+ Free(oiStart_item);
+ Free(oiNum_item);
+ //Free(ars_XI_alpha);
+ //Free(ars_XI_beta);
+ //Free(ars_XI_u);
+ //Free(ars_XI_v);
+ Free(temp);
+}
+
View
399 src/RLFM-ars-logistic/C/MCEM_EStep_logistic.c~
@@ -0,0 +1,399 @@
+/*
+To Compile:
+R CMD SHLIB util.c MCEM_EStep.c -o MCEM_EStep.so
+
+*/
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include "util.h"
+//#include "logistic.h"
+#include "arsspline.h"
+
+///// centered logistic ars code starts from here
+// FUNCTION: mainEffect_condMeanVar
+// thisEffIndex: nObs x 1 otherEffIndex: nObs x 1 y_minus_xb_minus_uv: nObs x 1
+// fittedEff: nThisEff x 1 otherEff: nOtherEff x 1
+// For alpha, call
+// mainEffect_condMeanVar(1, user, item, y-xb-uv, g0w, beta, var_y, var_alpha, nObs, nUsers, nItems, sample, NULL, NULL, ...);
+// For beta, call
+// mainEffect_condMeanVar(1, item, user, y-xb-uv, d0z, alpha, var_y, var_beta, nObs, nItems, nUsers, sample, NULL, NULL, ...);
+// Observation index (consider, say, user i) - R indices (starting from 1, NOT 0)
+// obsIndex[ oiStart[i]+0 ], ..., obsIndex[ oiStart[i]+oiNum[i]-1 ] are the indices of user i's observations
+// in y, x, user, item
+void mainEffect_condMeanVarSample_arscid(
+ // OUTPUT
+ double* outSample, //int* naccepts,
+ double* ars_XI, /* ars_ninit X ncov X nThisEff*/
+ //INPUT
+ const int* thisEffIndex /*user or item*/, const int* otherEffIndex /*item or user*/,
+ const double* y, const double* offset, /*xb+uv*/
+ const double* fittedEff /*g0w or d0z*/,
+ const double* otherEff /*beta or alpha*/,
+ const double* var_eff /*var_alpha or var_beta*/,
+ const int* nObs, const int* nThisEff, const int* nOtherEff,
+ const int* obsIndex, const int* oiStart, const int* oiNum,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug
+
+ ){
+ int i,j,k,m, thisIndex, otherIndex, oIndex;
+ //*naccepts = 0;
+ GetRNGstate();
+ int neval = 0;
+ for(i=0; i<*nThisEff; i++){
+ thisIndex = i+1;
+ if(oiNum[i]>0){
+ double* y_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the observations for this user
+ double* offset_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the offsets for this user
+ double* X = (double*)calloc(oiNum[i], sizeof(double)); // the design matrix nobs*ncov, ncov=1 for main effects
+ for (j=0;j<oiNum[i];j++)
+ X[j] = 1;
+ for(j=0; j<oiNum[i]; j++){
+ oIndex = obsIndex[R_VEC(oiStart[i]+j)];
+ otherIndex = otherEffIndex[R_VEC(oIndex)];
+ if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
+ if(*debug > 0) CHK_R_INDEX(otherIndex, *nOtherEff);
+ if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
+ y_thisEff[j] = y[R_VEC(oIndex)];
+ offset_thisEff[j] = offset[R_VEC(oIndex)] + otherEff[R_VEC(otherIndex)];
+ }
+ //double new_thisEff = 0;
+ //int accept = 0;
+ int ncov = 1;
+ //MHlogistic(&outSample[i], X, y_thisEff, offset_thisEff, &fittedEff[i], var_eff, &oiNum[i], &ncov, &new_thisEff, &accept);
+ ARSLOGISTICSPLINE(offset_thisEff, &fittedEff[i], var_eff, X, y_thisEff, &ncov, &oiNum[i], &outSample[i], ars_qcent, ars_ninit, ars_ninit, ars_xl, ars_xu, &ars_XI[i*(*ars_ninit)], ars_alpha, &neval);
+ //printf("%f\n",outSample[i]);
+ //outSample[i] = new_thisEff;
+ //if (accept==1) *naccepts = *naccepts + 1;
+ free(y_thisEff);
+ free(offset_thisEff);
+ free(X);
+ }
+ }
+ PutRNGstate();
+}
+void factor_condMeanVarSample_arscid(
+ // OUTPUT
+ double* outSample, //int* naccepts,
+ double* ars_XI, /* ars_ninit X ncov X nThisEff*/
+ // INPUT
+ const int* thisEffIndex /*user or item*/, const int* otherEffIndex /*item or user*/,
+ const double* y, const double* offset, /*xb+alpha+beta*/
+ const double* fittedEff /*Gw or Dz*/,
+ const double* otherEff /*v or u*/,
+ const double* var_eff /*var_u or var_v*/,
+ const int* nObs, const int* nThisEff, const int* nOtherEff, const int* nFactors,
+ const int* obsIndex, const int* oiStart, const int* oiNum,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug
+ ){
+ int i,j,k,m, thisIndex, otherIndex, oIndex;
+ int* neval = (int*)calloc(*nFactors,sizeof(double));
+ //*naccepts = 0;
+ GetRNGstate();
+ for(i=0; i<*nThisEff; i++){
+ thisIndex = i+1;
+ if(oiNum[i]>0)
+ {
+ double* y_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the observations for this user
+ double* offset_thisEff = (double*)calloc(oiNum[i], sizeof(double)); // the offsets for this user
+ double* vj = (double*)calloc(oiNum[i]*(*nFactors), sizeof(double)); // the design matrix nobs*ncov, ncov=nfactors
+ double* thisEff_i = (double*)calloc((*nFactors),sizeof(double));
+ //double* new_thisEff_i = (double*)calloc((*nFactors),sizeof(double));
+ double* fittedEff_i = (double*)calloc((*nFactors),sizeof(double));
+ double* var_this_eff = (double*)calloc((*nFactors),sizeof(double));
+
+ for (j=0;j<*nFactors;j++)
+ {
+ thisEff_i[j] = outSample[C_MAT(i,j,*nThisEff)];
+ fittedEff_i[j] = fittedEff[C_MAT(i,j,*nThisEff)];
+ var_this_eff[j] = var_eff[j];
+ }
+ for(j=0; j<oiNum[i]; j++){
+ oIndex = obsIndex[R_VEC(oiStart[i]+j)];
+ otherIndex = otherEffIndex[R_VEC(oIndex)];
+ if(*debug > 0) CHK_R_INDEX(oIndex, *nObs);
+ if(*debug > 0) CHK_R_INDEX(otherIndex, *nOtherEff);
+ if(*debug > 1) if(thisEffIndex[R_VEC(oIndex)] != i+1) error("error in obsIndex, oiStart, oiNum\n");
+ y_thisEff[j] = y[R_VEC(oIndex)];
+ offset_thisEff[j] = offset[R_VEC(oIndex)];
+ for(k=1; k<=*nFactors; k++) {
+ vj[R_MAT((j+1),k,oiNum[i])] = otherEff[R_MAT(otherIndex,k,*nOtherEff)];
+ }
+ }
+ //double new_thisEff = 0;
+ // int accept = 0;
+ // MHlogistic(thisEff_i, vj, y_thisEff, offset_thisEff, fittedEff_i, var_eff, &oiNum[i], nFactors, new_thisEff_i, &accept);
+ //for (j=0;j<*nFactors;j++)
+ // printf("%f ",thisEff_i[j]);
+ //printf("\n");
+ ARSLOGISTICSPLINE(offset_thisEff, fittedEff_i, var_this_eff, vj, y_thisEff, nFactors, &oiNum[i], thisEff_i, ars_qcent, ars_ninit, ars_ninit, ars_xl, ars_xu, &ars_XI[i*(*ars_ninit)*(*nFactors)], ars_alpha, neval);
+
+ for (j=0;j<*nFactors;j++)
+ {
+ //printf("%f ",thisEff_i[j]);
+ outSample[C_MAT(i,j,*nThisEff)] = thisEff_i[j];
+ //if (accept==1) *naccepts = *naccepts + 1;
+ }
+ //printf("\n");
+ free(y_thisEff);
+ free(offset_thisEff);
+ free(vj);
+ free(thisEff_i);
+ //free(new_thisEff_i);
+ free(fittedEff_i);
+ free(var_this_eff);
+ }
+ }
+ PutRNGstate();
+ free(neval);
+}
+
+
+void computeMeanSumvar_arsc(
+ // OUTPUT
+ double *mean, double *sumvar,
+ // INPUT
+ const double *sum, const double *sos /* sum of squares */, const int length, const int nSamples
+ ){
+ int k;
+ double ratio = ((double)nSamples) / (nSamples - 1.0);
+ sumvar[0] = 0;
+ for(k=0; k<length; k++){
+ mean[k] = sum[k] / nSamples;
+ sumvar[0] += (sos[k] / (nSamples - 1.0)) - (ratio * mean[k] * mean[k]);
+ }
+}
+
+void computeMeanSumvarFactor_arsc(
+ // OUTPUT
+ double *mean, double *sumvar /* nFactors x 1*/,
+ // INPUT
+ const double *sum, const double *sos /* sum of squares */, const int length, const int nFactors, const int nSamples
+ ){
+ int k,l;
+ double ratio = ((double)nSamples) / (nSamples - 1.0);
+ for (k=0;k<nFactors;k++)
+ sumvar[k] = 0;
+ for (l=0; l<nFactors;l++)
+ for(k=0; k<length; k++){
+ mean[C_MAT(k,l,length)] = sum[C_MAT(k,l,length)] / nSamples;
+ sumvar[l] += (sos[C_MAT(k,l,length)] / (nSamples - 1.0)) - (ratio * mean[C_MAT(k,l,length)] * mean[C_MAT(k,l,length)]);
+ }
+}
+
+void InitializeXI_arsc(double* xi, const int ninit, const int ncov, const int nthisEff, const double xu, const double xl)
+{
+ int i,j,k;
+ for (i=0;i<nthisEff;i++)
+ for (j=0;j<ncov;j++)
+ for (k=0;k<ninit;k++)
+ {
+ double tmp = xl + (k + 2.0)*(xu - xl)/(ninit + 1.0);
+ if (tmp>=xu) tmp = tmp - .1;
+ if (tmp<=xl) tmp = tmp + .1;
+ xi[i*(ncov*ninit)+j*ninit+k] = tmp;
+ }
+}
+void MCEM_EStep_arscid(
+ // OUTPUT
+ double* o_mean /*nObs x 1*/,
+ double* alpha_mean/*nUsers x 1*/, double* alpha_sumvar/*1x1*/,
+ double* beta_mean/*nItems x 1*/, double* beta_sumvar/*1x1*/,
+ double* u_mean/*nUsers x nFactors*/, double* u_sumvar/*nFactors x 1*/,
+ double* v_mean/*nItems x nFactors*/, double* v_sumvar/*nFactors x 1*/,
+ double* ars_XI_alpha, double* ars_XI_beta,
+ double* ars_XI_u, double* ars_XI_v,
+ //double* acceptrate_maineff, // the overall acceptance rate for main effects
+ //double* acceptrate_fact, // the overall acceptance rate for factors
+ // INPUT
+ const int* nSamples, const int* nBurnin,
+ const int* user/*nObs x 1*/, const int* item/*nObs x 1*/,
+ const double* y/*nObs x 1*/, const double* xb/*nObs x 1*/,
+ const double* g0w/*nUsers x 1*/, const double* d0z/*nItems x 1*/,
+ const double* Gw/*nUsers x nFactors*/, const double* Dz/*nItems x nFactors*/,
+ const double* alpha_in/*nUsers x 1*/, const double* beta_in/*nItems x 1*/,
+ const double* u_in/*nUsers x nFactors*/, const double* v_in/*nItems x nFactors*/,
+ const double* var_alpha, const double* var_beta,
+ const double* var_u, const double* var_v, /*nFactors x 1*/
+ const int* nObs, const int* nUsers, const int* nItems, const int* nFactors,
+ const int* ars_ninit, const double* ars_qcent, const double* ars_xl, const double* ars_xu, const double* ars_alpha,
+ // OTHER
+ const int* debug, const int* main_effects, const int* beta_int, const int* center
+
+ ){
+ int *obsIndex_user, *oiStart_user, *oiNum_user,
+ *obsIndex_item, *oiStart_item, *oiNum_item,
+ s, k, f, user_i, item_j, naccepts;
+ double *alpha, *beta, *u, *v,
+ *alpha_sum, *alpha_sos, *beta_sum, *beta_sos, *u_sum, *u_sos, *v_sum, *v_sos,
+ *temp, *xb_plus_uv, *xb_plus_alpha_beta, uv;
+ double *ars_xl_v = (double*)calloc(1, sizeof(double));
+ ars_xl_v[0] = 0;
+
+ // double *ars_XI_alpha, *ars_XI_beta, *ars_XI_u, *ars_XI_v;
+ //double accept_maineff_denom = 0, accept_maineff_numeri = 0,
+ // accept_fact_denom = 0, accept_fact_numeri = 0;
+
+
+ // Allocate space for sum and sum-of-squares
+ alpha_sum = (double*)calloc(*nUsers, sizeof(double));
+ beta_sum = (double*)calloc(*nItems, sizeof(double));
+ u_sum = (double*)calloc((*nUsers)*(*nFactors), sizeof(double));
+ v_sum = (double*)calloc((*nItems)*(*nFactors), sizeof(double));
+ alpha_sos = (double*)calloc(*nUsers, sizeof(double));
+ beta_sos = (double*)calloc(*nItems, sizeof(double));
+ u_sos = (double*)calloc((*nUsers)*(*nFactors), sizeof(double));
+ v_sos = (double*)calloc((*nItems)*(*nFactors), sizeof(double));
+ // Allocate space for the observation indices
+ obsIndex_user = (int*)calloc(*nObs, sizeof(int));
+ oiStart_user = (int*)calloc(*nUsers, sizeof(int));
+ oiNum_user = (int*)calloc(*nUsers, sizeof(int));
+ obsIndex_item = (int*)calloc(*nObs, sizeof(int));
+ oiStart_item = (int*)calloc(*nItems, sizeof(int));
+ oiNum_item = (int*)calloc(*nItems, sizeof(int));
+ // Allocate space for XI ars
+ //ars_XI_alpha = (double*)calloc((*ars_ninit)*(*nUsers), sizeof(double));
+ //ars_XI_beta = (double*)calloc((*ars_ninit)*(*nItems), sizeof(double));
+ //ars_XI_u = (double*)calloc((*ars_ninit)*(*nFactors)*(*nUsers), sizeof(double));
+ //ars_XI_v = (double*)calloc((*ars_ninit)*(*nFactors)*(*nItems), sizeof(double));
+
+ // Allocate temp space
+ temp = (double*)calloc(*nObs, sizeof(double));
+
+ memset(o_mean,0,(*nObs)*sizeof(double));
+ // Use the memory space of the output to store the current alpha, beta, u and v
+ alpha = alpha_mean; beta = beta_mean; u = u_mean; v = v_mean;
+ // Use the temp space for both xb_plus_uv and xb_plus_alpha_beta
+ xb_plus_uv = temp;
+ xb_plus_alpha_beta = temp;
+
+
+ // Create Observation indices for users and items
+ generateObsIndex(obsIndex_user, oiStart_user, oiNum_user, user, nObs, nUsers, debug);
+ generateObsIndex(obsIndex_item, oiStart_item, oiNum_item, item, nObs, nItems, debug);
+
+ // Initialize alpha, beta, u, v
+ for(k=0; k<*nUsers; k++) alpha[k] = alpha_in[k];
+ for(k=0; k<*nItems; k++) beta[k] = beta_in[k];
+ for(k=0; k<(*nUsers)*(*nFactors); k++) u[k] = u_in[k];
+ for(k=0; k<(*nItems)*(*nFactors); k++) v[k] = v_in[k];
+
+ // Initialize ars_XI_alpha, etc.
+ InitializeXI_arsc(ars_XI_alpha, *ars_ninit, 1, *nUsers, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_beta, *ars_ninit, 1, *nItems, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_u, *ars_ninit, *nFactors, *nUsers, *ars_xu, *ars_xl);
+ InitializeXI_arsc(ars_XI_v, *ars_ninit, *nFactors, *nItems, *ars_xu, *ars_xl_v);
+
+ for(s=0; s<(*nSamples+*nBurnin); s++){
+
+ // Compute xb+uv
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ uv = 0;
+ for(f=1; f<=*nFactors; f++) uv += u[R_MAT(user_i,f,*nUsers)] * v[R_MAT(item_j,f,*nItems)];
+ xb_plus_uv[k] = xb[k] + uv;
+ }
+ // Sample alpha
+ //mainEffect_condMeanVarSample_arsc(alpha, &naccepts, user, item, y, xb_plus_uv, g0w, beta, var_alpha, nObs, nUsers, nItems, obsIndex_user, oiStart_user, oiNum_user, debug);
+ //accept_maineff_denom += *nUsers; accept_maineff_numeri += naccepts;
+ mainEffect_condMeanVarSample_arscid(alpha, ars_XI_alpha, user, item, y, xb_plus_uv, g0w, beta, var_alpha, nObs, nUsers, nItems, obsIndex_user, oiStart_user, oiNum_user, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ //center
+ if(*center==1){ center_array(alpha, *nUsers); }
+ // Sample beta
+ //mainEffect_condMeanVarSample_arsc(beta, &naccepts, item, user, y, xb_plus_uv, d0z, alpha, var_beta, nObs, nItems, nUsers, obsIndex_item, oiStart_item, oiNum_item, debug);
+ //accept_maineff_denom += *nItems; accept_maineff_numeri += naccepts;
+ mainEffect_condMeanVarSample_arscid(beta, ars_XI_beta, item, user, y, xb_plus_uv, d0z, alpha, var_beta, nObs, nItems, nUsers, obsIndex_item, oiStart_item, oiNum_item, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ //subtract mean from the betas ...
+ if(*center==1 && *beta_int==0){ center_array(beta, *nItems); }
+
+ // Compute y - (xb + alpha + beta)
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ xb_plus_alpha_beta[k] = xb[k] + alpha[R_VEC(user_i)] + beta[R_VEC(item_j)];
+ }
+
+ if(*main_effects==0){
+ // Sample u
+ //factor_condMeanVarSample_arsc(u, &naccepts, user, item, y, xb_plus_alpha_beta, Gw, v, var_u, nObs, nUsers, nItems, nFactors, obsIndex_user, oiStart_user, oiNum_user, debug);
+ //accept_fact_denom += *nUsers; accept_fact_numeri += naccepts;
+ factor_condMeanVarSample_arscid(u, ars_XI_u, user, item, y, xb_plus_alpha_beta, Gw, v, var_u, nObs, nUsers, nItems, nFactors, obsIndex_user, oiStart_user, oiNum_user, ars_ninit, ars_qcent, ars_xl, ars_xu,ars_alpha,debug);
+ if(*center==1){ center_array_2d(u, *nUsers, *nFactors, 2); }
+ // Sample v
+ //factor_condMeanVarSample_arsc(v, &naccepts, item, user, y, xb_plus_alpha_beta, Dz, u, var_v, nObs, nItems, nUsers, nFactors, obsIndex_item, oiStart_item, oiNum_item, debug);
+ //accept_fact_denom += *nItems; accept_fact_numeri += naccepts;
+ factor_condMeanVarSample_arscid(v, ars_XI_v, item, user, y, xb_plus_alpha_beta, Dz, u, var_v, nObs, nItems, nUsers, nFactors, obsIndex_item, oiStart_item, oiNum_item, ars_ninit, ars_qcent, ars_xl_v, ars_xu,ars_alpha,debug);
+ //if(*center==1){ center_array_2d(v, *nItems, *nFactors, 2); }
+ }
+ // Ignore the first several samples
+ if(s >= *nBurnin){
+ // update o
+ for(k=0; k<*nObs; k++){
+ user_i = user[k]; item_j = item[k];
+ if(*debug > 0) CHK_R_INDEX(user_i, *nUsers);
+ if(*debug > 0) CHK_R_INDEX(item_j, *nItems);
+ uv = 0;
+ for(f=1; f<=*nFactors; f++) uv += u[R_MAT(user_i,f,*nUsers)] * v[R_MAT(item_j,f,*nItems)];
+ double o = alpha[R_VEC(user_i)] + beta[R_VEC(item_j)] + uv;
+ o_mean[k] += o/(*nSamples);
+ }
+ // update alpha
+ for(k=0; k<*nUsers; k++){
+ alpha_sum[k] += alpha[k];
+ alpha_sos[k] += alpha[k]*alpha[k];
+ }
+ // update beta
+ for(k=0; k<*nItems; k++){
+ beta_sum[k] += beta[k];
+ beta_sos[k] += beta[k]*beta[k];
+ }
+ // update u
+ for(k=0; k<(*nUsers)*(*nFactors); k++){
+ u_sum[k] += u[k];
+ u_sos[k] += u[k]*u[k];
+ }
+ // update v
+ for(k=0; k<(*nItems)*(*nFactors); k++){
+ v_sum[k] += v[k];
+ v_sos[k] += v[k]*v[k];
+ }
+ }
+ }
+
+ computeMeanSumvar_arsc(alpha_mean, alpha_sumvar, alpha_sum, alpha_sos, *nUsers, *nSamples);
+ computeMeanSumvar_arsc(beta_mean, beta_sumvar, beta_sum, beta_sos, *nItems, *nSamples);
+ computeMeanSumvarFactor_arsc(u_mean, u_sumvar, u_sum, u_sos, *nUsers, *nFactors, *nSamples);
+ computeMeanSumvarFactor_arsc(v_mean, v_sumvar, v_sum, v_sos, *nItems, *nFactors, *nSamples);
+ //*acceptrate_maineff = accept_maineff_numeri/accept_maineff_denom;
+ //*acceptrate_fact = accept_fact_numeri/accept_fact_denom;
+
+ Free(alpha_sos);
+ Free(beta_sos);
+ Free(u_sos);
+ Free(v_sos);
+ Free(alpha_sum);
+ Free(beta_sum);
+ Free(u_sum);
+ Free(v_sum);
+
+ Free(obsIndex_user);
+ Free(oiStart_user);
+ Free(oiNum_user);
+ Free(obsIndex_item);
+ Free(oiStart_item);
+ Free(oiNum_item);
+ //Free(ars_XI_alpha);
+ //Free(ars_XI_beta);
+ //Free(ars_XI_u);
+ //Free(ars_XI_v);
+ Free(temp);
+}
+
View
506 src/RLFM-ars-logistic/C/util.c
@@ -0,0 +1,506 @@
+/*
+ Copyright (c) 2011, Yahoo! Inc. All rights reserved.
+ Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
+
+ Author: Liang Zhang
+*/
+
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include <R_ext/Applic.h>
+#include "util.h"
+
+void sym_eigen(const double* x, const int *nrow, double *eigen_val, double *eigen_vec){
+ char jobz = 'V', uplo = 'L';
+ double work_size, *work;
+ int i,j, info, lwork=-1;
+
+ for(i=0; i<(*nrow)*(*nrow); i++) eigen_vec[i] = x[i];
+
+ F77_NAME(dsyev)(&jobz, &uplo, nrow, eigen_vec, nrow, eigen_val, &work_size, &lwork, &info);
+ if(info != 0) error("error in dsyev(...)");
+
+ lwork = work_size;
+ work = (double*)Calloc(lwork, double);
+ F77_NAME(dsyev)(&jobz, &uplo, nrow, eigen_vec, nrow, eigen_val, work, &lwork, &info);
+ if(info != 0) error("error in dsyev(...)");
+
+ Free(work);
+}
+
+void sym_inv_byCholesky(
+ double *A /* n x n matrix */, const int *n, const int *check_sym
+){
+ char uplo = 'L';
+ int info, i, j;
+
+ if(*check_sym > 0) CHK_SYMMETRIC(A, *n, i, j);
+
+ F77_NAME(dpotrf)(&uplo, n, A, n, &info);
+ if(info != 0) error("error in dpotrf(...): info=%d", info);
+ F77_NAME(dpotri)(&uplo, n, A, n, &info);
+ if(info != 0) error("error in dpotri(...): info=%d", info);
+ for(i=1; i<(*n); i++){
+ for(j=0; j<i; j++){
+ A[C_MAT(j,i,*n)] = A[C_MAT(i,j,*n)];
+ }
+ }
+}
+
+void sym_inv3DA_byCholesky(
+ double *invA /* k x n x n matrix */, const double *A /* k x n x n matrix */,
+ const int *k, const int *n, double *temp /* n x n */, const int *check_sym
+){
+ for(int i=0; i<*k; i++){
+ for(int j=0; j<*n; j++)
+ for(int m=0; m<*n; m++) temp[C_MAT(j,m,*n)] = A[C_3DA(i,j,m,*k,*n)];
+ sym_inv_byCholesky(temp, n, check_sym);
+
+ for(int j=0; j<*n; j++)
+ for(int m=0; m<*n; m++) invA[C_3DA(i,j,m,*k,*n)] = temp[C_MAT(j,m,*n)];
+ }
+}
+
+void sum_margin(
+ // OUTPUT
+ double *ans,
+ // INPUT
+ const double *A, const int *nrow, const int *ncol,
+ const int *side // side=1: Sum up each row and return a vector with length nrow
+ // side=2: Sum up each column and return a vector with length ncol
+){
+ int i, j, end;
+ if((*side) == 1){
+ end=(*nrow)*(*ncol);
+ for(i=0; i<*nrow; i++) ans[i] = 0;
+ for(j=0; j<end; j+=(*nrow))
+ for(i=0; i<*nrow; i++) ans[i] += A[i+j];
+ }else if((*side) == 2){
+ for(j=0; j<*ncol; j++){
+ end = (j+1)*(*nrow);
+ ans[j] = 0;
+ for(i=j*(*nrow); i<end; i++) ans[j] += A[i];
+ }
+ }else{
+ error("Unknown side=%d (please specify side=1 (for rows) or side=2 (for columns)");
+ }
+}
+
+void center_array(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n // array length
+ ){
+ // double sum = 0.0;
+ double mean = x[0];
+ // get sum
+ //for(int i=0; i<n; i++){sum += x[i];}
+ //mean = sum / (double) n;
+ //get mean ... more stable ...
+ for(int i=1; i<n; i++){
+ mean = ( (double) i / (double) (i+1) ) * mean + ( 1.0 / (double) (i+1) ) * x[i];
+ }
+
+ // Rprintf("mean = %f\n", mean);
+ // center
+ for(int i=0; i<n; i++){x[i] -= mean;}
+}
+
+void center_array_2d(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim // 1 for subtracting row mean, 2 for column
+){
+
+ double sum, mean;
+ int i,j;
+
+ if(dim == 1){
+ for(i = 0; i<n; i++){
+ sum = 0.0;
+ for(j = 0; j<m; j++){sum += x[C_MAT(i, j, n)];}
+ mean = sum / (double) m;
+ for(j = 0; j<m; j++){x[C_MAT(i, j, n)] -= mean;}
+ }
+ } else if(dim == 2){
+ for(j = 0; j<m; j++){
+ //sum = 0.0;
+ //for(i = 0; i<n; i++){sum += x[C_MAT(i, j, n)];}
+ //mean = sum / (double) n;
+
+ mean = x[C_MAT(0, j, n)];
+ for(i=1; i<n; i++){
+ mean = ((double)i/(double)(i+1)) * mean + (1.0/(double) (i+1))*x[C_MAT(i,j,n)];
+ }
+ // Rprintf("fac mean col %d = %f\n", j, mean);
+ for(i = 0; i<n; i++){x[C_MAT(i, j, n)] -= mean;}
+ }
+
+ }
+}
+
+void center_array_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // array length
+ int * oiNum
+ ){
+ // double sum = 0.0;
+ double mean = 0.0;
+ int nav = 0;
+ // get sum
+ //for(int i=0; i<n; i++){sum += x[i];}
+ //mean = sum / (double) n;
+ //get mean ... more stable ...
+ for(int i=0; i<n; i++){
+ if(oiNum[i]>0){
+ mean = ( (double) nav / (double) (nav+1) ) * mean + \
+ ( 1.0 / (double) (nav+1) ) * x[i];
+ nav += 1;
+ }
+ }
+
+ // Rprintf("mean = %f\n", mean);
+ // center
+ for(int i=0; i<n; i++){
+ if(oiNum[i]>0){
+ x[i] -= mean;
+ }
+ }
+}
+
+void center_array_2d_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim, // 1 for subtracting row mean, 2 for column
+ int * oiNum
+){
+
+ double sum, mean;
+ int i,j;
+ int nav;
+
+ // dim == 1 is broken!!!
+ if(dim == 1){
+ for(i = 0; i<n; i++){
+ if(oiNum[i]>0){
+ sum = 0.0;
+ for(j = 0; j<m; j++){sum += x[C_MAT(i, j, n)];}
+ mean = sum / (double) m;
+ for(j = 0; j<m; j++){x[C_MAT(i, j, n)] -= mean;}
+ }
+ }
+ } else if(dim == 2){
+ for(j = 0; j<m; j++){
+ //sum = 0.0;
+ //for(i = 0; i<n; i++){sum += x[C_MAT(i, j, n)];}
+ //mean = sum / (double) n;
+ nav = 0;
+ mean = 0.0;
+ for(i=0; i<n; i++){
+ if(oiNum[i]>0){
+ mean = ((double)nav/(double)(nav+1)) * mean + \
+ (1.0/(double) (nav+1))*x[C_MAT(i,j,n)];
+ nav += 1;
+ }
+ }
+ // Rprintf("fac mean col %d = %f\n", j, mean);
+ for(i = 0; i<n; i++){
+ if(oiNum[i]>0){
+ x[C_MAT(i, j, n)] -= mean;
+ }
+ }
+ }
+
+ }
+}
+
+
+void print_vector(char* prefix, double* vector, int length){
+ int i;
+ if(prefix != NULL) Rprintf("%s", prefix);
+ for(i=0; i<length; i++)
+ Rprintf("%f ", vector[i]);
+ Rprintf("\n");
+}
+
+void print_matrix(const char* prefix, const double* matrix, const int nrow, const int ncol){
+ int i,j;
+ for(i=0; i<nrow; i++){
+ if(prefix != NULL) Rprintf("%s", prefix);
+ for(j=0; j<ncol; j++)
+ Rprintf("%f\t", matrix[C_MAT(i,j,nrow)]);
+ Rprintf("\n");
+ }
+}
+
+// The intput/output are all R indices (start from 1, NOT 0)
+void generateObsIndex(
+ //OUTPUT
+ int* obsIndex, // E.g., consider user i
+ int* start, // y[ obsIndex[ start[i]+(0:(num[i]-1)) ] ]
+ int* num, // are the observations of user i
+ //INPUT
+ const int* effIndex /* user or item */,
+ const int* nObs, const int* nEff,
+ //OTHER
+ const int* debug
+){
+ int *ind, i, j, cIndex, eIndex, oIndex;
+ ind = (int*)Calloc(*nEff, int);
+
+ for(i=0; i<*nObs; i++){
+ eIndex = effIndex[i];
+ if(*debug > 0){
+ CHK_R_INDEX(eIndex, *nEff);
+ }
+ num[R_VEC(eIndex)]++;
+ }
+
+ start[0] = 1; ind[0] = 1;
+ for(i=1; i<*nEff; i++){
+ start[i] = start[i-1]+num[i-1];
+ ind[i] = start[i];
+ }
+
+ for(i=0; i<*nObs; i++){
+ cIndex = R_VEC(effIndex[i]);
+ obsIndex[R_VEC(ind[cIndex])] = i+1;
+ ind[cIndex]++;
+ }
+
+ if(*debug > 0){
+ for(i=0; i<*nEff; i++){
+ if(ind[i] != start[i]+num[i]) error("logical error (level 1)");
+ if(*debug > 1){
+ for(j=0; j<num[i]; j++){
+ oIndex = obsIndex[R_VEC(start[i]+j)];
+ if(effIndex[R_VEC(oIndex)] != i+1)
+ error("logical error (level 2)");
+ }
+ }
+ }
+ }
+
+ Free(ind);
+}
+
+#define stepredn 0.2
+#define acctol 0.0001
+#define reltest 10.0
+#define _(A) A
+
+static double * vect(int n)
+{
+ return (double *)R_alloc(n, sizeof(double));
+}
+
+int nParam_debug = 0;
+double* gr_debug = NULL;
+
+void my_cgmin(int n, double *Bvec, double *X, double *Fmin,
+ optimfn fminfn, optimgr fmingr, int *fail,
+ double abstol, double intol, void *ex, int type, int trace,
+ int *fncount, int *grcount, int maxit)
+{
+ Rboolean accpoint;
+ double *c, *g, *t;
+ int count, cycle, cyclimit;
+ double f;
+ double G1, G2, G3, gradproj;
+ int funcount=0, gradcount=0, i;
+ double newstep, oldstep, setstep, steplength=1.0;
+ double tol;
+
+ if (maxit <= 0) {
+ *Fmin = fminfn(n, Bvec, ex);
+ *fncount = *grcount = 0;
+ *fail = FALSE;
+ return;
+ }
+ if (trace) {
+ Rprintf(" Conjugate gradients function minimizer\n");
+ switch (type) {
+ case 1: Rprintf("Method: Fletcher Reeves\n"); break;
+ case 2: Rprintf("Method: Polak Ribiere\n"); break;
+ case 3: Rprintf("Method: Beale Sorenson\n"); break;
+ default:
+ error(_("unknown 'type' in CG method of optim"));
+ }
+ }
+ c = vect(n); g = vect(n); t = vect(n);
+
+ // DEBUG
+ nParam_debug = n;
+ gr_debug = g;
+ // END
+
+ print_vector("INITIAL: ", gr_debug, nParam_debug);
+
+ setstep = 1.7;
+ *fail = 0;
+ cyclimit = n;
+ tol = intol * n * sqrt(intol);
+
+ if (trace) Rprintf("tolerance used in gradient test=%g\n", tol);
+ f = fminfn(n, Bvec, ex);
+ print_vector("AFTER fminfn IN cgmin: ", gr_debug, nParam_debug);
+ if (!R_FINITE(f)) {
+ error(_("Function cannot be evaluated at initial parameters"));
+ } else {
+ *Fmin = f;
+ funcount = 1;
+ gradcount = 0;
+ do {
+ for (i = 0; i < n; i++) {
+ t[i] = 0.0;
+ c[i] = 0.0;
+ }
+ cycle = 0;
+ oldstep = 1.0;
+ count = 0;
+ do {
+ cycle++;
+ if (trace) {
+ Rprintf("gradcount=%d, funcount=%d, Fmin=%f\n", gradcount, funcount, *Fmin);
+ Rprintf("parameters:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", Bvec[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ gradcount++;
+ if (gradcount > maxit) {
+ *fncount = funcount;
+ *grcount = gradcount;
+ *fail = 1;
+ return;
+ }
+ if (trace) {
+ Rprintf("BEFORE: gradient:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", g[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ fmingr(n, Bvec, g, ex);
+ if (trace) {
+ Rprintf("AFTER: gradient:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", g[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ G1 = 0.0;
+ G2 = 0.0;
+ for (i = 0; i < n; i++) {
+ X[i] = Bvec[i];
+ switch (type) {
+
+ case 1: /* Fletcher-Reeves */
+ G1 += g[i] * g[i];
+ G2 += c[i] * c[i];
+ break;
+
+ case 2: /* Polak-Ribiere */
+ G1 += g[i] * (g[i] - c[i]);
+ G2 += c[i] * c[i];
+ break;
+
+ case 3: /* Beale-Sorenson */
+ G1 += g[i] * (g[i] - c[i]);
+ G2 += t[i] * (g[i] - c[i]);
+ break;
+
+ default:
+ error(_("unknown type in CG method of optim"));
+ }
+ c[i] = g[i];
+ }
+ if (G1 > tol) {
+ if (G2 > 0.0)
+ G3 = G1 / G2;
+ else
+ G3 = 1.0;
+ gradproj = 0.0;
+ for (i = 0; i < n; i++) {
+ t[i] = t[i] * G3 - g[i];
+ gradproj += t[i] * g[i];
+ }
+ steplength = oldstep;
+
+ accpoint = FALSE;
+ do {
+ count = 0;
+ for (i = 0; i < n; i++) {
+ Bvec[i] = X[i] + steplength * t[i];
+ if (reltest + X[i] == reltest + Bvec[i])
+ count++;
+ }
+ if (count < n) { /* point is different */
+ f = fminfn(n, Bvec, ex);
+ funcount++;
+ accpoint = (R_FINITE(f) &&
+ f <= *Fmin + gradproj * steplength * acctol);
+
+ if (!accpoint) {
+ steplength *= stepredn;
+ if (trace){
+ Rprintf("* %f\n", f);
+ }
+ } else *Fmin = f; /* we improved, so update value */
+ }
+ } while (!(count == n || accpoint));
+ if (count < n) {
+ newstep = 2 * (f - *Fmin - gradproj * steplength);
+ if (newstep > 0) {
+ newstep = -(gradproj * steplength * steplength / newstep);
+ for (i = 0; i < n; i++)
+ Bvec[i] = X[i] + newstep * t[i];
+ *Fmin = f;
+ f = fminfn(n, Bvec, ex);
+ funcount++;
+ if (f < *Fmin) {
+ *Fmin = f;
+ if (trace) Rprintf("i< ");
+ } else { /* reset Bvec to match lowest point */
+ if (trace) Rprintf("i> ");
+ for (i = 0; i < n; i++)
+ Bvec[i] = X[i] + steplength * t[i];
+ }
+ }
+ }
+ }
+ oldstep = setstep * steplength;
+ if (oldstep > 1.0)
+ oldstep = 1.0;
+ } while ((count != n) && (G1 > tol) && (cycle != cyclimit));
+
+ } while ((cycle != 1) ||
+ ((count != n) && (G1 > tol) && *Fmin > abstol));
+
+ }
+ if (trace) {
+ Rprintf("Exiting from conjugate gradients minimizer\n");
+ Rprintf(" %d function evaluations used\n", funcount);
+ Rprintf(" %d gradient evaluations used\n", gradcount);
+ }
+ *fncount = funcount;
+ *grcount = gradcount;
+}
View
500 src/RLFM-ars-logistic/C/util.c~
@@ -0,0 +1,500 @@
+
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include <R_ext/Applic.h>
+#include "util.h"
+
+void sym_eigen(const double* x, const int *nrow, double *eigen_val, double *eigen_vec){
+ char jobz = 'V', uplo = 'L';
+ double work_size, *work;
+ int i,j, info, lwork=-1;
+
+ for(i=0; i<(*nrow)*(*nrow); i++) eigen_vec[i] = x[i];
+
+ F77_NAME(dsyev)(&jobz, &uplo, nrow, eigen_vec, nrow, eigen_val, &work_size, &lwork, &info);
+ if(info != 0) error("error in dsyev(...)");
+
+ lwork = work_size;
+ work = (double*)Calloc(lwork, double);
+ F77_NAME(dsyev)(&jobz, &uplo, nrow, eigen_vec, nrow, eigen_val, work, &lwork, &info);
+ if(info != 0) error("error in dsyev(...)");
+
+ Free(work);
+}
+
+void sym_inv_byCholesky(
+ double *A /* n x n matrix */, const int *n, const int *check_sym
+){
+ char uplo = 'L';
+ int info, i, j;
+
+ if(*check_sym > 0) CHK_SYMMETRIC(A, *n, i, j);
+
+ F77_NAME(dpotrf)(&uplo, n, A, n, &info);
+ if(info != 0) error("error in dpotrf(...): info=%d", info);
+ F77_NAME(dpotri)(&uplo, n, A, n, &info);
+ if(info != 0) error("error in dpotri(...): info=%d", info);
+ for(i=1; i<(*n); i++){
+ for(j=0; j<i; j++){
+ A[C_MAT(j,i,*n)] = A[C_MAT(i,j,*n)];
+ }
+ }
+}
+
+void sym_inv3DA_byCholesky(
+ double *invA /* k x n x n matrix */, const double *A /* k x n x n matrix */,
+ const int *k, const int *n, double *temp /* n x n */, const int *check_sym
+){
+ for(int i=0; i<*k; i++){
+ for(int j=0; j<*n; j++)
+ for(int m=0; m<*n; m++) temp[C_MAT(j,m,*n)] = A[C_3DA(i,j,m,*k,*n)];
+ sym_inv_byCholesky(temp, n, check_sym);
+
+ for(int j=0; j<*n; j++)
+ for(int m=0; m<*n; m++) invA[C_3DA(i,j,m,*k,*n)] = temp[C_MAT(j,m,*n)];
+ }
+}
+
+void sum_margin(
+ // OUTPUT
+ double *ans,
+ // INPUT
+ const double *A, const int *nrow, const int *ncol,
+ const int *side // side=1: Sum up each row and return a vector with length nrow
+ // side=2: Sum up each column and return a vector with length ncol
+){
+ int i, j, end;
+ if((*side) == 1){
+ end=(*nrow)*(*ncol);
+ for(i=0; i<*nrow; i++) ans[i] = 0;
+ for(j=0; j<end; j+=(*nrow))
+ for(i=0; i<*nrow; i++) ans[i] += A[i+j];
+ }else if((*side) == 2){
+ for(j=0; j<*ncol; j++){
+ end = (j+1)*(*nrow);
+ ans[j] = 0;
+ for(i=j*(*nrow); i<end; i++) ans[j] += A[i];
+ }
+ }else{
+ error("Unknown side=%d (please specify side=1 (for rows) or side=2 (for columns)");
+ }
+}
+
+void center_array(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n // array length
+ ){
+ // double sum = 0.0;
+ double mean = x[0];
+ // get sum
+ //for(int i=0; i<n; i++){sum += x[i];}
+ //mean = sum / (double) n;
+ //get mean ... more stable ...
+ for(int i=1; i<n; i++){
+ mean = ( (double) i / (double) (i+1) ) * mean + ( 1.0 / (double) (i+1) ) * x[i];
+ }
+
+ // Rprintf("mean = %f\n", mean);
+ // center
+ for(int i=0; i<n; i++){x[i] -= mean;}
+}
+
+void center_array_2d(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim // 1 for subtracting row mean, 2 for column
+){
+
+ double sum, mean;
+ int i,j;
+
+ if(dim == 1){
+ for(i = 0; i<n; i++){
+ sum = 0.0;
+ for(j = 0; j<m; j++){sum += x[C_MAT(i, j, n)];}
+ mean = sum / (double) m;
+ for(j = 0; j<m; j++){x[C_MAT(i, j, n)] -= mean;}
+ }
+ } else if(dim == 2){
+ for(j = 0; j<m; j++){
+ //sum = 0.0;
+ //for(i = 0; i<n; i++){sum += x[C_MAT(i, j, n)];}
+ //mean = sum / (double) n;
+
+ mean = x[C_MAT(0, j, n)];
+ for(i=1; i<n; i++){
+ mean = ((double)i/(double)(i+1)) * mean + (1.0/(double) (i+1))*x[C_MAT(i,j,n)];
+ }
+ // Rprintf("fac mean col %d = %f\n", j, mean);
+ for(i = 0; i<n; i++){x[C_MAT(i, j, n)] -= mean;}
+ }
+
+ }
+}
+
+void center_array_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // array length
+ int * oiNum
+ ){
+ // double sum = 0.0;
+ double mean = 0.0;
+ int nav = 0;
+ // get sum
+ //for(int i=0; i<n; i++){sum += x[i];}
+ //mean = sum / (double) n;
+ //get mean ... more stable ...
+ for(int i=0; i<n; i++){
+ if(oiNum[i]>0){
+ mean = ( (double) nav / (double) (nav+1) ) * mean + \
+ ( 1.0 / (double) (nav+1) ) * x[i];
+ nav += 1;
+ }
+ }
+
+ // Rprintf("mean = %f\n", mean);
+ // center
+ for(int i=0; i<n; i++){
+ if(oiNum[i]>0){
+ x[i] -= mean;
+ }
+ }
+}
+
+void center_array_2d_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim, // 1 for subtracting row mean, 2 for column
+ int * oiNum
+){
+
+ double sum, mean;
+ int i,j;
+ int nav;
+
+ // dim == 1 is broken!!!
+ if(dim == 1){
+ for(i = 0; i<n; i++){
+ if(oiNum[i]>0){
+ sum = 0.0;
+ for(j = 0; j<m; j++){sum += x[C_MAT(i, j, n)];}
+ mean = sum / (double) m;
+ for(j = 0; j<m; j++){x[C_MAT(i, j, n)] -= mean;}
+ }
+ }
+ } else if(dim == 2){
+ for(j = 0; j<m; j++){
+ //sum = 0.0;
+ //for(i = 0; i<n; i++){sum += x[C_MAT(i, j, n)];}
+ //mean = sum / (double) n;
+ nav = 0;
+ mean = 0.0;
+ for(i=0; i<n; i++){
+ if(oiNum[i]>0){
+ mean = ((double)nav/(double)(nav+1)) * mean + \
+ (1.0/(double) (nav+1))*x[C_MAT(i,j,n)];
+ nav += 1;
+ }
+ }
+ // Rprintf("fac mean col %d = %f\n", j, mean);
+ for(i = 0; i<n; i++){
+ if(oiNum[i]>0){
+ x[C_MAT(i, j, n)] -= mean;
+ }
+ }
+ }
+
+ }
+}
+
+
+void print_vector(char* prefix, double* vector, int length){
+ int i;
+ if(prefix != NULL) Rprintf("%s", prefix);
+ for(i=0; i<length; i++)
+ Rprintf("%f ", vector[i]);
+ Rprintf("\n");
+}
+
+void print_matrix(const char* prefix, const double* matrix, const int nrow, const int ncol){
+ int i,j;
+ for(i=0; i<nrow; i++){
+ if(prefix != NULL) Rprintf("%s", prefix);
+ for(j=0; j<ncol; j++)
+ Rprintf("%f\t", matrix[C_MAT(i,j,nrow)]);
+ Rprintf("\n");
+ }
+}
+
+// The intput/output are all R indices (start from 1, NOT 0)
+void generateObsIndex(
+ //OUTPUT
+ int* obsIndex, // E.g., consider user i
+ int* start, // y[ obsIndex[ start[i]+(0:(num[i]-1)) ] ]
+ int* num, // are the observations of user i
+ //INPUT
+ const int* effIndex /* user or item */,
+ const int* nObs, const int* nEff,
+ //OTHER
+ const int* debug
+){
+ int *ind, i, j, cIndex, eIndex, oIndex;
+ ind = (int*)Calloc(*nEff, int);
+
+ for(i=0; i<*nObs; i++){
+ eIndex = effIndex[i];
+ if(*debug > 0){
+ CHK_R_INDEX(eIndex, *nEff);
+ }
+ num[R_VEC(eIndex)]++;
+ }
+
+ start[0] = 1; ind[0] = 1;
+ for(i=1; i<*nEff; i++){
+ start[i] = start[i-1]+num[i-1];
+ ind[i] = start[i];
+ }
+
+ for(i=0; i<*nObs; i++){
+ cIndex = R_VEC(effIndex[i]);
+ obsIndex[R_VEC(ind[cIndex])] = i+1;
+ ind[cIndex]++;
+ }
+
+ if(*debug > 0){
+ for(i=0; i<*nEff; i++){
+ if(ind[i] != start[i]+num[i]) error("logical error (level 1)");
+ if(*debug > 1){
+ for(j=0; j<num[i]; j++){
+ oIndex = obsIndex[R_VEC(start[i]+j)];
+ if(effIndex[R_VEC(oIndex)] != i+1)
+ error("logical error (level 2)");
+ }
+ }
+ }
+ }
+
+ Free(ind);
+}
+
+#define stepredn 0.2
+#define acctol 0.0001
+#define reltest 10.0
+#define _(A) A
+
+static double * vect(int n)
+{
+ return (double *)R_alloc(n, sizeof(double));
+}
+
+int nParam_debug = 0;
+double* gr_debug = NULL;
+
+void my_cgmin(int n, double *Bvec, double *X, double *Fmin,
+ optimfn fminfn, optimgr fmingr, int *fail,
+ double abstol, double intol, void *ex, int type, int trace,
+ int *fncount, int *grcount, int maxit)
+{
+ Rboolean accpoint;
+ double *c, *g, *t;
+ int count, cycle, cyclimit;
+ double f;
+ double G1, G2, G3, gradproj;
+ int funcount=0, gradcount=0, i;
+ double newstep, oldstep, setstep, steplength=1.0;
+ double tol;
+
+ if (maxit <= 0) {
+ *Fmin = fminfn(n, Bvec, ex);
+ *fncount = *grcount = 0;
+ *fail = FALSE;
+ return;
+ }
+ if (trace) {
+ Rprintf(" Conjugate gradients function minimizer\n");
+ switch (type) {
+ case 1: Rprintf("Method: Fletcher Reeves\n"); break;
+ case 2: Rprintf("Method: Polak Ribiere\n"); break;
+ case 3: Rprintf("Method: Beale Sorenson\n"); break;
+ default:
+ error(_("unknown 'type' in CG method of optim"));
+ }
+ }
+ c = vect(n); g = vect(n); t = vect(n);
+
+ // DEBUG
+ nParam_debug = n;
+ gr_debug = g;
+ // END
+
+ print_vector("INITIAL: ", gr_debug, nParam_debug);
+
+ setstep = 1.7;
+ *fail = 0;
+ cyclimit = n;
+ tol = intol * n * sqrt(intol);
+
+ if (trace) Rprintf("tolerance used in gradient test=%g\n", tol);
+ f = fminfn(n, Bvec, ex);
+ print_vector("AFTER fminfn IN cgmin: ", gr_debug, nParam_debug);
+ if (!R_FINITE(f)) {
+ error(_("Function cannot be evaluated at initial parameters"));
+ } else {
+ *Fmin = f;
+ funcount = 1;
+ gradcount = 0;
+ do {
+ for (i = 0; i < n; i++) {
+ t[i] = 0.0;
+ c[i] = 0.0;
+ }
+ cycle = 0;
+ oldstep = 1.0;
+ count = 0;
+ do {
+ cycle++;
+ if (trace) {
+ Rprintf("gradcount=%d, funcount=%d, Fmin=%f\n", gradcount, funcount, *Fmin);
+ Rprintf("parameters:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", Bvec[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ gradcount++;
+ if (gradcount > maxit) {
+ *fncount = funcount;
+ *grcount = gradcount;
+ *fail = 1;
+ return;
+ }
+ if (trace) {
+ Rprintf("BEFORE: gradient:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", g[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ fmingr(n, Bvec, g, ex);
+ if (trace) {
+ Rprintf("AFTER: gradient:\n");
+ for (i = 1; i <= n; i++) {
+ Rprintf("%10.5f ", g[i - 1]);
+ if (i / 7 * 7 == i && i < n)
+ Rprintf("\n");
+ }
+ Rprintf("\n");
+ }
+ G1 = 0.0;
+ G2 = 0.0;
+ for (i = 0; i < n; i++) {
+ X[i] = Bvec[i];
+ switch (type) {
+
+ case 1: /* Fletcher-Reeves */
+ G1 += g[i] * g[i];
+ G2 += c[i] * c[i];
+ break;
+
+ case 2: /* Polak-Ribiere */
+ G1 += g[i] * (g[i] - c[i]);
+ G2 += c[i] * c[i];
+ break;
+
+ case 3: /* Beale-Sorenson */
+ G1 += g[i] * (g[i] - c[i]);
+ G2 += t[i] * (g[i] - c[i]);
+ break;
+
+ default:
+ error(_("unknown type in CG method of optim"));
+ }
+ c[i] = g[i];
+ }
+ if (G1 > tol) {
+ if (G2 > 0.0)
+ G3 = G1 / G2;
+ else
+ G3 = 1.0;
+ gradproj = 0.0;
+ for (i = 0; i < n; i++) {
+ t[i] = t[i] * G3 - g[i];
+ gradproj += t[i] * g[i];
+ }
+ steplength = oldstep;
+
+ accpoint = FALSE;
+ do {
+ count = 0;
+ for (i = 0; i < n; i++) {
+ Bvec[i] = X[i] + steplength * t[i];
+ if (reltest + X[i] == reltest + Bvec[i])
+ count++;
+ }
+ if (count < n) { /* point is different */
+ f = fminfn(n, Bvec, ex);
+ funcount++;
+ accpoint = (R_FINITE(f) &&
+ f <= *Fmin + gradproj * steplength * acctol);
+
+ if (!accpoint) {
+ steplength *= stepredn;
+ if (trace){
+ Rprintf("* %f\n", f);
+ }
+ } else *Fmin = f; /* we improved, so update value */
+ }
+ } while (!(count == n || accpoint));
+ if (count < n) {
+ newstep = 2 * (f - *Fmin - gradproj * steplength);
+ if (newstep > 0) {
+ newstep = -(gradproj * steplength * steplength / newstep);
+ for (i = 0; i < n; i++)
+ Bvec[i] = X[i] + newstep * t[i];
+ *Fmin = f;
+ f = fminfn(n, Bvec, ex);
+ funcount++;
+ if (f < *Fmin) {
+ *Fmin = f;
+ if (trace) Rprintf("i< ");
+ } else { /* reset Bvec to match lowest point */
+ if (trace) Rprintf("i> ");
+ for (i = 0; i < n; i++)
+ Bvec[i] = X[i] + steplength * t[i];
+ }
+ }
+ }
+ }
+ oldstep = setstep * steplength;
+ if (oldstep > 1.0)
+ oldstep = 1.0;
+ } while ((count != n) && (G1 > tol) && (cycle != cyclimit));
+
+ } while ((cycle != 1) ||
+ ((count != n) && (G1 > tol) && *Fmin > abstol));
+
+ }
+ if (trace) {
+ Rprintf("Exiting from conjugate gradients minimizer\n");
+ Rprintf(" %d function evaluations used\n", funcount);
+ Rprintf(" %d gradient evaluations used\n", gradcount);
+ }
+ *fncount = funcount;
+ *grcount = gradcount;
+}
View
136 src/RLFM-ars-logistic/C/util.h
@@ -0,0 +1,136 @@
+/*
+ Copyright (c) 2011, Yahoo! Inc. All rights reserved.
+ Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
+
+ Author: Liang Zhang
+*/
+
+#include <R_ext/Applic.h>
+
+#define R_VEC(i) (i-1)
+#define R_MAT(i,j,nrow) (((j-1)*(nrow))+(i-1))
+#define R_3DA(i,j,k,nrow,ncol) ((((k-1)*(ncol))+(j-1))*(nrow) + (i-1))
+
+#define C_MAT(i,j,nrow) (((j)*(nrow))+(i))
+#define C_3DA(i,j,k,nrow,ncol) ((((k)*(ncol))+(j))*(nrow) + (i))
+
+#define CHK_C_INDEX(index, num) if(index < 0 || index >= num) error("index out of bound: index=%d, bound=%d", index, num)
+#define CHK_C_INDEX_2D(row_i, col_j, nrow, ncol) if(row_i < 0 || row_i >= nrow || col_j < 0 || col_j >= ncol) error("index out of bound: i=%d, j=%d nrow=%d ncol=%d", row_i, col_j, nrow, ncol)
+
+#define CHK_R_INDEX(index, num) if(index < 1 || index > num) error("index out of bound: index=%d, bound=%d", index, num)
+#define CHK_SYMMETRIC(A, nrow, i, j) for(i=0; i<nrow; i++) for(j=0; j<i; j++) if(abs(A[C_MAT(i,j,nrow)] - A[C_MAT(j,i,nrow)]) > 1e-10) error("A symmetric matrix is not symmetric: %f vs %f, diff=%e", A[C_MAT(i,j,nrow)], A[C_MAT(j,i,nrow)], A[C_MAT(i,j,nrow)] - A[C_MAT(j,i,nrow)])
+
+#define MAX(x,y) (x) > (y) ? (x) : (y)
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* xb/*nObs x 1*/;
+ double* g0w/*nUsers x 1*/; double* d0z/*nItems x 1*/;
+ double* Gw/*nUsers x nFactors*/; double* Dz/*nItems x nFactors*/;
+ double* var_y; double* var_alpha; double* var_beta; double* var_u; double* var_v;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* User_obsIndex; int* User_oiStart; int* User_oiNum;
+ int* Item_obsIndex; int* Item_oiStart; int* Item_oiNum;
+ int* debug;
+} RegParams;
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* xb/*nObs x 1*/;
+ double* g0w/*nUsers x 1*/; double* d0z/*nItems x 1*/;
+ double* Gw/*nUsers x nFactors*/; double* Dz/*nItems x nFactors*/;
+ double* var_y; double* var_alpha; double* var_beta; double* var_u; double* var_v;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* nVar_y; int* nVar_alpha; int* nVar_beta; int* nVar_u; int* nVar_v;
+ double* inv_var_u; double* inv_var_v;
+ int* User_obsIndex; int* User_oiStart; int* User_oiNum;
+ int* Item_obsIndex; int* Item_oiStart; int* Item_oiNum;
+ int* debug;
+} RegParams2;
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* x/*nObs x nJointFeatures*/;
+ double* w/*nUsers x nUserFeatures*/; double* z/*nItems x nItemFeatures*/;
+ double* lambda_b; double* lambda_g0; double* lambda_d0; double* lambda_G; double* lambda_D;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* nJointFeatures; int* nUserFeatures; int* nItemFeatures;
+ int* debug;
+} RegInputFeatureOnly;
+
+// Matrix inversion (only for a symmetric matrix)
+void sym_inv_byCholesky(double *A /* n x n matrix */, const int *n, const int *check_sym);
+
+// Matrix inversion (only for a symmetric matrix) 3-dim array
+void sym_inv3DA_byCholesky(
+ double *invA /* k x n x n matrix */, const double *A /* k x n x n matrix */,
+ const int *k, const int *n, double *temp /* n x n */, const int *check_sym
+);
+
+// Compute the eigenvalues and vectors of a symmetric matrix x
+void sym_eigen(const double* x, const int *nrow, double *eigen_val, double *eigen_vec);
+
+void sum_margin(
+ // OUTPUT
+ double *ans,
+ // INPUT
+ const double *A, const int *nrow, const int *ncol,
+ const int *side // side=1: Sum up each row and return a vector with length nrow
+ // side=2: Sum up each column and return a vector with length ncol
+);
+
+void center_array(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n // array length
+);
+
+void center_array_2d(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim // 1 for subtracting row mean, 2 for column
+);
+
+void center_array_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // array length
+ int * oiNum
+);
+
+void center_array_2d_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim, // 1 for subtracting row mean, 2 for column
+ int * oiNum
+);
+
+
+// The intput/output are all R indices (start from 1, NOT 0)
+void generateObsIndex(
+ //OUTPUT
+ int* obsIndex, // E.g., consider user i
+ int* start, // y[ obsIndex[ start[i]+(0:(num[i]-1)) ] ]
+ int* num, // are the observations of user i
+ //INPUT
+ const int* effIndex /* user or item */,
+ const int* nObs, const int* nEff,
+ //OTHER
+ const int* debug
+);
+
+void print_vector(char* prefix, double* vector, int length);
+void print_matrix(const char* prefix, const double* matrix, const int nrow, const int ncol);
+
+void my_cgmin(int n, double *Bvec, double *X, double *Fmin,
+ optimfn fminfn, optimgr fmingr, int *fail,
+ double abstol, double intol, void *ex, int type, int trace,
+ int *fncount, int *grcount, int maxit);
View
130 src/RLFM-ars-logistic/C/util.h~
@@ -0,0 +1,130 @@
+
+#include <R_ext/Applic.h>
+
+#define R_VEC(i) (i-1)
+#define R_MAT(i,j,nrow) (((j-1)*(nrow))+(i-1))
+#define R_3DA(i,j,k,nrow,ncol) ((((k-1)*(ncol))+(j-1))*(nrow) + (i-1))
+
+#define C_MAT(i,j,nrow) (((j)*(nrow))+(i))
+#define C_3DA(i,j,k,nrow,ncol) ((((k)*(ncol))+(j))*(nrow) + (i))
+
+#define CHK_C_INDEX(index, num) if(index < 0 || index >= num) error("index out of bound: index=%d, bound=%d", index, num)
+#define CHK_C_INDEX_2D(row_i, col_j, nrow, ncol) if(row_i < 0 || row_i >= nrow || col_j < 0 || col_j >= ncol) error("index out of bound: i=%d, j=%d nrow=%d ncol=%d", row_i, col_j, nrow, ncol)
+
+#define CHK_R_INDEX(index, num) if(index < 1 || index > num) error("index out of bound: index=%d, bound=%d", index, num)
+#define CHK_SYMMETRIC(A, nrow, i, j) for(i=0; i<nrow; i++) for(j=0; j<i; j++) if(abs(A[C_MAT(i,j,nrow)] - A[C_MAT(j,i,nrow)]) > 1e-10) error("A symmetric matrix is not symmetric: %f vs %f, diff=%e", A[C_MAT(i,j,nrow)], A[C_MAT(j,i,nrow)], A[C_MAT(i,j,nrow)] - A[C_MAT(j,i,nrow)])
+
+#define MAX(x,y) (x) > (y) ? (x) : (y)
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* xb/*nObs x 1*/;
+ double* g0w/*nUsers x 1*/; double* d0z/*nItems x 1*/;
+ double* Gw/*nUsers x nFactors*/; double* Dz/*nItems x nFactors*/;
+ double* var_y; double* var_alpha; double* var_beta; double* var_u; double* var_v;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* User_obsIndex; int* User_oiStart; int* User_oiNum;
+ int* Item_obsIndex; int* Item_oiStart; int* Item_oiNum;
+ int* debug;
+} RegParams;
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* xb/*nObs x 1*/;
+ double* g0w/*nUsers x 1*/; double* d0z/*nItems x 1*/;
+ double* Gw/*nUsers x nFactors*/; double* Dz/*nItems x nFactors*/;
+ double* var_y; double* var_alpha; double* var_beta; double* var_u; double* var_v;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* nVar_y; int* nVar_alpha; int* nVar_beta; int* nVar_u; int* nVar_v;
+ double* inv_var_u; double* inv_var_v;
+ int* User_obsIndex; int* User_oiStart; int* User_oiNum;
+ int* Item_obsIndex; int* Item_oiStart; int* Item_oiNum;
+ int* debug;
+} RegParams2;
+
+typedef struct {
+ int* user/*nObs x 1*/; int* item/*nObs x 1*/;
+ double* y/*nObs x 1*/; double* x/*nObs x nJointFeatures*/;
+ double* w/*nUsers x nUserFeatures*/; double* z/*nItems x nItemFeatures*/;
+ double* lambda_b; double* lambda_g0; double* lambda_d0; double* lambda_G; double* lambda_D;
+ int* nObs; int* nUsers; int* nItems; int* nFactors;
+ int* nJointFeatures; int* nUserFeatures; int* nItemFeatures;
+ int* debug;
+} RegInputFeatureOnly;
+
+// Matrix inversion (only for a symmetric matrix)
+void sym_inv_byCholesky(double *A /* n x n matrix */, const int *n, const int *check_sym);
+
+// Matrix inversion (only for a symmetric matrix) 3-dim array
+void sym_inv3DA_byCholesky(
+ double *invA /* k x n x n matrix */, const double *A /* k x n x n matrix */,
+ const int *k, const int *n, double *temp /* n x n */, const int *check_sym
+);
+
+// Compute the eigenvalues and vectors of a symmetric matrix x
+void sym_eigen(const double* x, const int *nrow, double *eigen_val, double *eigen_vec);
+
+void sum_margin(
+ // OUTPUT
+ double *ans,
+ // INPUT
+ const double *A, const int *nrow, const int *ncol,
+ const int *side // side=1: Sum up each row and return a vector with length nrow
+ // side=2: Sum up each column and return a vector with length ncol
+);
+
+void center_array(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n // array length
+);
+
+void center_array_2d(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim // 1 for subtracting row mean, 2 for column
+);
+
+void center_array_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // array length
+ int * oiNum
+);
+
+void center_array_2d_online(
+ // Array to be centered
+ double *x,
+ // INPUT
+ int n, // number of rows
+ int m, // number of columns
+ int dim, // 1 for subtracting row mean, 2 for column
+ int * oiNum
+);
+
+
+// The intput/output are all R indices (start from 1, NOT 0)
+void generateObsIndex(
+ //OUTPUT
+ int* obsIndex, // E.g., consider user i
+ int* start, // y[ obsIndex[ start[i]+(0:(num[i]-1)) ] ]
+ int* num, // are the observations of user i
+ //INPUT
+ const int* effIndex /* user or item */,
+ const int* nObs, const int* nEff,
+ //OTHER
+ const int* debug
+);
+
+void print_vector(char* prefix, double* vector, int length);
+void print_matrix(const char* prefix, const double* matrix, const int nrow, const int ncol);
+
+void my_cgmin(int n, double *Bvec, double *X, double *Fmin,
+ optimfn fminfn, optimgr fmingr, int *fail,
+ double abstol, double intol, void *ex, int type, int trace,
+ int *fncount, int *grcount, int maxit);
View
59 src/RLFM-ars-logistic/R/bounded.logistic.R
@@ -0,0 +1,59 @@
+### Copyright (c) 2011, Yahoo! Inc. All rights reserved.
+### Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
+###
+### Author: Liang Zhang
+
+# doesn't work yet
+logexp <- function(eta){
+ pos <- eta > 0
+ ans <- rep(NA,length(eta))
+ ans[pos] <- eta[pos] + log(1 + exp(-eta[pos]))
+ ans[!pos] <- log(1 + exp(eta[!pos]))
+ ans
+}
+
+getg <- function(b,x,y,offset,spline,rho,penalize.firstcolumn){
+ eta <- x%*%b+offset
+ p = 2*spline/(1+exp(-2*(1-spline)*eta));
+ p[eta>=0] = 2*spline - 1 + 2*(1-spline)/(1+exp(-2*spline*eta[eta>=0]));
+ s = - sum(y*log(p)+(1-y)*log(1-p));
+ s = s + rho*sum(b^2)
+ if (!penalize.firstcolumn) s = s - rho*b[1]*b[1]
+ s
+}
+
+gradg <- function(b,x,y,offset,spline,rho,penalize.firstcolumn){
+ eta <- x%*%b+offset
+ q <- y*2*(1-spline)/(1+exp(2*(1-spline)*eta)) + (y-1)*4*spline*(1-spline)/(1+exp(2*(1-spline)*eta))/(1-2*spline+exp(-2*(1-spline)*eta));
+ q[eta>=0] <- y[eta>=0]*4*spline*(1-spline)/(1+(2*spline-1)*exp(-2*spline*eta[eta>=0]))/(1+exp(2*spline*eta[eta>=0]))+(y[eta>=0]-1)*spline/(1-spline)/(1+exp(-2*spline*eta[eta>=0]));
+ grad = -t(x)%*%q
+ t = 2*rho*b;
+ if (!penalize.firstcolumn) t[1] = 0;
+ grad + t
+ }
+
+# y is the binary response
+# x is the design matrix including intercept
+# beta saves the initial values of the coefficients
+# offset is the offset of the observation, default NULL (no offset), length = length(y)
+# spline is the spline knots for logistic regression, default 0.5, 0<spline<1
+# rho is the L2 penalty parameter, the Loss function = - Loglik + rho*sum(beta^2), default 0 (no penalty)
+# penalize.firstcolumn indicates whether we want to penalize the first column of beta (usually intercept). default TRUE.
+# lower is the lower bound of each coefficient including intercept, If a scalar, lower = rep(lower,length(y))
+# upper is the upper bound of each coefficient including intercept. Can also be a scalar
+# maxit is the maximum number of iterations for lbfgs
+bounded.logistic<-function(y,x,beta,offset=NULL,spline=0.5,rho=0,penalize.firstcolumn=T,lower=-Inf,upper=Inf,maxit=1000)
+{
+ if (!is.vector(y)) stop("The response y has to be a vector!");
+ if (!is.matrix(x)) stop("The design matrix x has to be a matrix!");
+ if (!is.vector(beta)) stop("The initial values of beta has to be a vector!");
+ if (ncol(x)!=length(beta)) stop("ncol(x)!=length(beta)");
+ if (is.null(offset)) {
+ offset = rep(0,length(y));
+ } else {
+ if (length(offset)!=length(y)) stop("length(offset)!=length(y)");
+ }
+ if (rho<0) stop ("rho<0");
+ optim(par=beta,fn=getg,gr=gradg,x=x,y=y,offset=offset,spline=spline,rho=rho,penalize.firstcolumn=penalize.firstcolumn,method="L-BFGS-B",lower=lower,upper=upper, control=list(maxit=maxit))
+}
+
View
1,252 src/RLFM-ars-logistic/R/c_funcs.R
@@ -0,0 +1,1252 @@
+### Copyright (c) 2011, Yahoo! Inc. All rights reserved.
+### Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
+###
+### Author: Liang Zhang
+
+###
+### Inversion of a symmetric metrix using Cholesky
+### This function is implemented to make sure the C implementation generates exactly
+### the same output as the R implementation
+###
+sym_inv.cholesky <- function(x){
+ if(!is.double(x)) stop("x is not a double matrix");
+ if(!is.matrix(x)) stop("x is not a double matrix");
+ if(nrow(x) != ncol(x)) stop("x is not a square matrix");
+ n = as.integer(nrow(x));
+ ans = .C("sym_inv_byCholesky", x, n, as.integer(1), DUP=TRUE);
+ return(ans[[1]]);
+}
+
+###
+### Eigen decomposition of a symmetric metrix
+### This function is implemented to make sure the C implementation generates exactly
+### the same output as the R implementation
+###
+sym_eigen <- function(x){
+ if(!is.matrix(x)) stop("x is not matrix");
+ if(nrow(x) != ncol(x)) stop("nrow(x) != ncol(x)");
+ if(!is.double(x)) stop("x is not a double-precision matrix");
+ output = list(values = rep(0.0, nrow(x)),
+ vectors = matrix(0.0, nrow=nrow(x), ncol=ncol(x)));
+ .C("sym_eigen", x, nrow(x), output$values, output$vectors, DUP=FALSE);
+ return(output);
+}
+
+###
+### Sum up each row (or column) of a matrix
+### side=1: Sum up each row and return a vector with length nrow
+### side=2: Sum up each column and return a vector with length ncol
+###
+sum_margin <- function(A, side){
+ if(!is.matrix(A)) stop("'A' is not a matrix");
+ if(!is.double(A)) stop("'A' is not double");
+ if(side == 1) out = double(nrow(A))
+ else if(side == 2) out = double(ncol(A))
+ else stop("Unknown side=",side," (side=1 for rows, side=2 for columns)");
+ ans = .C("sum_margin", out, A, as.integer(nrow(A)), as.integer(ncol(A)), as.integer(side), DUP=FALSE);
+ return(out);
+}
+
+###
+### Draw a random sample from a multivariate normal distribtuion
+### This function is implemented to make sure the C implementation generates exactly
+### the same output as the R implementation
+###
+my_rmvnorm <- function(n, mu, Sigma, debug=10, tol=1e-8){
+ p <- length(mu)
+ if (!all(dim(Sigma) == c(p, p)))
+ stop("incompatible arguments")
+ eS <- sym_eigen(solve(Sigma))
+ ev <- 1/eS$values
+
+ if(debug >= 3){
+ if(max(abs(eS$vectors %*% diag(ev, p) %*% t(eS$vectors) - Sigma)) > tol * abs(ev[1]))
+ stop("sym_eigen(Sigma) seems to have some problems!!");
+ }
+
+ if (!all(ev >= -tol * abs(ev[1])))
+ stop("'Sigma' is not positive definite")
+ Z <- matrix(rnorm(p * n), n)
+
+ # cat("temp:\n");
+ # print(eS$vectors %*% diag(sqrt(pmax(ev, 0)), p));
+ # cat("rnd: ");
+ # print(drop(Z));
+
+ X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(Z)
+ nm <- names(mu)
+ if (is.null(nm) && !is.null(dn <- dimnames(Sigma)))
+ nm <- dn[[1]]
+ dimnames(X) <- list(nm, NULL)
+ if (n == 1)
+ return(drop(X))
+ else return(t(X))
+}
+
+###
+### The E-Step for the Monte-Carlo EM algorithm (version 1)
+###
+MC_EStep.C <- function(
+ nSamples,
+ user, item, y, x, w, z,
+ alpha, beta, u, v,
+ b, g0, G, d0, D,
+ var_y, var_alpha, var_beta, var_u, var_v=1, debug=0
+){
+ if(debug >= 1) check.input(user, item, y, x, w, z, alpha, beta, u, v, b, g0, G, d0, D, var_y, var_alpha, var_beta, var_u, var_v);
+ nObs = as.integer(length(y));
+ nUsers = as.integer(length(alpha));
+ nItems = as.integer(length(beta));
+ nFactors = as.integer(ncol(u));
+
+ xb = x %*% b;
+ g0w = w %*% g0;
+ d0z = z %*% d0;
+ Gw = w %*% G;
+ Dz = z %*% D;
+
+ output = list(
+ o.mean = rep(0.0, nObs),
+ o.sumvar = 0.0,
+ alpha.mean = rep(0.0, nUsers),
+ alpha.sumvar = 0.0,
+ beta.mean = rep(0.0, nItems),
+ beta.sumvar = 0.0,
+ u.mean = matrix(0.0, nrow=nUsers, ncol=nFactors),
+ u.sumvar = 0.0,
+ v.mean = matrix(0.0, nrow=nItems, ncol=nFactors),
+ v.sumvar = 0.0,
+ pred.y.square = rep(0.0, nObs)
+ );
+
+ if(!is.integer(user)) stop("!is.integer(user)");
+ if(!is.integer(item)) stop("!is.integer(item)");
+ if(!is.double(y)) stop("!is.double(y)");
+ if(!is.double(xb)) stop("!is.double(xb)");
+ if(!is.double(g0w)) stop("!is.double(g0w)");
+ if(!is.double(d0z)) stop("!is.double(d0z)");
+ if(!is.double(Gw)) stop("!is.double(Gw)");
+ if(!is.double(Dz)) stop("!is.double(Dz)");
+ if(!is.double(alpha)) stop("!is.double(alpha)");
+ if(!is.double(beta)) stop("!is.double(beta)");
+ if(!is.double(u)) stop("!is.double(u)");
+ if(!is.double(v)) stop("!is.double(v)");
+ if(!is.double(var_y)) stop("!is.double(var_y)");
+ if(!is.double(var_alpha)) stop("!is.double(var_alpha)");
+ if(!is.double(var_beta)) stop("!is.double(var_beta)");
+ if(!is.double(var_u)) stop("!is.double(var_u)");
+ if(!is.double(var_v)) stop("!is.double(var_v)");
+
+ .C("MCEM_EStep",
+ # OUTPUT
+ output$o.mean, output$o.sumvar,
+ output$alpha.mean, output$alpha.sumvar,
+ output$beta.mean, output$beta.sumvar,
+ output$u.mean, output$u.sumvar,
+ output$v.mean, output$v.sumvar,
+ output$pred.y.square,
+ # INPUT
+ as.integer(nSamples),
+ user, item, y, xb, g0w, d0z, Gw, Dz,
+ alpha, beta, u, v,
+ var_y, var_alpha, var_beta, var_u, var_v,
+ nObs, nUsers, nItems, nFactors, as.integer(length(var_y)),
+ # OTHER
+ as.integer(debug),
+ DUP=FALSE
+ );
+ return(output);
+}
+
+###
+### The E-Step for the Monte-Carlo EM algorithm (version 1)
+###
+MC_EStep_logistic.C <- function(
+ nSamples, nBurnin,
+ user, item, y, x, w, z,
+ alpha, beta, u, v,
+ b, g0, G, d0, D,
+ var_alpha, var_beta, var_u, var_v=1, debug=0
+){
+ if(debug >= 1) check.input.logistic(user, item, y, x, w, z, alpha, beta, u, v, b, g0, G, d0, D, var_alpha, var_beta, var_u, var_v);
+ nObs = as.integer(length(y));
+ nUsers = as.integer(length(alpha));
+ nItems = as.integer(length(beta));
+ nFactors = as.integer(ncol(u));
+
+ xb = x %*% b;
+ g0w = w %*% g0;
+ d0z = z %*% d0;
+ Gw = w %*% G;
+ Dz = z %*% D;
+
+ output = list(
+ o.mean = rep(0.0, nObs),
+ alpha.mean = rep(0.0, nUsers),
+ alpha.sumvar = 0.0,
+ beta.mean = rep(0.0, nItems),
+ beta.sumvar = 0.0,
+ u.mean = matrix(0.0, nrow=nUsers, ncol=nFactors),
+ u.sumvar = 0.0,
+ v.mean = matrix(0.0, nrow=nItems, ncol=nFactors),
+ v.sumvar = 0.0,
+ acceptrate.maineff = 0.0,
+ acceptrate.fact = 0.0
+ );
+
+ if(!is.integer(user)) stop("!is.integer(user)");
+ if(!is.integer(item)) stop("!is.integer(item)");
+ if(!is.double(y)) stop("!is.double(y)");
+ if(!is.double(xb)) stop("!is.double(xb)");
+ if(!is.double(g0w)) stop("!is.double(g0w)");
+ if(!is.double(d0z)) stop("!is.double(d0z)");
+ if(!is.double(Gw)) stop("!is.double(Gw)");
+ if(!is.double(Dz)) stop("!is.double(Dz)");
+ if(!is.double(alpha)) stop("!is.double(alpha)");
+ if(!is.double(beta)) stop("!is.double(beta)");
+ if(!is.double(u)) stop("!is.double(u)");
+ if(!is.double(v)) stop("!is.double(v)");
+ if(!is.double(var_alpha)) stop("!is.double(var_alpha)");
+ if(!is.double(var_beta)) stop("!is.double(var_beta)");
+ if(!is.double(var_u)) stop("!is.double(var_u)");
+ if(!is.double(var_v)) stop("!is.double(var_v)");
+
+ .C("MCEM_EStep_logistic",
+ # OUTPUT
+ output$o.mean,
+ output$alpha.mean, output$alpha.sumvar,
+ output$beta.mean, output$beta.sumvar,
+ output$u.mean, output$u.sumvar,
+ output$v.mean, output$v.sumvar,
+ output$acceptrate.maineff, output$acceptrate.fact,
+ # INPUT
+ as.integer(nSamples), as.integer(nBurnin),
+ user, item, y, xb, g0w, d0z, Gw, Dz,
+ alpha, beta, u, v,
+ var_alpha, var_beta, var_u, var_v,
+ nObs, nUsers, nItems, nFactors,
+ # OTHER
+ as.integer(debug),
+ DUP=FALSE
+ );
+ return(output);
+}
+
+###
+### The E-Step for the Monte-Carlo EM algorithm (version 2)
+###
+### Options:
+###
+### outputPerUserVar: [TRUE/FALSE] Whether to output per-user variance-covariance matrix
+### outputPerItemVar: [TRUE/FALSE] Whether to output per-item variance-covariance matrix
+###
+### isOldUser[i]: [TRUE/FALSE] Whehter the ith user is an old user; default: all FALSE
+### For old users, the input alpha[i], u[i,] will be used instead of g0w[i] and Gw[i]
+### isOldItem[j]: [TRUE/FALSE] Whehter the jth item is an old item; default: all FALSE
+### For old items, the input beta[j], v[j,] will be used instead of d0z[j] and Dz[j]
+###
+### var_y (1x1 or nObs x 1) specifies global or per-observation variance
+### var_alpha (1x1 or nUsers x 1) specifies global or per-user variance
+### var_beta (1x1 or nItems x 1) specifies global or per_item variance
+### var_u (1x1 or nUsers x nFactors x nFactors) specifies global variance or per-user variance-covariance matrix
+### var_v (1x1 or nItems x nFactors x nFactors) specifies global variance or per-item variance-covariance matrix
+###
+### Output:
+###
+### output$o.mean[k] = E[ o[k] = alpha[user[k]] + beta[item[k]] + t(u[user[k]]) %*% v[item[k]] ]
+### output$o.sumvar = sum_k Var(o[k])
+###
+### output$alpha.mean[i] = E[ alpha[i] ]
+### output$alpha.sumvar = sum_i Var( alpha[i] )
+### output$alpha.var[i] = Var( alpha[i] ), which is available only when outputPerUserVar=TRUE
+###
+### output$beta.mean[i] = E[ beta[i] ]
+### output$beta.sumvar = sum_i Var( beta[i] )
+### output$beta.var[i] = Var( beta[i] ), which is available only when outputPerItemVar=TRUE
+###
+### output$u.mean[i,] = E[ u[i] ]; u.mean is nUsers x nFactors
+### output$u.sumvar = sum_i Var( u[i] )
+### output$u.var[i,,] = VarCov( u[i] ), which is available only when outputPerUserVar=TRUE
+### u.var is nUsers x nFactors x nFactors
+###
+### output$v.mean[i,] = E[ v[i] ]; v.mean is nItems x nFactors
+### output$v.sumvar = sum_i Var( v[i] )
+### output$v.var[i,,] = VarCov( v[i] ), which is available only when outputPerItemVar=TRUE
+### v.var is nItems x nFactors x nFactors
+###
+MC_EStep2.C <- function(
+ nSamples,
+ user, item, y, x, w, z,
+ alpha, beta, u, v,
+ b, g0, G, d0, D,
+ var_y, var_alpha, var_beta, var_u, var_v=1,
+ outputPerUserVar=FALSE,
+ outputPerItemVar=FALSE,
+ isOldUser=NULL,
+ isOldItem=NULL,
+ debug=0
+){
+ if(debug >= 1) check.input(user, item, y, x, w, z, alpha, beta, u, v, b, g0, G, d0, D, var_y, var_alpha, var_beta, var_u, var_v, version=2);
+ nObs = as.integer(length(y));
+ nUsers = as.integer(length(alpha));
+ nItems = as.integer(length(beta));
+ nFactors = as.integer(ncol(u));
+
+ nVar_y = as.integer(length(var_y));
+ nVar_alpha = as.integer(length(var_alpha));
+ nVar_beta = as.integer(length(var_beta));
+ nVar_u = as.integer(1); if(length(var_u) > 1) nVar_u = as.integer(dim(var_u)[1]);
+ nVar_v = as.integer(1); if(length(var_v) > 1) nVar_v = as.integer(dim(var_v)[1]);
+
+ xb = x %*% b;
+ if(is.null(isOldUser)){
+ g0w = w %*% g0;
+ Gw = w %*% G;
+ }else{
+ if(length(isOldUser) != nUsers) stop("length(isOldUser) != nUsers");
+ w.new = w[!isOldUser,,drop=FALSE];
+ g0w = alpha;
+ Gw = u;
+ g0w[!isOldUser] = w.new %*% g0;
+ Gw[!isOldUser,] = w.new %*% G;
+ }
+
+ if(is.null(isOldItem)){
+ d0z = z %*% d0;
+ Dz = z %*% D;
+ }else{
+ if(length(isOldItem) != nItems) stop("length(isOldItem) != nItems");
+ z.new = z[!isOldItem,,drop=FALSE];
+ d0z = beta;
+ Dz = v;
+ d0z[!isOldItem] = z.new %*% d0;
+ Dz[!isOldItem,] = z.new %*% D;
+ }
+
+ output = list(
+ o.mean = rep(0.0, nObs),
+ o.sumvar = 0.0,
+ alpha.mean = rep(0.0, nUsers),
+ alpha.sumvar = 0.0,
+ beta.mean = rep(0.0, nItems),
+ beta.sumvar = 0.0,
+ u.mean = matrix(0.0, nrow=nUsers, ncol=nFactors),
+ u.sumvar = 0.0,
+ v.mean = matrix(0.0, nrow=nItems, ncol=nFactors),
+ v.sumvar = 0.0
+ );
+
+ if(outputPerUserVar){
+ perUserVar = as.integer(1);
+ output$alpha.var = rep(0.0, nUsers);
+ output$u.var = array(0.0, dim=c(nUsers, nFactors, nFactors));
+ }else{
+ perUserVar = as.integer(0);
+ }
+
+ if(outputPerItemVar){
+ perItemVar = as.integer(1);
+ output$beta.var = rep(0.0, nItems);
+ output$v.var = array(0.0, dim=c(nItems, nFactors, nFactors));
+ }else{
+ perItemVar = as.integer(0);
+ }
+
+ if(!is.integer(user)) stop("!is.integer(user)");
+ if(!is.integer(item)) stop("!is.integer(item)");
+ if(!is.double(y)) stop("!is.double(y)");
+ if(!is.double(xb)) stop("!is.double(xb)");
+ if(!is.double(g0w)) stop("!is.double(g0w)");
+ if(!is.double(d0z)) stop("!is.double(d0z)");