{
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;
}
}
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)];
}
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
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;
#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++)
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++)
{
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++)
{
}
}
- //[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);
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++)
}
}
free(rhoLambda);
- //~ pi(:,lambdaIndex) = piLambda;
for (int r=0; r<k; r++)
pi[mi(r,lambdaIndex,k,L)] = piLambda[r];
free(piLambda);
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)
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]
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]);
}
//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;
#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"
************************/
$(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)
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)
{
#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)
-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))
}
+++ /dev/null
-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
}
}
}
- return (list(phi=phi, llh = llh))
+ return (list("phi"=phi, "llh" = llh))
}
#include "constructionModelesLassoMLE.h"
#include "test_utils.h"
+#include <stdlib.h>
+
+#include <stdio.h>
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);
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);
projects / valse.git / commitdiff