Skip to content

Commit

Permalink
update C to avoid "Strict" errors/issues
Browse files Browse the repository at this point in the history 
R_Alloc instead of Alloc etc.
  • Loading branch information
fabian-s committed Sep 3, 2024
1 parent ac67c66 commit 2474d6f
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 86 deletions.
94 changes: 47 additions & 47 deletions src/sampler.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,40 +93,40 @@ SEXP sampler(
save = 0, keep = *burnin, nrv =1, info=0,
nSamp=(*totalLength-*burnin)/(*thin), oneInt = 1, zeroInt = 0;

double *p1 =Calloc(*pPen, double);
double *p1 =R_Calloc(*pPen, double);

double infV = 100000, oneV = 1.0, zeroV = 0.0, minusOneV =-1.0;
double *one=&oneV, *zero=&zeroV, *minusOne=&minusOneV, *inf=&infV, acceptance=0;
double invSigma2 = 1 / *sigma2, sqrtInvSigma2 = R_pow(invSigma2, 0.5);
double *penAlphaSq, *alphaLong, *varAlpha, *priorMeanAlpha, *modeAlpha, *offsetAlpha;;
penAlphaSq = Calloc(*pPen, double);
penAlphaSq = R_Calloc(*pPen, double);
for(int i=*p-*pPen; i<*p; i++) penAlphaSq[i- *p + *pPen] = R_pow(alpha[i], 2.0);
alphaLong = Calloc(*q, double);
alphaLong = R_Calloc(*q, double);
F77_CALL(dgemm)("N","N", q, &oneInt, p, one, G, q, alpha, p, zero, alphaLong, q FCONE FCONE);
varAlpha = Calloc(*p, double);
varAlpha = R_Calloc(*p, double);
for(int i=0; i<startPen; i++) varAlpha[i] = *inf; /*unpenalized*/
for(int i=startPen; i<*p; i++) varAlpha[i] = tau2[i-startPen]*gamma[i-startPen]; /*penalized*/
priorMeanAlpha = Calloc(*p, double);
priorMeanAlpha = R_Calloc(*p, double);
setToZero(priorMeanAlpha, *p);
modeAlpha = Calloc(*p, double);
modeAlpha = R_Calloc(*p, double);
F77_CALL(dcopy)(p, alpha, &oneInt, modeAlpha, &oneInt);
offsetAlpha = Calloc(*n, double);
offsetAlpha = R_Calloc(*n, double);
F77_CALL(dcopy)(n, offset, &oneInt, offsetAlpha, &oneInt);


double *ksiUpdate, *priorMeanKsi, *modeKsi, *offsetKsi;
int safeQKsiUpdate = imax2(1, *qKsiUpdate);
//ksiUpdate contains the last qKsiUpdate elements in ksi
ksiUpdate = Calloc(safeQKsiUpdate, double);
ksiUpdate = R_Calloc(safeQKsiUpdate, double);
F77_CALL(dcopy)(&safeQKsiUpdate, &ksi[*q-safeQKsiUpdate], &oneInt, ksiUpdate, &oneInt);
priorMeanKsi = Calloc(safeQKsiUpdate, double);
priorMeanKsi = R_Calloc(safeQKsiUpdate, double);
setToZero(priorMeanKsi, safeQKsiUpdate);
for(int i=0; i<*qKsiUpdate; i++) priorMeanKsi[i] = 1.0;
modeKsi = Calloc(safeQKsiUpdate, double);
modeKsi = R_Calloc(safeQKsiUpdate, double);
setToZero(modeKsi, safeQKsiUpdate);
for(int i=0; i<*qKsiUpdate; i++) modeKsi[i] = ksi[i+qKsiNoUpdate];
// offsetKsi = offset + X_d=1*alpha : use lin.predictor of grps with ksi==1 as offset
offsetKsi = Calloc(*n, double);
offsetKsi = R_Calloc(*n, double);
F77_CALL(dcopy)(n, offset, &oneInt, offsetKsi, &oneInt);
if(qKsiNoUpdate < *q){
if(qKsiNoUpdate > 0){
Expand All @@ -135,31 +135,31 @@ SEXP sampler(
}

double *eta, *resid, rss, *XAlpha, *XKsiUpdate, *etaOffset;
eta = Calloc(*n, double);
eta = R_Calloc(*n, double);
F77_CALL(dgemm)("N","N", n, &oneInt, q, one, X, n, beta, q, zero, eta, n FCONE FCONE);
resid = Calloc(*n, double);
resid = R_Calloc(*n, double);
rss = 0;
for(int i=0; i<*n; i++) {
resid[i] = y[i]-eta[i] - offset[i];
rss += R_pow(resid[i], 2.0);
}
XAlpha = Calloc(*p * (*n), double);
XAlpha = R_Calloc(*p * (*n), double);
updateXAlpha(XAlpha, X, G, ksi, q, qKsiUpdate, p, n);
XKsiUpdate = Calloc( *n * safeQKsiUpdate, double);
XKsiUpdate = R_Calloc( *n * safeQKsiUpdate, double);
setToZero(XKsiUpdate, *n * safeQKsiUpdate);
if(qKsiNoUpdate < *q){
updateXKsi(XKsiUpdate, X, alphaLong, q, &qKsiNoUpdate, n);
}
etaOffset = Calloc(*n, double);
etaOffset = R_Calloc(*n, double);
for(int i=0; i<*n; i++) etaOffset[i] = eta[i]+offset[i];


// ############################################################ //
// ######## set up blocks for blockwise updates ############### //
// ############################################################ //

XBlockQR *AlphaBlocks = Calloc(*blocksAlpha, XBlockQR);
XBlockQR *KsiBlocks = Calloc(*blocksKsi, XBlockQR);
XBlockQR *AlphaBlocks = R_Calloc(*blocksAlpha, XBlockQR);
XBlockQR *KsiBlocks = R_Calloc(*blocksKsi, XBlockQR);


for(int i=0; i < *blocksAlpha; i++){
Expand All @@ -169,27 +169,27 @@ SEXP sampler(
(AlphaBlocks[i]).qA = (AlphaBlocks[i]).indA2 - (AlphaBlocks[i]).indA1 + 1;
(AlphaBlocks[i]).qI = *p - (AlphaBlocks[i]).qA;

(AlphaBlocks[i]).qraux = Calloc((AlphaBlocks[i]).qA, double);
(AlphaBlocks[i]).qraux = R_Calloc((AlphaBlocks[i]).qA, double);
setToZero((AlphaBlocks[i]).qraux, (AlphaBlocks[i]).qA);
(AlphaBlocks[i]).work = Calloc((AlphaBlocks[i]).qA, double);
(AlphaBlocks[i]).work = R_Calloc((AlphaBlocks[i]).qA, double);
setToZero((AlphaBlocks[i]).work, (AlphaBlocks[i]).qA);
(AlphaBlocks[i]).pivots = Calloc((AlphaBlocks[i]).qA, int);
(AlphaBlocks[i]).pivots = R_Calloc((AlphaBlocks[i]).qA, int);
for(int j=0; j < (AlphaBlocks[i]).qA; j++) (AlphaBlocks[i]).pivots[j] = 0;

(AlphaBlocks[i]).coefI = Calloc((AlphaBlocks[i]).qI, double);
(AlphaBlocks[i]).coefI = R_Calloc((AlphaBlocks[i]).qI, double);
setToZero((AlphaBlocks[i]).coefI, (AlphaBlocks[i]).qI);

(AlphaBlocks[i]).Xa = Calloc(((AlphaBlocks[i]).qA + *n) * (AlphaBlocks[i]).qA, double);
(AlphaBlocks[i]).Xa = R_Calloc(((AlphaBlocks[i]).qA + *n) * (AlphaBlocks[i]).qA, double);
setToZero((AlphaBlocks[i]).Xa, ((AlphaBlocks[i]).qA + *n) * (AlphaBlocks[i]).qA);
(AlphaBlocks[i]).Xi = Calloc(*n * (AlphaBlocks[i]).qI, double);
(AlphaBlocks[i]).Xi = R_Calloc(*n * (AlphaBlocks[i]).qI, double);
setToZero((AlphaBlocks[i]).Xi, *n * (AlphaBlocks[i]).qI );
(AlphaBlocks[i]).ya = Calloc(((AlphaBlocks[i]).qA + *n), double);
(AlphaBlocks[i]).ya = R_Calloc(((AlphaBlocks[i]).qA + *n), double);
F77_CALL(dcopy)(n, y, &nrv, (AlphaBlocks[i]).ya, &nrv);
setToZero((AlphaBlocks[i]).ya + *n, (AlphaBlocks[i]).qA);

(AlphaBlocks[i]).m = Calloc((AlphaBlocks[i]).qA, double);
(AlphaBlocks[i]).m = R_Calloc((AlphaBlocks[i]).qA, double);
setToZero((AlphaBlocks[i]).m, (AlphaBlocks[i]).qA);
(AlphaBlocks[i]).err = Calloc((AlphaBlocks[i]).qA, double);
(AlphaBlocks[i]).err = R_Calloc((AlphaBlocks[i]).qA, double);
setToZero((AlphaBlocks[i]).err, (AlphaBlocks[i]).qA);

}
Expand All @@ -204,27 +204,27 @@ SEXP sampler(
(KsiBlocks[i]).qA = (KsiBlocks[i]).indA2 - (KsiBlocks[i]).indA1 + 1;
(KsiBlocks[i]).qI = *qKsiUpdate - (KsiBlocks[i]).qA;

(KsiBlocks[i]).qraux = Calloc((KsiBlocks[i]).qA, double);
(KsiBlocks[i]).qraux = R_Calloc((KsiBlocks[i]).qA, double);
setToZero((KsiBlocks[i]).qraux, (KsiBlocks[i]).qA);
(KsiBlocks[i]).work = Calloc((KsiBlocks[i]).qA, double);
(KsiBlocks[i]).work = R_Calloc((KsiBlocks[i]).qA, double);
setToZero((KsiBlocks[i]).work, (KsiBlocks[i]).qA);
(KsiBlocks[i]).pivots = Calloc((KsiBlocks[i]).qA, int);
(KsiBlocks[i]).pivots = R_Calloc((KsiBlocks[i]).qA, int);
for(int j=0; j < (KsiBlocks[i]).qA; j++) (KsiBlocks[i]).pivots[j] = 0;

(KsiBlocks[i]).coefI = Calloc((KsiBlocks[i]).qI, double);
(KsiBlocks[i]).coefI = R_Calloc((KsiBlocks[i]).qI, double);
setToZero((KsiBlocks[i]).coefI, (KsiBlocks[i]).qI);

(KsiBlocks[i]).Xa = Calloc(((KsiBlocks[i]).qA + *n) * (KsiBlocks[i]).qA, double);
(KsiBlocks[i]).Xa = R_Calloc(((KsiBlocks[i]).qA + *n) * (KsiBlocks[i]).qA, double);
setToZero((KsiBlocks[i]).Xa, ((KsiBlocks[i]).qA + *n) * (KsiBlocks[i]).qA);
(KsiBlocks[i]).Xi = Calloc(*n * (KsiBlocks[i]).qI, double);
(KsiBlocks[i]).Xi = R_Calloc(*n * (KsiBlocks[i]).qI, double);
setToZero((KsiBlocks[i]).Xi, *n * (KsiBlocks[i]).qI );
(KsiBlocks[i]).ya = Calloc(((KsiBlocks[i]).qA + *n), double);
(KsiBlocks[i]).ya = R_Calloc(((KsiBlocks[i]).qA + *n), double);
F77_CALL(dcopy)(n, y, &nrv, (KsiBlocks[i]).ya, &nrv);
setToZero((KsiBlocks[i]).ya + *n, (KsiBlocks[i]).qA);

(KsiBlocks[i]).m = Calloc((KsiBlocks[i]).qA, double);
(KsiBlocks[i]).m = R_Calloc((KsiBlocks[i]).qA, double);
setToZero((KsiBlocks[i]).m, (KsiBlocks[i]).qA);
(KsiBlocks[i]).err = Calloc((KsiBlocks[i]).qA, double);
(KsiBlocks[i]).err = R_Calloc((KsiBlocks[i]).qA, double);
setToZero((KsiBlocks[i]).err, (KsiBlocks[i]).qA);
}
initializeBlocksQR(KsiBlocks, XKsiUpdate, *n, *blocksKsi, *qKsiUpdate, varKsi, scale);
Expand Down Expand Up @@ -405,17 +405,17 @@ SEXP sampler(
}
}

Free(etaOffset); Free(XKsiUpdate); Free(XAlpha); Free(resid); Free(eta);
Free(offsetKsi); Free(modeKsi); Free(priorMeanKsi); Free(ksiUpdate);
Free(offsetAlpha);
Free(modeAlpha);
Free(priorMeanAlpha);
Free(varAlpha);
Free(alphaLong);
Free(penAlphaSq);
freeXBlockQR(AlphaBlocks, *blocksAlpha);
if(qKsiNoUpdate < *q) freeXBlockQR(KsiBlocks, *blocksKsi);
Free(p1);
R_Free(etaOffset); R_Free(XKsiUpdate); R_Free(XAlpha); R_Free(resid); R_Free(eta);
R_Free(offsetKsi); R_Free(modeKsi); R_Free(priorMeanKsi); R_Free(ksiUpdate);
R_Free(offsetAlpha);
R_Free(modeAlpha);
R_Free(priorMeanAlpha);
R_Free(varAlpha);
R_Free(alphaLong);
R_Free(penAlphaSq);
R_FreeXBlockQR(AlphaBlocks, *blocksAlpha);
if(qKsiNoUpdate < *q) R_FreeXBlockQR(KsiBlocks, *blocksKsi);
R_Free(p1);
return(R_NilValue);
}/* end sampler ()*/

52 changes: 26 additions & 26 deletions src/updaters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ static void updateXAlpha(double *XAlpha, double *X, double *G, double *ksi, int
//calculates X*diag(ksi)*G
int i,j, nqNoUpdate=(*q - *qKsiUpdate) * (*n), nrv=1;
double one=1.0, zero=0.0;
double *temp = Calloc(*q * (*n), double);
double *temp = R_Calloc(*q * (*n), double);

//copy first qKsiNoUpdate columns of X to temp
if(nqNoUpdate > 0) F77_CALL(dcopy)(&nqNoUpdate, X, &nrv, temp, &nrv);
Expand All @@ -29,7 +29,7 @@ static void updateXAlpha(double *XAlpha, double *X, double *G, double *ksi, int
}
//XAlpha = temp %*%G = (X*ksi) %*% G
F77_CALL(dgemm)("N","N", n, p, q, &one, temp, n, G, q, &zero, XAlpha, n FCONE FCONE);
Free(temp);
R_Free(temp);
}

static void updateXKsi(double *XKsiUpdate, double *X, double *alphaLong, int *q, int *qKsiNoUpdate, int *n){
Expand Down Expand Up @@ -197,13 +197,13 @@ static void updateCoefBlockQR_IWLS(int j, int n, double *X, double *y,
int i, k, l, info=0, oneInt = 1, switchMode = 0;
int qA = (Xblocks[j]).qA, qI = (Xblocks[j]).qI, rows = n + qA;
double minusOne = -1.0, one=1.0, likC, logAccProb, proposalDens;
double *coefC = Calloc(qA, double);
double *coefChange = Calloc(qA, double);
double *var = Calloc(n, double);
double *w = Calloc(n, double);
double *mu = Calloc(n, double);
double *etaOffset = Calloc(n, double);
double *etaMode= Calloc(n, double);
double *coefC = R_Calloc(qA, double);
double *coefChange = R_Calloc(qA, double);
double *var = R_Calloc(n, double);
double *w = R_Calloc(n, double);
double *mu = R_Calloc(n, double);
double *etaOffset = R_Calloc(n, double);
double *etaMode= R_Calloc(n, double);

//0. update Xa, Xi (see updateBlocks)
l = 0;
Expand Down Expand Up @@ -416,13 +416,13 @@ static void updateCoefBlockQR_IWLS(int j, int n, double *X, double *y,
}
}

Free(etaMode);
Free(etaOffset);
Free(mu);
Free(w);
Free(var);
Free(coefChange);
Free(coefC);
R_Free(etaMode);
R_Free(etaOffset);
R_Free(mu);
R_Free(w);
R_Free(var);
R_Free(coefChange);
R_Free(coefC);
}

static void updateCoefQR(
Expand All @@ -443,9 +443,9 @@ static void updateCoefQR(
){
// FIXME: DIRTYDIRTY- scale[0] is 1/sqrt(sigma2) for gaussian, else the real scale parameter (e.g. n_i for binomial)
int j, k, length=blocksCoef;
double *yOffset = Calloc(n, double);
int *shuffle = Calloc(blocksCoef, int);
int *iterator = Calloc(blocksCoef, int);
double *yOffset = R_Calloc(n, double);
int *shuffle = R_Calloc(blocksCoef, int);
int *iterator = R_Calloc(blocksCoef, int);

if(family == 0) { // offset for IWLS incorporated into determining working obs, only need this for non-IWLS
for(j = 0; j < n; j++) yOffset[j] = y[j] - offset[j];
Expand All @@ -471,19 +471,19 @@ static void updateCoefQR(
}

}
Free(iterator);
Free(shuffle);
Free(yOffset);
R_Free(iterator);
R_Free(shuffle);
R_Free(yOffset);
}

static void rescaleKsiAlpha(double *ksi, double *alpha, double *varKsi, double *tau2, double *G, int *d, int p, int q,
int qKsiNoUpdate, int pPen, int scaleMode, double *modeAlpha, double *modeKsi, int family){
/*rescale ksi groupwise so that sum(abs(ksi_g))/dim(g)=1, change alphas accordingly*/
/* also tried rescaling by 1/max(abs(ksi_g)) --> stronger regularization*/
int i ,j, oneInt =1;
double *scale = Calloc(p, double);
double *scale = R_Calloc(p, double);
setToZero(scale, p);
double *scaleLong = Calloc(q, double);
double *scaleLong = R_Calloc(q, double);
setToZero(scaleLong, q);
double one = 1.0, zero=0.0;

Expand Down Expand Up @@ -523,8 +523,8 @@ static void rescaleKsiAlpha(double *ksi, double *alpha, double *varKsi, double
//varKsi[j-qKsiNoUpdate] *= R_pow(scaleLong[j], 2.0);
}
}
Free(scaleLong);
Free(scale);
R_Free(scaleLong);
R_Free(scale);
}

static double updateLogPost(
Expand Down
27 changes: 14 additions & 13 deletions src/utils.h
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,14 @@
#ifndef MYUTILS_H_
#define MYUTILS_H_

/* R_alloc*/

#include <R.h>
/* FlushConsole*/
#include <Rinternals.h>
/* rng's, log, etc..*/
#include <Rmath.h>
/* R_alloc*/
#include <R_ext/RS.h>
/* BLAS */
#include <R_ext/BLAS.h>
/* Lapack*/
Expand All @@ -15,7 +17,6 @@
#include <R_ext/Linpack.h>
/* CheckUserInterrupt*/
#include <R_ext/Utils.h>

/*Fortran routines DQRLS etc...*/
#include <R_ext/Applic.h>

Expand Down Expand Up @@ -44,22 +45,22 @@ typedef struct {
} XBlockQR;


void freeXBlockQR(XBlockQR *blocks, int n){
void R_FreeXBlockQR(XBlockQR *blocks, int n){
for (int i =0; i < n; i++) {
Free(blocks[i].ya);
Free(blocks[i].Xa);
Free(blocks[i].Xi);
Free(blocks[i].coefI);
R_Free(blocks[i].ya);
R_Free(blocks[i].Xa);
R_Free(blocks[i].Xi);
R_Free(blocks[i].coefI);

Free(blocks[i].qraux);
Free(blocks[i].work);
Free(blocks[i].pivots);
R_Free(blocks[i].qraux);
R_Free(blocks[i].work);
R_Free(blocks[i].pivots);

Free(blocks[i].m);
Free(blocks[i].err);
R_Free(blocks[i].m);
R_Free(blocks[i].err);

}
Free(blocks);
R_Free(blocks);
}


Expand Down

0 comments on commit 2474d6f

Please sign in to comment.