fix memory leaks on EMGLLF, test OK for EMGrank
[valse.git] / src / test / generate_test_data / helpers / EMGLLF.R
1 EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau){
2 #matrix dimensions
3 n = dim(X)[1]
4 p = dim(phiInit)[1]
5 m = dim(phiInit)[2]
6 k = dim(phiInit)[3]
7
8 #init outputs
9 phi = phiInit
10 rho = rhoInit
11 Pi = piInit
12 LLF = rep(0, maxi)
13 S = array(0, dim=c(p,m,k))
14
15
16 gam = gamInit
17 Gram2 = array(0, dim=c(p,p,k))
18 ps2 = array(0, dim=c(p,m,k))
19 b = rep(0, k)
20 pen = matrix(0, maxi, k)
21 X2 = array(0, dim=c(n,p,k))
22 Y2 = array(0, dim=c(n,m,k))
23 dist = 0
24 dist2 = 0
25 ite = 1
26 Pi2 = rep(0, k)
27 ps = matrix(0, m,k)
28 nY2 = matrix(0, m,k)
29 ps1 = array(0, dim=c(n,m,k))
30 nY21 = array(0, dim=c(n,m,k))
31 Gam = matrix(0, n,k)
32 EPS = 1E-15
33
34 while(ite <= mini || (ite<= maxi && (dist>= tau || dist2 >= sqrt(tau)))){
35 Phi = phi
36 Rho = rho
37 PI = Pi
38 #calcul associé à Y et X
39 for(r in 1:k){
40 for(mm in 1:m){
41 Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm] ##bon calcul ? idem pour X2 ??...
42 }
43 for(i in 1:n){
44 X2[i,,r] = X[i,] *sqrt(gam[i,r])
45 }
46 for(mm in 1:m){
47 ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
48 }
49 for(j in 1:p){
50 for(s in 1:p){
51 Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
52 }
53 }
54 }
55
56 ##########
57 #Etape M #
58 ##########
59
60 #pour pi
61 for(r in 1:k){
62 b[r] = sum(sum(abs(phi[,,r])))
63 }
64 gam2 = colSums(gam)
65 a = sum(gam%*%(log(Pi)))
66
67 #tant que les props sont negatives
68 kk = 0
69 pi2AllPositive = FALSE
70 while(pi2AllPositive == FALSE){
71 Pi2 = Pi + 0.1^kk * ((1/n)*gam2 - Pi)
72 pi2AllPositive = TRUE
73 for(r in 1:k){
74 if(Pi2[r] < 0){
75 pi2AllPositive = false;
76 break
77 }
78 }
79 kk = kk+1
80 }
81
82 #t[m]la plus grande valeur dans la grille O.1^k tel que ce soit
83 #décroissante ou constante
84 while((-1/n*a+lambda*((Pi^gamma)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(b)) && kk<1000){
85 Pi2 = Pi+0.1^kk*(1/n*gam2-Pi)
86 kk = kk+1
87 }
88 t = 0.1^(kk)
89 Pi = (Pi+t*(Pi2-Pi)) / sum(Pi+t*(Pi2-Pi))
90
91 #Pour phi et rho
92 for(r in 1:k){
93 for(mm in 1:m){
94 for(i in 1:n){
95 ps1[i,mm,r] = Y2[i,mm,r] * (X2[i,,r]%*%(phi[,mm,r]))
96 nY21[i,mm,r] = (Y2[i,mm,r])^2
97 }
98 ps[mm,r] = sum(ps1[,mm,r])
99 nY2[mm,r] = sum(nY21[,mm,r])
100 rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r]))
101 }
102 }
103 for(r in 1:k){
104 p1 = p-1
105 for(j in 1:p1){
106 for(mm in 1:m){
107 j1 = j-1
108 j2 = j+1
109 v1 = c(1:j1)
110 v2 = c(j2:p)
111 S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r]) + phi[v2,mm,r]%*%(Gram2[j,v2,r]) #erreur indice
112 if(abs(S[j,mm,r]) <= n*lambda*(Pi[r]^gamma)){
113 phi[j,mm,r]=0
114 }else{
115 if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)){
116 phi[j,mm,r] = (n*lambda*(Pi[r]^gamma)-S[j,mm,r])/Gram2[j,j,r]
117 }else{
118 phi[j,mm,r] = -(n*lambda*(Pi[r]^gamma)+S[j,mm,r])/Gram2[j,j,r]
119 }
120 }
121 }
122 }
123 }
124
125 ##########
126 #Etape E #
127 ##########
128 sumLogLLF2 = 0
129 for(i in 1:n){
130 #precompute dot products to numerically adjust their values
131 dotProducts = rep(0,k)
132 for(r in 1:k){
133 dotProducts[r] = tcrossprod(Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])
134 }
135 shift = 0.5*min(dotProducts)
136
137 #compute Gam(:,:) using shift determined above
138 sumLLF1 = 0.0;
139 for(r in 1:k){
140 Gam[i,r] = Pi[r]*det(rho[,,r])*exp(-0.5*dotProducts[r] + shift)
141 sumLLF1 = sumLLF1 + Gam[i,r]/(2*pi)^(m/2)
142 }
143 sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
144 sumGamI = sum(Gam[i,])
145 if(sumGamI > EPS)
146 gam[i,] = Gam[i,] / sumGamI
147 else
148 gam[i,] = rep(0,k)
149 }
150
151
152 sumPen = 0
153 for(r in 1:k){
154 sumPen = sumPen + Pi[r]^gamma^b[r]
155 }
156 LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
157
158 if(ite == 1)
159 dist = LLF[ite]
160 else
161 dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
162
163 Dist1=max(max(max((abs(phi-Phi))/(1+abs(phi)))))
164 Dist2=max(max(max((abs(rho-Rho))/(1+abs(rho)))))
165 Dist3=max(max((abs(Pi-PI))/(1+abs(PI))))
166 dist2=max(c(Dist1,Dist2,Dist3))
167
168 ite=ite+1
169 }
170
171 Pi = t(Pi)
172 return(list("phi"=phi, "rho"=rho, "pi"=Pi, "LLF"=LLF, "S"=S))
173 }