/src/C/util.h

  1/*
  2	Copyright (c) 2011, Yahoo! Inc.  All rights reserved.
  3	Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
  4
  5    Author: Bee-Chung Chen
  6*/
  7
  8
  9#include <R_ext/Applic.h>
 10
 11#define R_VEC(i) (i-1)
 12#define R_MAT(i,j,nrow) (((j-1)*(nrow))+(i-1))
 13#define R_3DA(i,j,k,nrow,ncol) ((((k-1)*(ncol))+(j-1))*(nrow) + (i-1))
 14#define R_4DA(i,j,k,m,dim1,dim2,dim3) (((((m-1)*(dim3)+(k-1))*(dim2))+(j-1))*(dim1) + (i-1))
 15
 16#define C_MAT(i,j,nrow) (((j)*(nrow))+(i))
 17#define C_3DA(i,j,k,nrow,ncol) ((((k)*(ncol))+(j))*(nrow) + (i))
 18#define C_4DA(i,j,k,m,dim1,dim2,dim3) (((((m)*(dim3)+(k))*(dim2))+(j))*(dim1) + (i))
 19
 20#define CHK_C_INDEX(index, num) if(index < 0 || index >= num) error("index out of bound: index=%d, bound=%d (file: %s, line: %d)", index, num, __FILE__, __LINE__)
 21#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 (file: %s, line: %d)", row_i, col_j, nrow, ncol, __FILE__, __LINE__)
 22
 23#define CHK_R_INDEX(index, num) if(index < 1 || index > num) error("index out of bound: index=%d, bound=%d (file: %s, line: %d)", index, num, __FILE__, __LINE__)
 24#define CHK_SYMMETRIC(A, nrow) for(int iii=0; iii<nrow; iii++) for(int jjj=0; jjj<iii; jjj++) if(fabs(A[C_MAT(iii,jjj,nrow)] - A[C_MAT(jjj,iii,nrow)]) > 1e-10) error("A symmetric matrix is not symmetric: %f vs %f, diff=%e (file: %s, line: %d)", A[C_MAT(iii,jjj,nrow)], A[C_MAT(jjj,iii,nrow)], A[C_MAT(iii,jjj,nrow)] - A[C_MAT(jjj,iii,nrow)], __FILE__, __LINE__)
 25
 26#define CHK_SAME_NUMBER(msg, x, y) if((x != 0 || y != 0) && (fabs(x-y) / fmax(fabs(x), fabs(y)) > 1e-8)) error("Error: %s. The two number should be the same: %f vs %f (file: %s, line: %d)", msg, x, y, __FILE__, __LINE__)
 27
 28#define MY_REAL(x) (TYPEOF(x) == NILSXP ? NULL : REAL(x))
 29#define MY_INTEGER(x) (TYPEOF(x) == NILSXP ? NULL : INTEGER(x))
 30
 31#define CHK_MAT_SYM(msg, A, nrow) for(int CHK_MAT_SYM_i=0; CHK_MAT_SYM_i<nrow; CHK_MAT_SYM_i++) for(int CHK_MAT_SYM_j=0; CHK_MAT_SYM_j<CHK_MAT_SYM_i; CHK_MAT_SYM_j++) CHK_SAME_NUMBER(msg, A[C_MAT(CHK_MAT_SYM_i,CHK_MAT_SYM_j,nrow)], A[C_MAT(CHK_MAT_SYM_j,CHK_MAT_SYM_i,nrow)])
 32
 33#define MAX(x,y) ((x) > (y) ? (x) : (y))
 34#define SQR(x) ((x) * (x))
 35
 36#define STOP(msg) error("Error at %s on line %d: %s\n", __FILE__, __LINE__, msg)
 37#define STOP1(msg,x) error2(__FILE__, __LINE__, msg, x)
 38#define STOP2(msg,x1,x2) error2(__FILE__, __LINE__, msg, x1, x2)
 39#define STOP3(msg,x1,x2,x3) error2(__FILE__, __LINE__, msg, x1, x2, x3)
 40#define STOP4(msg,x1,x2,x3,x4) error2(__FILE__, __LINE__, msg, x1, x2, x3, x4)
 41#define DIE_HERE error("Error in file: %s, at line: %d", __FILE__, __LINE__)
 42
 43void error2(const char *filename, int lineno, const char *fmt, ...);
 44
 45// Matrix inversion (only for a symmetric matrix)
 46void sym_inv_byCholesky(double *A /* n x n matrix */, const int *n, const int *check_sym);
 47
 48// Matrix inversion (only for a symmetric matrix) 3-dim array
 49void sym_inv3DA_byCholesky(
 50    double *invA /* k x n x n matrix */, const double *A /* k x n x n matrix */,
 51    const int *k, const int *n, double *temp /* n x n */, const int *check_sym
 52);
 53
 54// Compute the eigenvalues and vectors of a symmetric matrix x
 55void sym_eigen(const double* x, const int *nrow, double *eigen_val, double *eigen_vec);
 56
 57// Compute the eigenvalues and vectors of a symmetric matrix x
 58void sym_eigen2(const double* x, const int *nrow, double *eigen_val, double *eigen_vec, double *workspace, const int *workspace_size, const int *check_sym);
 59int sym_eigen2_workspace_size(const int *nrow);
 60
 61void sum_margin(
 62    // OUTPUT
 63    double *ans,
 64    // INPUT
 65    const double *A, const int *nrow, const int *ncol,
 66    const int *side // side=1: Sum up each row and return a vector with length nrow
 67                    // side=2: Sum up each column and return a vector with length ncol
 68);
 69
 70// The intput/output are all R indices (start from 1, NOT 0)
 71void generateObsIndex(
 72    //OUTPUT
 73    int* obsIndex, // E.g., consider user i (i starts from 0)
 74    int* start,    //  y[ obsIndex[ start[i]+(0:(num[i]-1)) ] ]
 75    int* num,      //  are the observations of user i
 76    //INPUT
 77    const int* effIndex /* user or item */,
 78    const int* nObs, const int* nEff,
 79    //OTHER
 80    const int* debug
 81);
 82
 83
 84void normalizeToSumUpToOne2(double *output, const double *input, const int length);
 85
 86void normalizeToSumUpToOne(double *vector, const int length);
 87
 88void indexWithQuantities(int *output, int *vector, const int *length);
 89
 90void print_vector(const char* prefix, const double* vector, const int length);
 91void print_intVector(const char* prefix, const int* vector, const int length);
 92void print_matrix(const char* prefix, const double* matrix, const int nrow, const int ncol);
 93
 94void compute_uBv_dense(
 95    // Output
 96    double *score, // nObs x 1
 97    // Input
 98    const int *userIndex, const int *itemIndex, // Both nObs x 1; index starts from 1 (not 0)
 99    const double *u, // nUsers x nUserFactors
100    const double *B, // nUserFeatures x nItemFeatures
101    const double *v, // nItems x nItemFactors
102    const int *nObs, const int *nUsers, const int *nItems,
103    const int *nUserFactors, const int *nItemFactors
104);
105
106double compute_szuBv_single_dense(
107	int i, // user i (start from 1)
108	int j, // item j (start from 1)
109	const double *s, // nUsers x nUserClusters
110	const double *z, // nItems x nItemClusters
111    const double *u, // nUsers x nFactorsPerUser
112	const double *B, // nUserClusters x nItemClusters x nFeaturesPerUser x nFactorsPerItem
113    const double *v, // nItems x nFactorsPerItem
114    const int nUsers, const int nItems,
115    const int nUserClusters, const int nItemClusters,
116    const int nFactorsPerUser, const int nFactorsPerItem
117);
118
119void compute_szuBv_dense(
120	double *ans, // nObs x 1
121	const int* userIndex, // nObs x 1 (start from 1)
122	const int* itemIndex, // nObs x 1 (start from 1)
123	const double *s, // nUsers x nUserClusters
124	const double *z, // nItems x nItemClusters
125	const double *u, // nUsers x nFactorsPerUser
126	const double *B, // nUserClusters x nItemClusters x nFeaturesPerUser x nFactorsPerItem
127	const double *v, // nItems x nFactorsPerItem
128	const int *nObs, const int *nUsers, const int *nItems,
129	const int *nUserClusters, const int *nItemClusters,
130	const int *nFactorsPerUser, const int *nFactorsPerItem,
131	const int *debug
132);
133
134int rdiscrete(double* probabilities, int nOutcomes, double *rnd_out);
135
136void bayesian_gaussian_regression(
137	double *postMean,     // (output)   nFeatures x 1
138	double *postVar,      // (output)   nFeatures x nFeatures
139	const double *y,      // (response) nObs x 1
140	const double *X,      // (features) nObs x nFeatures
141	const double *offset, // nObs x 1
142	const double *prior_mean, // nFeatures x 1
143	const double *prior_var,  // nFeatures x nFeatures
144	const double *var_y,      // nObs x 1
145	const int *nObs, const int *nFeatures, const int *nOffset,
146	const int *debug, const int *verbose
147);
148
149void computeMeanSumvar(
150    // OUTPUT
151    double *mean   /* length x 1 */,
152    double *sumvar /* 1x1 */,
153    // INPUT
154    const double *sum /* length x 1 */,
155    const double *sos /* sum of squares: length x 1 */,
156    const int length, const int nSamples
157);
158
159void computeMeanVar(
160    // OUTPUT
161    double *mean   /* nLevels x nFactors */,
162    double *sumvar /* 1x1 */,
163    double *outputVar /*IN:sum-of-products; OUT:variance; nLevels x nFactors x nFactors*/,
164    // INPUT
165    const double *sum /* nLevels x nFactors */,
166    const int nLevels, const int nFactors, const int nSamples
167);
168
169void computeCov(
170    // OUTPUT
171    double *outputCov /*IN:sum-of-products; OUT:covariance; size: num1 x num2*/,
172    // INPUT
173    const double *mean1, const double *mean2,
174    const int num1, const int num2, const int nSamples
175);
176
177/**
178 *  output[m] = u[src_id[m], , src_ctx[m]]' v[dst_id[m, , dst_ctx[m]]
179 */
180void computeMultiResponseUV(
181	// OUTPUT
182	double *output,
183	// INPUT
184	const double *u, // nSrcNodes x nFactors x nSrcContexts
185	const double *v, // nDstNodes x nFactors x nDstContexts
186	const int *src_id,  const int *dst_id, // nObs x 1
187	const int *src_ctx, const int *dst_ctx, // nObs x 1
188	const int *nObs, const int *nFactors,
189	const int *nSrcNodes, const int *nSrcContexts,
190	const int *nDstNodes, const int *nDstContexts,
191	const int *debug
192);
193
194
195/**
196 * Draw a multivariate Gaussian vector
197 * option == 0 : A    is the variance-covariance matrix
198 * option == 1 : A^-1 is the variance-covariance matrix
199 * option == 2 : eigen_val and eigen_vec are the eigen value and vector of
200 *               the variance-covariance matrix (A is not used)
201 * In any case, the output eigen_val and eigen_vec are the eigen value and vector
202 * of the variance-covariance matrix.
203 *
204 * GetRNGstate() and PutRNGstate() are NOT called inside draw_multivar_gaussian().
205 */
206void draw_multivar_gaussian(
207	// OUTPUT
208	double* out /* nDim x 1 */,
209	// INPUT/OUTPUT
210	double* eigen_val /* nDim x 1 */,
211	double* eigen_vec /* nDim x nDim */,
212	// INPUT
213	const double* mean /* nDim x 1*/,
214	const double* A    /* nDim x nDim  or  NULL */,
215	const int* nDim,
216	const int* option, const int* debug,
217	// WORK/TEMP SPACE (see workspace_size_for_draw_multivar_gaussian(nDim))
218	double *workspace /* workspace_size x 1 */,
219	const int* workspace_size /* 0 or at least 3*nDim or call workspace_size_for_draw_multivar_gaussian(nDim) */,
220	double *temp1 /* nDim x 1 */,
221	double *temp2 /* nDim x nDim */
222);
223int workspace_size_for_draw_multivar_gaussian(int nDim);