summaryrefslogtreecommitdiffstats
path: root/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m')
-rwxr-xr-xSD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m187
1 files changed, 187 insertions, 0 deletions
diff --git a/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m
new file mode 100755
index 0000000..3ca1046
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m
@@ -0,0 +1,187 @@
1% function [v,s] = quickNcutHardBiais(W,U,nbEigenValues,dataGraphCut)
2%ligne 35 : changement tim
3%[v,s] = ncut(W,nbEigenValues,[],offset);
4%devient :
5%[v,s] = tim_ncut_rapide(W,nbEigenValues,[],offset);
6%et eigs devient tim_eigs
7
8% Input:
9% W = affinity matrix
10% U = hard constraint matrix, could be a cell of partial grouping
11% nbEigenValues = number of eigenvectors
12% offset = regularization factor, default = 0
13% Output:
14% v = eigenvector
15% s = eigenvalue of (W,d), s.t. U' * x = 0.
16
17% call eigs using my own * operation
18
19% Stella X. Yu, Jan 2002.
20
21function [v,s] = quickNcutHardBiais2(W,U,nbEigenValues,dataGraphCut)%voir : rajouter sigma
22n = size(W,1);
23nbEigenValues = min(nbEigenValues,n);
24
25offset = dataGraphCut.offset;
26%offset = 2;
27
28% degrees and regularization
29d = sum(abs(W),2);
30dr = 0.5 * (d - sum(W,2));
31d = d + offset * 2;
32dr = dr + offset;
33W = W + spdiags(dr,0,n,n);
34clear dr
35
36% normalize
37Dinvsqrt = 1./sqrt(d+eps);
38P = spmtimesd(W,Dinvsqrt,Dinvsqrt);
39clear W;
40
41% if max(max(P-P')) < 1e-10 %ou eps
42% %S = sparse(1:n,1:n,0.5);
43% P =max(P,P');
44% % P=S*(P+P');
45% %P=0.5*(P+P');
46% options.issym = 1;
47% end
48P = sparsifyc(P,dataGraphCut.valeurMin);
49options.issym = 1;
50
51Ubar = spmtimesd(U,Dinvsqrt,[]);
52%Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin); %voir
53
54options.disp = dataGraphCut.verbose;
55options.maxit = dataGraphCut.maxiterations;
56options.tol = dataGraphCut.eigsErrorTolerance;
57
58options.v0 = ones(size(P,1),1);%voir
59
60options.p = max(35,2*nbEigenValues); %voir
61options.p = min(options.p , n);
62
63% nouvelle idee : factorisation de Cholesky
64C=Ubar'*Ubar;
65%permutation = symamd(C);
66%R = cholinc(C(permutation,permutation));
67t_chol_Ubar = cputime;
68[R,ooo] = cholinc(C,'0');
69t_chol_Ubar = cputime - t_chol_Ubar;
70%if error occurs, check if Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin);
71%sparsifies too much
72
73
74% compute H = (Ubar'*Ubar)^(-1)
75% t_inv_H = cputime;
76% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer
77% t_inv_H = cputime - t_inv_H;
78% H = sparsifyc(H,dataGraphCut.valeurMin);
79% tEigs = cputime;
80% if options.issym & max(max(H-H')) < 1e-10
81% [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H));
82% else
83% [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H);
84% end
85% tEigs = cputime - tEigs;
86%
87
88
89
90R = sparsifyc(R,dataGraphCut.valeurMin);
91tEigs = cputime;
92if options.issym
93 [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,tril(P),Ubar,R);
94else
95 [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R);
96end
97tEigs = cputime - tEigs;
98
99
100%afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs),dataGraphCut.verbose,2);
101%afficheTexte(sprintf('\nTemps ecoule pendant chol(Ubar''*Ubar) : %g',t_chol_Ubar),dataGraphCut.verbose);
102if convergence~=0
103 afficheTexte(sprintf(' (Non-convergence)'),dataGraphCut.verbose);
104end
105
106
107%disp(sprintf('nnz(P) : %f\n',nnz(P)));
108%disp(sprintf('nnz(Ubar) : %f\n',nnz(Ubar)));
109%disp(sprintf('nnz(R) : %f\n',nnz(R)));
110%disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 4*nnz(R)));
111
112
113
114s = real(diag(s));
115[x,y] = sort(-s);
116s = -x;
117vbar = vbar(:,y);
118
119
120v = spdiags(Dinvsqrt,0,n,n) * vbar;
121
122for i=1:size(v,2)
123 %v(:,i) = v(:,i) / max(abs(v(:,i)));
124 v(:,i) = (v(:,i) / norm(v(:,i)) )*norm(ones(n,1));
125 if v(1,i)~=0
126 v(:,i) = - v(:,i) * sign(v(1,i));
127 end
128end
129
130% % nouvelle idee : factorisation de Cholesky
131% t_chol_Ubar = cputime;
132% R = chol(Ubar' * Ubar);
133% t_chol_Ubar = cputime - t_chol_Ubar;
134% R = sparsifyc(R,dataGraphCut.valeurMin);
135% tEigs = cputime;
136% if options.issym
137% [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R);
138% else
139% [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R);
140% end
141% tEigs = cputime - tEigs;
142
143
144% % compute H = (Ubar'*Ubar)^(-1)
145% t_inv_H = cputime;
146% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer
147% t_inv_H = cputime - t_inv_H;
148% H = sparsifyc(H,dataGraphCut.valeurMin);
149% tEigs = cputime;
150% if options.issym & max(max(H-H')) < 1e-10
151% [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H));
152% else
153% [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H);
154% end
155% tEigs = cputime - tEigs;
156
157
158
159% % idee de mon rapport... semble pas marcher car R = qr(Ubar,0) est plus
160% % lent que H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin));
161% t_qr_Ubar = cputime;
162% R = qr(Ubar,0);
163% t_qr_Ubar = cputime - t_qr_Ubar;
164% R = sparsifyc(R,dataGraphCut.valeurMin);
165% tEigs2 = cputime;
166% if options.issym
167% [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R);
168% else
169% [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R);
170% end
171% tEigs2 = cputime - tEigs2;
172
173
174
175% idee de Jianbo... semble pas marcher car on a besoin de prendre k maximal
176% dans [A,S,B] = svds(Ubar,k);
177%
178% [A,S,B] = svds(Ubar,300);
179% A = sparsifyc(A,dataGraphCut.valeurMin);
180% tEigs = cputime;
181% [vbar,s,convergence] = tim_eigs(@mex_projection_svd,n,nbEigenValues,'lm',options,P,A);
182
183% afficheTexte(sprintf('\ninv(H) : %g',t_inv_H),dataGraphCut.verbose);
184% afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs2),dataGraphCut.verbose,2);
185% afficheTexte(sprintf('\nqr(Ubar) : %g',t_qr_Ubar),dataGraphCut.verbose);
186% disp(sprintf('nnz(H) : %f\n',nnz(H)));
187% disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 2*nnz(H)));
generated by cgit v1.2.2 (