Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

remove wrong files

  • Loading branch information...
commit 3815cbb311da8819b686661ce7007a7cb62e0f7a 1 parent 622fa89
Liang Zhang authored
View
399 src/RLFM-ars-logistic/C/MCEM_EStep_logistic.c~
@@ -1,399 +0,0 @@
-/*
-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
500 src/RLFM-ars-logistic/C/util.c~
@@ -1,500 +0,0 @@
-
-
-#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
130 src/RLFM-ars-logistic/C/util.h~
@@ -1,130 +0,0 @@
-
-#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);
Please sign in to comment.
Something went wrong with that request. Please try again.