/** * @file utilities.h * @brief header files for utility functions * @author Wayne Zhang */ #ifndef CPLM_UTILS_H #define CPLM_UTILS_H /* inline functions */ /** * cumulative sum of a vector of double * * @param x double vector to be summed * @param n number of elements * * @return cumulative sum */ static R_INLINE double dcumsum(double *x, int n){ double sm = 0.0 ; for (int i = 0; i < n; i++) sm += x[i] ; return sm ; } /** * get the maximum value of a double vector * * @param x a double vector * @param n length of the vector * * @return the maximum value */ static R_INLINE double dmax (double *x, int n){ double s = x[0] ; for (int i = 1; i < n; i++) if (x[i] > s) s = x[i] ; return s ; } /** * get the minmimum value of a double vector * * @param x a double vector * @param n length of the vector * * @return the minimum value */ static R_INLINE double dmin(double *x, int n){ double s = x[0]; for (int i = 1; i < n; i++){ if (x[i] < s) s = x[i]; } return s; } /** * get the maximum value of an int vector * * @param x an int vector * @param n length of the vector * * @return the maximum value */ static R_INLINE int imax (int *x, int n){ int s = x[0] ; for (int i = 1; i < n; i++) if (x[i] > s) s = x[i] ; return s ; } /** * Compute sample mean * * @param n number of samples * @param x samples in long vector * * @return mean */ static R_INLINE double mean(double *x, int n){ return dcumsum(x, n) / n ; } /** * inverse link function * * @param eta linear predictors * @param lp link power * * @return mean */ static R_INLINE double link_inv(double eta, double lp){ return (lp == 0) ? exp(eta) : pow(eta, 1.0 / lp); } /** * derivative d(mu) / d(eta) * * @param eta linear predictors * @param lp link power * * @return d(mu)/d(eta) */ static R_INLINE double mu_eta(double eta, double lp){ return (lp == 0) ? exp(eta) : pow(eta, 1.0 / lp - 1.0) / lp ; } /** * Return the sum of squares of the first n elements of x * * @param n * @param x * * @return sum of squares */ static R_INLINE double sqr_length(double *x, int n) { double ans = 0; for (int i = 0; i < n; i++) ans += x[i] * x[i]; return ans; } /** * Return the index of the term associated with parameter index ind * * @param ind an index in [0, Gp[nt] - 1] * @param nt total number of terms * @param Gp group pointers, a vector of length nt+1 with Gp[0] = 0 * * @return sum of squares */ static R_INLINE int Gp_grp(int ind, int nt, const int *Gp) { for (int i = 0; i < nt; i++) if (ind < Gp[i + 1]) return i; error(("invalid row index %d (max is %d)"), ind, Gp[nt]); return -1; /* -Wall */ } /* prototypes defined in utilities.c */ double var(double *x, int n) ; void cov(int n, int p, double *x, double *ans); // matrix computation void solve_po(int d, double *v, double *iv) ; void mult_xtx(int m, int n, double *x, double *out) ; void mult_mv(char *trans, int m, int n, double *A, double *x, double *out) ; void chol(int d, double *v, double *iv); void solve_po(int d, double *v, double *iv); // numerical derivatives void grad(int n, double *x, double (*myfunc)(double *x, void *data), void *data, double *ans) ; void hess(int n, double *x, double (*myfunc)(double *x, void *data), void *data, double *ans) ; // wishart simulation void rwishart(int d, double nu, double *scal, double *out); #endif