0001 function [IDX,TODSS,SCORE]=nt_cluster_jd2(x,dsr,flags,init)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 if nargin<2; error('!'); end
0017 if nargin<3 ||isempty(flags); flags=[]; end
0018 if nargin<4; init=[]; end
0019
0020 if ~nargout; disp('entering nt_cluster_jd...'); end
0021
0022
0023 [xx,ind]=nt_xprod(x,'lower',dsr);
0024
0025
0026 if find(strcmp(flags,'norm'))
0027 xx=nt_normrow(xx);
0028 end
0029 if find(strcmp(flags,'norm2'))
0030 xx=norm2(xx,size(x,2),ind);
0031 end
0032
0033
0034 if isempty(init)
0035 [C,A,score]=nt_cluster1D(xx);
0036 [~,idx]=min(score);
0037 A=A(:,idx);
0038
0039 A=repmat(A',[dsr,1]);
0040 A=A(:);
0041 A(end:size(x,1))=A(end);
0042 IDX{1}=find(A==1);
0043 IDX{2}=find(A==2);
0044 else
0045 IDX{1}=setdiff(1:size(x,1),init);
0046 IDX{2}=init;
0047 end
0048
0049
0050 c0=nt_cov(x);
0051 c1=nt_cov(x(IDX{2},:));
0052 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0053 todss2=todss(:,[1 end]);
0054 z=nt_mmat(x,todss2);
0055
0056 PLOT_FIG2=0;
0057 if PLOT_FIG2
0058 figure(2); clf; set(gcf, 'name','in nt_cluster_jd');
0059 A=zeros(size(x,1),1); A(IDX{1})=1;
0060 subplot 511; plot(x); title('data');
0061 subplot 512; plot(A,'.-'); title('initial cluster map');
0062 subplot 513; plot(z(:,1)); title('initial DSS1');
0063 subplot 514; plot(z(:,2)); title('initial DSS2');
0064 drawnow;
0065 end
0066
0067
0068 old_IDX=IDX;
0069 for k=1:10
0070
0071 [zz,ind]=nt_xprod(z,[],dsr);
0072
0073
0074 if find(strcmp(flags,'norm'))
0075 zz=nt_normrow(zz);
0076 end
0077 if find(strcmp(flags,'norm2'))
0078 zz=norm2(zz,size(z,2),ind);
0079 end
0080
0081 if 0
0082 EXPONENT=0.5;
0083 zz(1,1)=(zz(1,1)).^EXPONENT;
0084 zz(2,2)=(zz(2,2)).^EXPONENT;
0085 end
0086
0087
0088 [C,A,score]=nt_cluster1D(zz);
0089 [~,idx]=min(score);
0090 A=A(:,idx);
0091 A=repmat(A',[dsr,1]);
0092 A=A(:);
0093 A(end:size(x,1))=A(end);
0094 IDX{1}=find(A==1);
0095 IDX{2}=find(A==2);
0096
0097 if ~nargout; disp(['cluster sizes: ', num2str([numel(IDX{1}), numel(IDX{2})])]); end
0098
0099
0100 c1=nt_cov(x(IDX{1},:));
0101 c0=c1+nt_cov(x(IDX{2},:));
0102 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0103 z=nt_mmat(x,todss(:,[1 end]));
0104
0105 if PLOT_FIG2
0106 figure(2);
0107 subplot 515; plot(A,'.-'); title('final cluster map');
0108 end
0109 if all(size(old_IDX{1})==size(IDX{1})) && all(old_IDX{1}==IDX{1}); break; end
0110 old_IDX=IDX;
0111 end
0112
0113
0114 c1=nt_cov(x(IDX{1},:));
0115 c0=c1+nt_cov(x(IDX{2},:));
0116 [TODSS{1},pwr0,pwr1]=nt_dss0(c0,c1);
0117 SCORE(1)=pwr1(end)/pwr0(end);
0118
0119 c1=nt_cov(x(IDX{2},:));
0120 c0=c1+nt_cov(x(IDX{1},:));
0121 [TODSS{2},pwr0,pwr1]=nt_dss0(c0,c1);
0122 SCORE(2)=pwr1(end)/pwr0(end);
0123
0124 SCORE=SCORE(:);
0125
0126 if nargout==0;
0127
0128
0129 disp(['cluster1: ',num2str(100*numel(find(A==1))/numel(A)), '%']);
0130
0131 figure(100); clf
0132 subplot 311;
0133 plot(x); hold on
0134 xx=x;
0135 xx(IDX{1},:)=nan;
0136 plot(xx,'k');
0137 axis tight
0138 title('black: cluster B');
0139
0140 z1=nt_mmat(x,todss(:,1));
0141 z2=nt_mmat(x,todss(:,end));
0142
0143 figure(101); clf ;
0144 subplot 221;
0145 plot(pwr1./pwr0,'.-'); xlabel('component'); ylabel('score'); title('DSS cluster A vs all');
0146 subplot 222;
0147 wsize=min(1024,size(z1,1));
0148 nt_spect_plot(z1,wsize,[],[],1);
0149 hold on
0150 wsize=min(1024,size(z2,1));
0151 nt_spect_plot(z2,wsize,[],[],1);
0152 xlim([0 .5])
0153 nt_linecolors([],[1 3]);
0154 legend('first','last'); legend boxoff
0155 hold off
0156
0157 z=nt_mmat(x,todss);
0158 z=nt_normcol(z);
0159 subplot 223; imagescc(nt_cov(z(IDX{1},:))); title('cluster A');
0160 subplot 224; imagescc(nt_cov(z(IDX{2},:))); title('cluster B');
0161
0162
0163 figure(100);
0164 subplot 312;
0165 plot(z1); axis tight
0166 title('first DSS component')
0167 subplot 313;
0168 plot(z2); axis tight
0169 title('last DSS component');
0170
0171
0172
0173
0174
0175
0176
0177 if 0
0178 figure(105); clf
0179 subplot 211;
0180 nt_sgram(z1,1024,32,[],1);
0181 title('first');
0182 subplot 212;
0183 nt_sgram(z2,1024,32,[],1);
0184 title('last');
0185 end
0186 clear IDX SCORE TODSS
0187
0188 end
0189
0190 function y=norm2(x,n,ind)
0191 [I,J]=ind2sub([n,n],ind);
0192 for k=1:size(x,1)
0193 a=x(k,1:n);
0194 b=sqrt(a(I).*a(J));
0195 y(k,:)=x(k,:)./b;
0196 end
0197
0198
0199
0200