From: Benjamin Auder Date: Thu, 12 Jan 2017 19:39:50 +0000 (+0100) Subject: C-part of tests ready [TODO: R part] X-Git-Url: https://git.auder.net/images/%7B%7B%20asset%28%27mixstore/css/store/upsert.css%27%29%20%7D%7D?a=commitdiff_plain;h=afa07d41c7592ac0ccd55d7af23c3bfef213291e;p=valse.git C-part of tests ready [TODO: R part] --- diff --git a/.gitignore b/.gitignore index b5f5640..767ee1a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ +#global ignores .Rhistory .RData -*.o -*.so +*.swp +*~ diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..d8d0ef0 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,5 @@ +#ignore object files, library and test executables +*.o +*.so +test.* +!test.*.c diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 42419ac..1439416 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -5,23 +5,23 @@ // TODO: don't recompute indexes every time...... void EMGLLF_core( // IN parameters - const double* phiInit, // parametre initial de moyenne renormalisé - const double* rhoInit, // parametre initial de variance renormalisé - const double* piInit, // parametre initial des proportions - const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon + const float* phiInit, // parametre initial de moyenne renormalisé + const float* rhoInit, // parametre initial de variance renormalisé + const float* piInit, // parametre initial des proportions + const float* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon int mini, // nombre minimal d'itérations dans l'algorithme EM int maxi, // nombre maximal d'itérations dans l'algorithme EM - double gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif - double lambda, // valeur du paramètre de régularisation du Lasso - const double* X, // régresseurs - const double* Y, // réponse - double tau, // seuil pour accepter la convergence + float gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif + float lambda, // valeur du paramètre de régularisation du Lasso + const float* X, // régresseurs + const float* Y, // réponse + float tau, // seuil pour accepter la convergence // OUT parameters (all pointers, to be modified) - double* phi, // parametre de moyenne renormalisé, calculé par l'EM - double* rho, // parametre de variance renormalisé, calculé par l'EM - double* pi, // parametre des proportions renormalisé, calculé par l'EM - double* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres - double* S, + float* phi, // parametre de moyenne renormalisé, calculé par l'EM + float* rho, // parametre de variance renormalisé, calculé par l'EM + float* pi, // parametre des proportions renormalisé, calculé par l'EM + float* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres + float* S, // additional size parameters int n, // nombre d'echantillons int p, // nombre de covariables @@ -37,32 +37,32 @@ void EMGLLF_core( //Other local variables //NOTE: variables order is always [maxi],n,p,m,k - double* gam = (double*)malloc(n*k*sizeof(double)); + float* gam = (float*)malloc(n*k*sizeof(float)); copyArray(gamInit, gam, n*k); - double* b = (double*)malloc(k*sizeof(double)); - double* Phi = (double*)malloc(p*m*k*sizeof(double)); - double* Rho = (double*)malloc(m*m*k*sizeof(double)); - double* Pi = (double*)malloc(k*sizeof(double)); - double* gam2 = (double*)malloc(k*sizeof(double)); - double* pi2 = (double*)malloc(k*sizeof(double)); - double* Gram2 = (double*)malloc(p*p*k*sizeof(double)); - double* ps = (double*)malloc(m*k*sizeof(double)); - double* nY2 = (double*)malloc(m*k*sizeof(double)); - double* ps1 = (double*)malloc(n*m*k*sizeof(double)); - double* ps2 = (double*)malloc(p*m*k*sizeof(double)); - double* nY21 = (double*)malloc(n*m*k*sizeof(double)); - double* Gam = (double*)malloc(n*k*sizeof(double)); - double* X2 = (double*)malloc(n*p*k*sizeof(double)); - double* Y2 = (double*)malloc(n*m*k*sizeof(double)); + float* b = (float*)malloc(k*sizeof(float)); + float* Phi = (float*)malloc(p*m*k*sizeof(float)); + float* Rho = (float*)malloc(m*m*k*sizeof(float)); + float* Pi = (float*)malloc(k*sizeof(float)); + float* gam2 = (float*)malloc(k*sizeof(float)); + float* pi2 = (float*)malloc(k*sizeof(float)); + float* Gram2 = (float*)malloc(p*p*k*sizeof(float)); + float* ps = (float*)malloc(m*k*sizeof(float)); + float* nY2 = (float*)malloc(m*k*sizeof(float)); + float* ps1 = (float*)malloc(n*m*k*sizeof(float)); + float* ps2 = (float*)malloc(p*m*k*sizeof(float)); + float* nY21 = (float*)malloc(n*m*k*sizeof(float)); + float* Gam = (float*)malloc(n*k*sizeof(float)); + float* X2 = (float*)malloc(n*p*k*sizeof(float)); + float* Y2 = (float*)malloc(n*m*k*sizeof(float)); gsl_matrix* matrix = gsl_matrix_alloc(m, m); gsl_permutation* permutation = gsl_permutation_alloc(m); - double* YiRhoR = (double*)malloc(m*sizeof(double)); - double* XiPhiR = (double*)malloc(m*sizeof(double)); - double dist = 0.; - double dist2 = 0.; + float* YiRhoR = (float*)malloc(m*sizeof(float)); + float* XiPhiR = (float*)malloc(m*sizeof(float)); + float dist = 0.; + float dist2 = 0.; int ite = 0; - double EPS = 1e-15; - double* dotProducts = (double*)malloc(k*sizeof(double)); + float EPS = 1e-15; + float* dotProducts = (float*)malloc(k*sizeof(float)); while (ite < mini || (ite < maxi && (dist >= tau || dist2 >= sqrt(tau)))) { @@ -90,7 +90,7 @@ void EMGLLF_core( //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r); for (int u=0; u - double dotProduct = 0.0; + float dotProduct = 0.0; for (int u=0; udata[u*m+v] = rho[ai(u,v,r,m,m,k)]; } gsl_linalg_LU_decomp(matrix, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrix, signum); + float detRhoR = gsl_linalg_LU_det(matrix, signum); Gam[mi(i,r,n,k)] = pi[r] * detRhoR * exp(-0.5*dotProducts[r] + shift); sumLLF1 += Gam[mi(i,r,n,k)] / pow(2*M_PI,m/2.0); @@ -327,7 +327,7 @@ void EMGLLF_core( } //sum(pen(ite,:)) - double sumPen = 0.0; + float sumPen = 0.0; for (int r=0; r Dist1) Dist1 = tmpDist; @@ -353,14 +353,14 @@ void EMGLLF_core( } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - double Dist2 = 0.0; + float Dist2 = 0.0; for (int u=0; u Dist2) Dist2 = tmpDist; @@ -368,12 +368,12 @@ void EMGLLF_core( } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - double Dist3 = 0.0; + float Dist3 = 0.0; for (int u=0; u Dist3) Dist3 = tmpDist; } diff --git a/src/sources/EMGLLF.h b/src/sources/EMGLLF.h index d09d27c..005c05b 100644 --- a/src/sources/EMGLLF.h +++ b/src/sources/EMGLLF.h @@ -3,23 +3,23 @@ void EMGLLF_core( // IN parameters - const double* phiInit, - const double* rhoInit, - const double* piInit, - const double* gamInit, + const float* phiInit, + const float* rhoInit, + const float* piInit, + const float* gamInit, int mini, int maxi, - double gamma, - double lambda, - const double* X, - const double* Y, - double tau, + float gamma, + float lambda, + const float* X, + const float* Y, + float tau, // OUT parameters - double* phi, - double* rho, - double* pi, - double* LLF, - double* S, + float* phi, + float* rho, + float* pi, + float* LLF, + float* S, // additional size parameters int n, int p, diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c index 0791892..3d5a9c7 100644 --- a/src/sources/EMGrank.c +++ b/src/sources/EMGrank.c @@ -3,13 +3,13 @@ #include "utils.h" // Compute pseudo-inverse of a square matrix -static double* pinv(const double* matrix, int dim) +static float* pinv(const float* matrix, int dim) { gsl_matrix* U = gsl_matrix_alloc(dim,dim); gsl_matrix* V = gsl_matrix_alloc(dim,dim); gsl_vector* S = gsl_vector_alloc(dim); gsl_vector* work = gsl_vector_alloc(dim); - double EPS = 1e-10; //threshold for singular value "== 0" + float EPS = 1e-10; //threshold for singular value "== 0" //copy matrix into U copyArray(matrix, U->data, dim*dim); @@ -19,12 +19,12 @@ static double* pinv(const double* matrix, int dim) gsl_vector_free(work); // Obtain pseudo-inverse by V*S^{-1}*t(U) - double* inverse = (double*)malloc(dim*dim*sizeof(double)); + float* inverse = (float*)malloc(dim*dim*sizeof(float)); for (int i=0; idata[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j]; inverse[i*dim+ii] = dotProduct; @@ -40,17 +40,17 @@ static double* pinv(const double* matrix, int dim) // TODO: comment EMGrank purpose void EMGrank_core( // IN parameters - const double* Pi, // parametre de proportion - const double* Rho, // parametre initial de variance renormalisé + const float* Pi, // parametre de proportion + const float* Rho, // parametre initial de variance renormalisé int mini, // nombre minimal d'itérations dans l'algorithme EM int maxi, // nombre maximal d'itérations dans l'algorithme EM - const double* X, // régresseurs - const double* Y, // réponse - double tau, // seuil pour accepter la convergence + const float* X, // régresseurs + const float* Y, // réponse + float tau, // seuil pour accepter la convergence const int* rank, // vecteur des rangs possibles // OUT parameters - double* phi, // parametre de moyenne renormalisé, calculé par l'EM - double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres + float* phi, // parametre de moyenne renormalisé, calculé par l'EM + float* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres // additional size parameters int n, // taille de l'echantillon int p, // nombre de covariables @@ -58,20 +58,20 @@ void EMGrank_core( int k) // nombre de composantes { // Allocations, initializations - double* Phi = (double*)calloc(p*m*k,sizeof(double)); - double* hatBetaR = (double*)malloc(p*m*sizeof(double)); + float* Phi = (float*)calloc(p*m*k,sizeof(float)); + float* hatBetaR = (float*)malloc(p*m*sizeof(float)); int signum; - double invN = 1.0/n; + float invN = 1.0/n; int deltaPhiBufferSize = 20; - double* deltaPhi = (double*)malloc(deltaPhiBufferSize*sizeof(double)); + float* deltaPhi = (float*)malloc(deltaPhiBufferSize*sizeof(float)); int ite = 0; - double sumDeltaPhi = 0.0; - double* YiRhoR = (double*)malloc(m*sizeof(double)); - double* XiPhiR = (double*)malloc(m*sizeof(double)); - double* Xr = (double*)malloc(n*p*sizeof(double)); - double* Yr = (double*)malloc(n*m*sizeof(double)); - double* tXrXr = (double*)malloc(p*p*sizeof(double)); - double* tXrYr = (double*)malloc(p*m*sizeof(double)); + float sumDeltaPhi = 0.0; + float* YiRhoR = (float*)malloc(m*sizeof(float)); + float* XiPhiR = (float*)malloc(m*sizeof(float)); + float* Xr = (float*)malloc(n*p*sizeof(float)); + float* Yr = (float*)malloc(n*m*sizeof(float)); + float* tXrXr = (float*)malloc(p*p*sizeof(float)); + float* tXrYr = (float*)malloc(p*m*sizeof(float)); gsl_matrix* matrixM = gsl_matrix_alloc(p, m); gsl_matrix* matrixE = gsl_matrix_alloc(m, m); gsl_permutation* permutation = gsl_permutation_alloc(m); @@ -115,7 +115,7 @@ void EMGrank_core( { for (int jj=0; jjdata[j*m+jj] = dotProduct; @@ -164,7 +164,7 @@ void EMGrank_core( { for (int jj=0; jjdata[u] * V->data[jj*m+u]; hatBetaR[mi(j,jj,p,m)] = dotProduct; @@ -176,7 +176,7 @@ void EMGrank_core( { for (int jj=0; jjdata[j*m+jj] = Rho[ai(j,jj,r,m,m,k)]; } gsl_linalg_LU_decomp(matrixE, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrixE, signum); + float detRhoR = gsl_linalg_LU_det(matrixE, signum); //compute Y(i,:)*Rho(:,:,r) for (int j=0; j - double dotProduct = 0.0; + float dotProduct = 0.0; for (int u=0; u maxLogGamIR) @@ -247,14 +247,14 @@ void EMGrank_core( *LLF = -invN * sumLogLLF2; //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi)))); - double newDeltaPhi = 0.0; + float newDeltaPhi = 0.0; for (int j=0; j newDeltaPhi) newDeltaPhi = tmpDist; diff --git a/src/sources/EMGrank.h b/src/sources/EMGrank.h index 8123e68..92c3f73 100644 --- a/src/sources/EMGrank.h +++ b/src/sources/EMGrank.h @@ -3,17 +3,17 @@ void EMGrank_core( // IN parameters - const double* Pi, - const double* Rho, + const float* Pi, + const float* Rho, int mini, int maxi, - const double* X, - const double* Y, - double tau, + const float* X, + const float* Y, + float tau, const int* rank, // OUT parameters - double* phi, - double* LLF, + float* phi, + float* LLF, // additional size parameters int n, int p, diff --git a/src/sources/constructionModelesLassoMLE.c b/src/sources/constructionModelesLassoMLE.c index 4ab62ad..a48f7d0 100644 --- a/src/sources/constructionModelesLassoMLE.c +++ b/src/sources/constructionModelesLassoMLE.c @@ -7,25 +7,25 @@ // TODO: comment on constructionModelesLassoMLE purpose void constructionModelesLassoMLE_core( // IN parameters - const double* phiInit, // parametre initial de moyenne renormalisé - const double* rhoInit, // parametre initial de variance renormalisé - const double* piInit,// parametre initial des proportions - const double* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon + const float* phiInit, // parametre initial de moyenne renormalisé + const float* rhoInit, // parametre initial de variance renormalisé + const float* piInit,// parametre initial des proportions + const float* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon int mini,// nombre minimal d'itérations dans l'algorithme EM int maxi,// nombre maximal d'itérations dans l'algorithme EM - double gamma,// valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif - const double* glambda, // valeur des paramètres de régularisation du Lasso - const double* X, // régresseurs - const double* Y, // réponse - double seuil,// seuil pour prendre en compte une variable - double tau,// seuil pour accepter la convergence + float gamma,// valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif + const float* glambda, // valeur des paramètres de régularisation du Lasso + const float* X, // régresseurs + const float* Y, // réponse + float seuil,// seuil pour prendre en compte une variable + float tau,// seuil pour accepter la convergence const int* A1, // matrice des coefficients des parametres selectionnes const int* A2, // matrice des coefficients des parametres non selectionnes // OUT parameters - double* phi,// estimateur ainsi calculé par le Lasso - double* rho,// estimateur ainsi calculé par le Lasso - double* pi, // estimateur ainsi calculé par le Lasso - double* lvraisemblance, // estimateur ainsi calculé par le Lasso + float* phi,// estimateur ainsi calculé par le Lasso + float* rho,// estimateur ainsi calculé par le Lasso + float* pi, // estimateur ainsi calculé par le Lasso + float* lvraisemblance, // estimateur ainsi calculé par le Lasso // additional size parameters int n, // taille de l'echantillon int p, // nombre de covariables @@ -58,7 +58,7 @@ void constructionModelesLassoMLE_core( continue; //Xa = X(:,a) - double* Xa = (double*)malloc(n*lengthA*sizeof(double)); + float* Xa = (float*)malloc(n*lengthA*sizeof(float)); for (int i=0; idata[u*m+v] = rho[ai4(u,v,r,lambdaIndex,m,m,k,L)]; } gsl_linalg_LU_decomp(matrix, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrix, signum); + float detRhoR = gsl_linalg_LU_det(matrix, signum); //compute Y(i,:)*rho(:,:,r,lambdaIndex) for (int u=0; u - double dotProduct = 0.0; + float dotProduct = 0.0; for (int u=0; u 0 - double* phiLambda = (double*)malloc(longueurActive*m*k*sizeof(double)); - double LLF; + float* phiLambda = (float*)malloc(longueurActive*m*k*sizeof(float)); + float LLF; for (int j=0; j maxPhi) diff --git a/src/sources/selectiontotale.h b/src/sources/selectiontotale.h index fa2ba84..d225592 100644 --- a/src/sources/selectiontotale.h +++ b/src/sources/selectiontotale.h @@ -4,23 +4,23 @@ // Main job on raw inputs (after transformation from mxArray) void selectiontotale_core( // IN parameters - const double* phiInit, - const double* rhoInit, - const double* piInit, - const double* gamInit, + const float* phiInit, + const float* rhoInit, + const float* piInit, + const float* gamInit, int mini, int maxi, - double gamma, - const double* glambda, - const double* X, - const double* Y, - double seuil, - double tau, + float gamma, + const float* glambda, + const float* X, + const float* Y, + float seuil, + float tau, // OUT parameters int* A1, int* A2, - double* Rho, - double* Pi, + float* Rho, + float* Pi, // additional size parameters int n, int p, diff --git a/src/test/.Makefile.swp b/src/test/.Makefile.swp deleted file mode 100644 index bfa54d3..0000000 Binary files a/src/test/.Makefile.swp and /dev/null differ diff --git a/src/test/Makefile b/src/test/Makefile index 796b424..459a5cb 100644 --- a/src/test/Makefile +++ b/src/test/Makefile @@ -1,8 +1,8 @@ CC = gcc CFLAGS = -g -std=gnu99 -Wno-implicit-function-declaration LDFLAGS = -lm -lgsl -lcblas -lgomp -LDFLAGS_TEST = -Lobj/ -lvalse_core -LIB = valse_core.so +TEST_LDFLAGS = -L. libvalse_core.so +LIB = libvalse_core.so LIB_SRC = $(wildcard ../sources/*.c) LIB_OBJ = $(LIB_SRC:.c=.o) INCLUDES = -I../sources @@ -12,23 +12,23 @@ all: $(LIB) test.EMGLLF test.EMGrank test.constructionModelesLassoMLE test.EMGra $(LIB): $(LIB_OBJ) $(CC) -shared -o $@ $^ $(LDFLAGS) -test.EMGLLF: test.EMGLLF.o +test.EMGLLF: test.EMGLLF.o utils.o $(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS) -test.constructionModelesLassoMLE: test.constructionModelesLassoMLE.o +test.constructionModelesLassoMLE: test.constructionModelesLassoMLE.o utils.o $(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS) -test.EMGrank: test.EMGrank.o +test.EMGrank: test.EMGrank.o utils.o $(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS) -test.constructionModelesLassoRank: test.constructionModelesLassoRank.o +test.constructionModelesLassoRank: test.constructionModelesLassoRank.o utils.o $(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS) -test.selectionTotale: test.selectionTotale.o +test.selectionTotale: test.selectionTotale.o utils.o $(CC) -o $@ $^ $(LDFLAGS) $(TEST_LDFLAGS) %.o: %.c - $(CC) -o $@ -c $< $(CFLAGS) $(INCLUDES) + $(CC) -fPIC -o $@ -c $< $(CFLAGS) $(INCLUDES) clean: rm -f *.o ../sources/*.o diff --git a/src/test/OLD_TEST_MATLAB/TODO b/src/test/OLD_TEST_MATLAB/TODO index d4b05cf..73301dd 100644 --- a/src/test/OLD_TEST_MATLAB/TODO +++ b/src/test/OLD_TEST_MATLAB/TODO @@ -1,4 +1,3 @@ 1) transformer les .m en .R -2) sauvegarder les résultats + entrée/sorties à partir du code R uniquement -3) dans test.Truc.c, prendre en entrée les paramètres de data/ juste sauvegardés, - puis comparer les sorties aux sorties enregistrées +2) sauvegarder les résultats + entrée/sorties à partir du code R uniquement, + dans le dossier src/test/data (** sous forme de vecteurs, sep=" " **) diff --git a/src/test/test.EMGLLF.c b/src/test/test.EMGLLF.c index e938b3a..1720641 100644 --- a/src/test/test.EMGLLF.c +++ b/src/test/test.EMGLLF.c @@ -68,8 +68,8 @@ int main(int argc, char** argv) free(pi); free(ref_pi); - float* ref_LLF = readArray_real("LLF", maxi); - compareArray_real("LLF", LLF, ref_LLF); + float* ref_LLF = readArray_real("LLF"); + compareArray_real("LLF", LLF, ref_LLF, maxi); free(LLF); free(ref_LLF); diff --git a/src/test/utils.c b/src/test/utils.c index 8895cda..a6370b0 100644 --- a/src/test/utils.c +++ b/src/test/utils.c @@ -1,3 +1,8 @@ +#include +#include +#include +#include + // Check if array == refArray void compareArray(const char* ID, const void* array, const void* refArray, int size, int isinteger) @@ -44,7 +49,7 @@ void* readArray(const char* fileName, int isinteger) strcat(command, " | wc -l"); FILE *countSpaces = popen(command, "r"); char* buffer = (char*)calloc(32, sizeof(char)); - fgets(buffer, sizeof(buffer), command); + fgets(buffer, sizeof(buffer), countSpaces); int n = atoi(buffer) + 1; free(buffer); pclose(countSpaces); @@ -92,7 +97,7 @@ int* readArray_int(const char* fileName) float* readArray_real(const char* fileName) { - return (int*)readArray(fileName, 0); + return (float*)readArray(fileName, 0); } int read_int(const char* fileName)