From 46a2e676b2d85eef5a1811a6e623b65327fc451d Mon Sep 17 00:00:00 2001 From: Benjamin Auder 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 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 Real dotProduct = 0.0; for (int u=0; u 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 + +#include 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