From 46a2e676b2d85eef5a1811a6e623b65327fc451d Mon Sep 17 00:00:00 2001 From: Benjamin Auder <benjamin.auder@somewhere> Date: Sun, 12 Feb 2017 23:31:28 +0100 Subject: [PATCH] fix test.constructionModelesLassoMLE; TODO: det(rho[,,r]) ? numerical errors ? ...and selectiontotale.m --- src/sources/EMGLLF.c | 4 +- src/sources/constructionModelesLassoMLE.c | 60 +++++++-------- src/sources/constructionModelesLassoRank.c | 2 + src/sources/utils.h | 4 - src/test/Makefile | 2 +- ...eRunSaveTest_constructionModelesLassoMLE.R | 2 +- src/test/generate_test_data/helpers/EMGLLF.R | 2 +- .../helpers/constructionModelesLassoMLE.R | 74 +++++++++---------- .../helpers/constructionModelesLassoMLE.m | 58 --------------- .../helpers/constructionModelesLassoRank.R | 2 +- ...E.c => test.constructionModelesLassoMLE.c} | 13 ++-- 11 files changed, 80 insertions(+), 143 deletions(-) delete mode 100644 src/test/generate_test_data/helpers/constructionModelesLassoMLE.m rename src/test/{test.ConstructionModelesLassoMLE.c => test.constructionModelesLassoMLE.c} (88%) diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 7e7a3d1..86b6060 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -92,7 +92,7 @@ void EMGLLF_core( { Real dotProduct = 0.; for (int v=0; v<n; v++) - dotProduct += X2[ai(v,u,r,n,m,k)] * Y2[ai(v,mm,r,n,m,k)]; + dotProduct += X2[ai(v,u,r,n,p,k)] * Y2[ai(v,mm,r,n,m,k)]; ps2[ai(u,mm,r,p,m,k)] = dotProduct; } } @@ -317,7 +317,7 @@ void EMGLLF_core( Real detRhoR = gsl_linalg_LU_det(matrix, signum); //FIXME: det(rho[,,r]) too small(?!). See EMGLLF.R - Gam[mi(i,r,n,k)] = pi[r] * exp(-0.5*sqNorm2[r] + shift) * detRhoR; + Gam[mi(i,r,n,k)] = pi[r] * exp(-0.5*sqNorm2[r] + shift) ; //* detRhoR; sumLLF1 += Gam[mi(i,r,n,k)] / gaussConstM; sumGamI += Gam[mi(i,r,n,k)]; } diff --git a/src/sources/constructionModelesLassoMLE.c b/src/sources/constructionModelesLassoMLE.c index ad2a718..34e5808 100644 --- a/src/sources/constructionModelesLassoMLE.c +++ b/src/sources/constructionModelesLassoMLE.c @@ -13,7 +13,8 @@ void constructionModelesLassoMLE_core( const Real* 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 - Real gamma,// valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif + Real gamma,// valeur de gamma : puissance des proportions dans la pénalisation + //pour un Lasso adaptatif const Real* glambda, // valeur des paramètres de régularisation du Lasso const Real* X, // régresseurs const Real* Y, // réponse @@ -33,9 +34,15 @@ void constructionModelesLassoMLE_core( int k, // nombre de composantes int L) // taille de glambda { - //preparation: phi = 0 + //preparation: phi,rho,pi = 0, llh=+Inf for (int u=0; u<p*m*k*L; u++) - phi[u] = 0.0; + phi[u] = 0.; + for (int u=0; u<m*m*k*L; u++) + rho[u] = 0.; + for (int u=0; u<k*L; u++) + pi[u] = 0.; + for (int u=0; u<L*2; u++) + llh[u] = INFINITY; //initiate parallel section int lambdaIndex; @@ -45,8 +52,7 @@ void constructionModelesLassoMLE_core( #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) { - //~ a = A1(:,1,lambdaIndex); - //~ a(a==0) = []; + //a = A1[,1,lambdaIndex] ; a = a[a!=0] int* a = (int*)malloc(p*sizeof(int)); int lengthA = 0; for (int j=0; j<p; j++) @@ -55,9 +61,12 @@ void constructionModelesLassoMLE_core( a[lengthA++] = A1[ai(j,0,lambdaIndex,p,m+1,L)] - 1; } if (lengthA == 0) + { + free(a); continue; + } - //Xa = X(:,a) + //Xa = X[,a] Real* Xa = (Real*)malloc(n*lengthA*sizeof(Real)); for (int i=0; i<n; i++) { @@ -65,7 +74,7 @@ void constructionModelesLassoMLE_core( Xa[mi(i,j,n,lengthA)] = X[mi(i,a[j],n,p)]; } - //phia = phiInit(a,:,:) + //phia = phiInit[a,,] Real* phia = (Real*)malloc(lengthA*m*k*sizeof(Real)); for (int j=0; j<lengthA; j++) { @@ -76,14 +85,13 @@ void constructionModelesLassoMLE_core( } } - //[phiLambda,rhoLambda,piLambda,~,~] = EMGLLF(... - // phiInit(a,:,:),rhoInit,piInit,gamInit,mini,maxi,gamma,0,X(:,a),Y,tau); + //Call to EMGLLF Real* phiLambda = (Real*)malloc(lengthA*m*k*sizeof(Real)); Real* rhoLambda = (Real*)malloc(m*m*k*sizeof(Real)); Real* piLambda = (Real*)malloc(k*sizeof(Real)); Real* LLF = (Real*)malloc((maxi+1)*sizeof(Real)); Real* S = (Real*)malloc(lengthA*m*k*sizeof(Real)); - EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau, + EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.,Xa,Y,tau, phiLambda,rhoLambda,piLambda,LLF,S, n,lengthA,m,k); free(Xa); @@ -91,19 +99,16 @@ void constructionModelesLassoMLE_core( free(LLF); free(S); - //~ for j=1:length(a) - //~ phi(a(j),:,:,lambdaIndex) = phiLambda(j,:,:); - //~ end + //Assign results to current variables for (int j=0; j<lengthA; j++) { for (int mm=0; mm<m; mm++) { for (int r=0; r<k; r++) - phi[ai4(a[j],mm,r,lambdaIndex,p,m,k,L)] = phiLambda[ai(j,mm,r,p,m,k)]; + phi[ai4(a[j],mm,r,lambdaIndex,p,m,k,L)] = phiLambda[ai(j,mm,r,lengthA,m,k)]; } } free(phiLambda); - //~ rho(:,:,:,lambdaIndex) = rhoLambda; for (int u=0; u<m; u++) { for (int v=0; v<m; v++) @@ -113,7 +118,6 @@ void constructionModelesLassoMLE_core( } } free(rhoLambda); - //~ pi(:,lambdaIndex) = piLambda; for (int r=0; r<k; r++) pi[mi(r,lambdaIndex,k,L)] = piLambda[r]; free(piLambda); @@ -122,29 +126,24 @@ void constructionModelesLassoMLE_core( int* b = (int*)malloc(m*sizeof(int)); for (int j=0; j<p; j++) { - //~ b = A2(j,2:end,lambdaIndex); - //~ b(b==0) = []; + //b = A2[j,2:dim(A2)[2],lambdaIndex] ; b = b[b!=0] int lengthB = 0; for (int mm=0; mm<m; mm++) { if (A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) b[lengthB++] = A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] - 1; } - //~ if length(b) > 0 - //~ phi(A2(j,1,lambdaIndex),b,:,lambdaIndex) = 0.0; - //~ end if (lengthB > 0) { + //phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0. for (int mm=0; mm<lengthB; mm++) { for (int r=0; r<k; r++) - phi[ai4( A2[ai(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.; + phi[ai4(A2[ai(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.; } } - //~ c = A1(j,2:end,lambdaIndex); - //~ c(c==0) = []; - //~ dimension = dimension + length(c); + //c = A1[j,2:dim(A1)[2],lambdaIndex] ; dimension = dimension + sum(c!=0) for (int mm=0; mm<m; mm++) { if (A1[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) @@ -162,11 +161,6 @@ void constructionModelesLassoMLE_core( Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); for (int i=0; i<n; i++) { - //~ for r=1:k - //~ delta = Y(i,:)*rho(:,:,r,lambdaIndex) - (X(i,a)*(phi(a,:,r,lambdaIndex))); - //~ densite(i,lambdaIndex) = densite(i,lambdaIndex) +... - //~ pi(r,lambdaIndex)*det(rho(:,:,r,lambdaIndex))/(sqrt(2*PI))^m*exp(-dot(delta,delta)/2.0); - //~ end for (int r=0; r<k; r++) { //compute det(rho(:,:,r,lambdaIndex)) [TODO: avoid re-computations] @@ -193,14 +187,16 @@ void constructionModelesLassoMLE_core( for (int v=0; v<lengthA; v++) XiPhiR[u] += X[mi(i,a[v],n,p)] * phi[ai4(a[v],u,r,lambdaIndex,p,m,k,L)]; } - // On peut remplacer X par Xa dans ce dernier calcul, mais je ne sais pas si c'est intéressant ... + // NOTE: On peut remplacer X par Xa dans ce dernier calcul, + // mais je ne sais pas si c'est intéressant ... // compute dotProduct < delta . delta > Real dotProduct = 0.0; for (int u=0; u<m; u++) dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); - densite[mi(lambdaIndex,i,L,n)] += (pi[mi(r,lambdaIndex,k,L)]*detRhoR/pow(sqrt(2.0*M_PI),m))*exp(-dotProduct/2.0); + densite[mi(lambdaIndex,i,L,n)] += + (pi[mi(r,lambdaIndex,k,L)]*detRhoR/pow(sqrt(2.0*M_PI),m))*exp(-dotProduct/2.0); } sumLogDensit += log(densite[lambdaIndex*n+i]); } diff --git a/src/sources/constructionModelesLassoRank.c b/src/sources/constructionModelesLassoRank.c index 59be2f7..1115311 100644 --- a/src/sources/constructionModelesLassoRank.c +++ b/src/sources/constructionModelesLassoRank.c @@ -55,6 +55,8 @@ void constructionModelesLassoRank_core( //Initialize phi to zero, because unactive variables won't be assigned for (int i=0; i<p*m*k*L*Size; i++) phi[i] = 0.0; + for (int i=0; i<L*Size*2; i++) + llh[i] = INFINITY; //initiate parallel section int lambdaIndex; diff --git a/src/sources/utils.h b/src/sources/utils.h index f8cf59e..f6d122f 100644 --- a/src/sources/utils.h +++ b/src/sources/utils.h @@ -37,10 +37,6 @@ typedef double Real; #define ai4(i,j,k,m,d1,d2,d3,d4)\ m*d1*d2*d3 + k*d1*d2 + j*d1 + i -// Array5 Index ; TODO? ... -#define ai5(i,j,k,m,n,d1,d2,d3,d4,d5)\ - n*d1*d2*d3*d4 + m*d1*d2*d3 + k*d1*d2 + j*d1 + i - /************************* * Array copy & "zeroing" ************************/ diff --git a/src/test/Makefile b/src/test/Makefile index 136b1d2..d582826 100644 --- a/src/test/Makefile +++ b/src/test/Makefile @@ -33,7 +33,7 @@ test.selectionTotale: $(LIB) test.selectionTotale.o test_utils.o $(CC) -fPIC -o $@ -c $< $(CFLAGS) $(INCLUDES) clean: - rm -f *.o ../sources/*.o + rm -f *.o ../sources/*.o ../adapters/*.o cclean: clean rm -f $(LIB) $(TESTS) diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R index 95aaf1e..210fe3c 100644 --- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R +++ b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R @@ -34,7 +34,7 @@ generateRunSaveTest_constructionModelesLassoMLE = function(n=200, p=15, m=10, k= row.names=F, col.names=F) write.table(as.double(gamma), paste(testFolder,"gamma",sep=""), row.names=F, col.names=F) - write.table(as.double(lambda), paste(testFolder,"lambda",sep=""), + write.table(as.double(glambda), paste(testFolder,"glambda",sep=""), row.names=F, col.names=F) write.table(as.double(xy$X), paste(testFolder,"X",sep=""), row.names=F, col.names=F) diff --git a/src/test/generate_test_data/helpers/EMGLLF.R b/src/test/generate_test_data/helpers/EMGLLF.R index 7c80081..a624156 100644 --- a/src/test/generate_test_data/helpers/EMGLLF.R +++ b/src/test/generate_test_data/helpers/EMGLLF.R @@ -135,7 +135,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau) { #FIXME: numerical problems, because 0 < det(Rho[,,r] < EPS; what to do ?! # consequence: error in while() at line 77 - Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r] + shift) * det(rho[,,r]) + Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r] + shift) #* det(rho[,,r]) sumLLF1 = sumLLF1 + Gam[i,r] / (2*base::pi)^(m/2) } sumLogLLF2 = sumLogLLF2 + log(sumLLF1) diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R index 359bada..ed05b2a 100644 --- a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R +++ b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R @@ -1,58 +1,56 @@ -constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2){ - #get matrix sizes +constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda, + X,Y,seuil,tau,A1,A2) +{ n = dim(X)[1]; p = dim(phiInit)[1] m = dim(phiInit)[2] - k = dim(phiInit)[3] + k = dim(phiInit)[3] L = length(glambda) - + #output parameters phi = array(0, dim=c(p,m,k,L)) rho = array(0, dim=c(m,m,k,L)) - Pi = matrix(0, k, L) + pi = matrix(0, k, L) llh = matrix(0, L, 2) #log-likelihood - for(lambdaIndex in 1:L){ - a = A1[, 1, lambdaIndex] - a[a==0] = c() - if(length(a)==0){ + for(lambdaIndex in 1:L) + { + a = A1[,1,lambdaIndex] + a = a[a!=0] + if(length(a)==0) next - } - EMGLLf = EMGLLF(phiInit[a,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0,X[,a],Y,tau) - - phiLambda = EMGLLf$phi - rhoLambda = EMGLLf$rho - piLambda = EMGLLf$Pi - - for(j in 1:length(a)){ - phi[a[j],,,lambdaIndex] = phiLambda[j,,] - } - rho[,,,lambdaIndex] = rhoLambda - Pi[,lambdaIndex] = piLambda - + + res = EMGLLF(phiInit[a,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0.,X[,a],Y,tau) + + for (j in 1:length(a)) + phi[a[j],,,lambdaIndex] = res$phi[j,,] + rho[,,,lambdaIndex] = res$rho + pi[,lambdaIndex] = res$pi + dimension = 0 - for(j in 1:p){ - vec = c(2, dim(A2)[2]) - b = A2[j,vec,lambdaIndex] - b[b==0] = c() - if(length(b) > 0){ - phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0 - } - c = A1[j,vec,lambdaIndex] - c[c==0] = c() - dimension = dimension + length(c) + for (j in 1:p) + { + b = A2[j,2:dim(A2)[2],lambdaIndex] + b = b[b!=0] + if (length(b) > 0) + phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0. + c = A1[j,2:dim(A1)[2],lambdaIndex] + dimension = dimension + sum(c!=0) } - + #on veut calculer l'EMV avec toutes nos estimations - densite = matrix(0, n, L) - for(i in 1:n){ - for( r in 1:k){ + densite = matrix(0, nrow=n, ncol=L) + for (i in 1:n) + { + for (r in 1:k) + { delta = Y[i,]%*%rho[,,r,lambdaIndex] - (X[i,a]%*%phi[a,,r,lambdaIndex]); - densite[i,lambdaIndex] = densite[i,lambdaIndex] + Pi[r,lambdaIndex]*det(rho[,,r,lambdaIndex])/(sqrt(2*pi))^m*exp(-tcrossprod(delta)/2.0) + densite[i,lambdaIndex] = densite[i,lambdaIndex] + pi[r,lambdaIndex] * + det(rho[,,r,lambdaIndex])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0) } } llh[lambdaIndex,1] = sum(log(densite[,lambdaIndex])) llh[lambdaIndex,2] = (dimension+m+1)*k-1 } - return(list(phi=phi, rho=rho, Pi=Pi, llh = llh)) + return (list("phi"=phi, "rho"=rho, "pi"=pi, "llh" = llh)) } diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m deleted file mode 100644 index 3e852a6..0000000 --- a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m +++ /dev/null @@ -1,58 +0,0 @@ -function[phi,rho,pi,llh] = constructionModelesLassoMLE(... - phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2) - - PI = 4.0 * atan(1.0); - - %get matrix sizes - n = size(X, 1); - [p,m,k] = size(phiInit); - L = length(glambda); - - %output parameters - phi = zeros(p,m,k,L); - rho = zeros(m,m,k,L); - pi = zeros(k,L); - llh = zeros(L,2); - - for lambdaIndex=1:L - % Procedure Lasso-MLE - a = A1(:,1,lambdaIndex); - a(a==0) = []; - if length(a) == 0 - continue; - end - [phiLambda,rhoLambda,piLambda,~,~] = EMGLLF(... - phiInit(a,:,:),rhoInit,piInit,gamInit,mini,maxi,gamma,0,X(:,a),Y,tau); - - for j=1:length(a) - phi(a(j),:,:,lambdaIndex) = phiLambda(j,:,:); - end - rho(:,:,:,lambdaIndex) = rhoLambda; - pi(:,lambdaIndex) = piLambda; - - dimension = 0; - for j=1:p - b = A2(j,2:end,lambdaIndex); - b(b==0) = []; - if length(b) > 0 - phi(A2(j,1,lambdaIndex),b,:,lambdaIndex) = 0.0; - end - c = A1(j,2:end,lambdaIndex); - c(c==0) = []; - dimension = dimension + length(c); - end - - %on veut calculer l'EMV avec toutes nos estimations - densite = zeros(n,L); - for i=1:n - for r=1:k - delta = Y(i,:)*rho(:,:,r,lambdaIndex) - (X(i,a)*(phi(a,:,r,lambdaIndex))); - densite(i,lambdaIndex) = densite(i,lambdaIndex) +... - pi(r,lambdaIndex)*det(rho(:,:,r,lambdaIndex))/(sqrt(2*PI))^m*exp(-dot(delta,delta)/2.0); - end - end - llh(lambdaIndex,1) = sum(log(densite(:,lambdaIndex))); - llh(lambdaIndex,2) = (dimension+m+1)*k-1; - end - -end diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoRank.R b/src/test/generate_test_data/helpers/constructionModelesLassoRank.R index 15305d9..9254473 100644 --- a/src/test/generate_test_data/helpers/constructionModelesLassoRank.R +++ b/src/test/generate_test_data/helpers/constructionModelesLassoRank.R @@ -45,5 +45,5 @@ constructionModelesLassoRank = function(pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rang } } } - return (list(phi=phi, llh = llh)) + return (list("phi"=phi, "llh" = llh)) } diff --git a/src/test/test.ConstructionModelesLassoMLE.c b/src/test/test.constructionModelesLassoMLE.c similarity index 88% rename from src/test/test.ConstructionModelesLassoMLE.c rename to src/test/test.constructionModelesLassoMLE.c index 45402db..b607467 100644 --- a/src/test/test.ConstructionModelesLassoMLE.c +++ b/src/test/test.constructionModelesLassoMLE.c @@ -1,5 +1,8 @@ #include "constructionModelesLassoMLE.h" #include "test_utils.h" +#include <stdlib.h> + +#include <stdio.h> int main(int argc, char** argv) { @@ -39,7 +42,7 @@ int main(int argc, char** argv) ///////////////////////////////////////// // Call to constructionModelesLassoMLE // - constructionModelesLassoMLE( + constructionModelesLassoMLE_core( phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2, phi,rho,pi,llh, n,p,m,k,L); @@ -56,22 +59,22 @@ int main(int argc, char** argv) free(glambda); // Compare to reference outputs - Real* ref_phi = readArray_real("phi",dimPhi,4); + Real* ref_phi = readArray_real("phi"); compareArray_real("phi", phi, ref_phi, p*m*k*L); free(phi); free(ref_phi); - Real* ref_rho = readArray_real("rho",dimRho,4); + Real* ref_rho = readArray_real("rho"); compareArray_real("rho", rho, ref_rho, m*m*k*L); free(rho); free(ref_rho); - Real* ref_pi = readArray_real("pi",dimPi,2); + Real* ref_pi = readArray_real("pi"); compareArray_real("pi", pi, ref_pi, k*L); free(pi); free(ref_pi); - Real* ref_llh = readArray_real("llh",dimllh,2); + Real* ref_llh = readArray_real("llh"); compareArray_real("llh", llh, ref_llh, L*2); free(llh); free(ref_llh); -- 2.44.0