clear; clc; ntime = 40; nsample = 500; ntest = 50; nq = 3; np = 4; ng =2; nsen = ng*np; changepoint = [0,ntime/2,ntime]; changepointset = [19,20,21]; repnum = 20; %% generate Y l1 = 0.01*log(nsample); l2 = 0.1*log(nsample).*sqrt(nsample); stdnos = 0.025; % construct the bases for every segment TrueA = zeros(ntime,nq*ng); for c = 1:length(changepoint)-1 ind = changepoint(c)+1:changepoint(c+1); tmplen = length(ind); sele = 1:2:2*(nq-1)+1; knot = max(sele)+2; % [psi,xval] = wavefun('morl',10); tmp = bspline_basismatrix(2,linspace(0,1,knot),linspace(0,1,tmplen+2)); % tmp = bspline_basismatrix(3,linspace(0,1,nq+3),linspace(0,1,ntime)); tmp1 = tmp(:,sele); tmp1(1,:) = []; tmp1(end,:) = []; TrueA1 = tmp1; % construct A fourier TrueA2 = fourierbasis(nq,tmplen); % construct wavelet bas % FA = WTortho(tmplen,'Vaidyanathan',4,2); % FA = WTortho(16,'Haar',32,4); % TrueA3 = FA(1:nq,:)'; % tmpTrueA = [TrueA1,TrueA2,TrueA3]; tmpTrueA = [TrueA1,TrueA2]; for i = 1:nq*ng tmpTrueA(:,i) = tmpTrueA(:,i)/norm(tmpTrueA(:,i)); % subplot(2,3,i); % plot(TrueA(:,i)); end TrueA(ind,:) = tmpTrueA; end CorrA = corr(TrueA); % imagesc(CorrA) torthI = TrueA'*TrueA; % for i = 1:ngroup % % segbase = RandOrthMat(seglen); % indr = (i-1)*seglen+1:i*seglen; % indc = (i-1)*qlen+1:i*qlen; % % A(indr,indc) = segbase(:,1:qlen); % tmp = bspline_basismatrix(3,linspace(0,1,knot),linspace(0,1,length(indr))); % A(indr,indc) = bspline_basis(,3,[0:nknot,1],linspace(0,1,length(indr))); % end icmean = zeros(1,np); covMat = eye(np); % covQ = eye(4,3); % for k = 1 : np-1 % for j = (k+1) : np % covQ(j,k) = 0.1^(sqrt(abs(k-j))); % covQ(k,j) = covQ(j,k); % end % end covQ = [1,0.5;0.5,1]; % covQ = 1; cind = [4,3,2,1;3,4,1,2;2,1,4,3]; TrueB = cell(length(changepoint)-1,3); ncause = 2; covcause = covQ; Y = zeros(nsample,ntime,nsen); rootcaus = cell(length(changepoint)-1,ng); Trueautocor = eye(ntime,ntime); stepauto = cell(1,length(changepoint)-1); for c = 1:length(changepoint)-1 ind = changepoint(c)+1:changepoint(c+1); steplen = length(ind); stepauto{c} = eye(steplen,steplen); for k = 1 : steplen-1 for j = (k+1) : steplen stepauto{c}(j,k) = 0.5^(sqrt(abs(k-j))); stepauto{c}(k,j) = stepauto{c}(j,k); end end end % get every sensor's neighbour for i = 1:ng for j = 1:np neigh{np*(i-1)+j} = np*(i-1)+1:np*i; end end invTrueautocor = Trueautocor^-0.5; aautocorr = cell(1,repnum); BBBB = cell(1,repnum); sumresid = zeros(repnum,nsen); ffalse = zeros(repnum,nsen); ffalsec = zeros(repnum,nsen); failc = zeros(repnum,nsen); %parpool('local',12) for iter =1:repnum for c = 1:length(changepoint)-1 ind = changepoint(c)+1:changepoint(c+1); for j = 1:ng indq = (j-1)*nq+1:j*nq; indp = (j-1)*np+1:j*np; % flag = 0; % while flag==0 % tmp1 = zeros(np,np); % for k = 1:np % tmp1(k,cind(c,k)) = 1; % end % TrueB{c,j} = tmp1; % TrueB{c,j} = ones(np,np)*1/(np-1); % TrueB{c,j} = TrueB{c,j} - diag(diag(TrueB{c,j})); % TrueB{c,j} = TrueB{c,j}./repmat(sum(TrueB{c,j},2),1,np); tmp = randn(ncause,np); rootcaus{c,j} = tmp./repmat(sqrt(sum(tmp.^2)),ncause,1); % [V,Q] = eig(TrueB{c,j}); % indd = find(diag(Q)==1); % Vcand = V(:,indd); % end for s = 1:nsample Q = mvnrnd(zeros(1,ncause),covcause,nq)*rootcaus{c,j}; Y(s,ind,indp) = TrueA(ind,indq)*Q; end for s = 1:ntest Q = mvnrnd(zeros(1,ncause),covcause,nq)*rootcaus{c,j}; Ytest(s,ind,indp) = TrueA(ind,indq)*Q; end end end %% calculate the theoritical correlation % covvY = cell(length(changepoint)-1,ng); % invcovY = cell(length(changepoint)-1,ng); % BB = cell(length(changepoint)-1,ng); % % for c = 1:length(changepoint)-1 % for g = 1:ng % covvY{c,g} = rootcaus{c,g}'*covcause*rootcaus{c,g}; % end % end % reshape tureB into matrix form % TrueBBB = zeros(nsen,nsen,ntime); % for c = 1:length(changepoint)-1 % ind = changepoint(c)+1:changepoint(c+1); % TrueBBB(:,:,ind) = repmat(blkdiag(TrueB{c,:}),[1,1,length(ind)]); % end for s = 1:nsample Y(s,:,:)= squeeze(Y(s,:,:)) + mvnrnd(zeros(1,ntime),Trueautocor,nsen)'.*stdnos; end for s = 1:ntest Ytest(s,:,:)= squeeze(Ytest(s,:,:)) + mvnrnd(zeros(1,ntime),Trueautocor,nsen)'.*stdnos; end % figure('position',[40,100,1600,800]); % for i = 1:nsen % subplot(3,4,i); % for k = 1:10 % plot(Y(k,:,i)); % hold on; % xlim([1 ntime]); % end % end % %% esimate B using all the samples icY = Y; [nIC,ntime,nsen] = size(icY); sennorm = mean(squeeze(sqrt(sum(icY.^2,2))),1); icY = reshape(icY,[],nsen); % newdata3 = bsxfun(@minus, newdata3,mean(newdata3)); icY = icY./repmat(sennorm,[size(icY,1),1]); icY = reshape(icY,nIC,ntime,nsen); data_denoise_ICdd = []; data_denoise_ICddt = []; for i = 1:nIC tmp = squeeze(icY(i,:,:)); tmpdata = []; for p = 1:nsen tmpdata = [tmpdata,diag(tmp(:,p))]; end data_denoise_ICddt = [data_denoise_ICddt;tmpdata]; end % DD = zeros((nsen-1)*(ntime-1),(nsen-1)*ntime); % core = zeros(ntime-1,ntime); for i = 1:ntime-1 core(i,i) = -1; core(i,i+1) = 1; end core = {core}; core = repmat(core,1,nsen-1); DD = blkdiag(core{:}); % % BBB = zeros(nsen,nsen,ntime); nstep = 500; l0 =1; tic; % autocorr = repmat(eye(ntime,ntime),[1,1,nsen]); BBB0 = zeros(nsen,nsen,ntime); % matlabpool open 12; BBB = zeros(nsen,nsen,ntime); parfor i = 1:nsen % curautocov = squeeze(autocorr(:,:,i)); % aa = all(abs(eig(curautocov))>=1e-6); % if aa ==1 % sqrtinvcov = curautocov^-0.5; % else % curautocov = squeeze(autocorr(:,:,i))+1e-6*eye(ntime); % sqrtinvcov = curautocov^-0.5; % end % b0 = squeeze(BBB0(i,:,:)); b0(i,:) = []; b0 = reshape(b0',[],1); tmp = zeros(nsen,ntime); icYY = permute(icY,[2,1,3]); icYY = reshape(icYY,ntime,[]); % icYY = sqrtinvcov*icYY; icYY = reshape(icYY,ntime,nIC,nsen); % data_denoise_ICdd = []; for s = 1:nIC tmp = squeeze(icYY(:,s,:)); tmpdata = []; for p = 1:nsen tmpdata = [tmpdata,diag(tmp(:,p))]; end data_denoise_ICdd = [data_denoise_ICdd;tmpdata]; end % cand0 = 1:nsen; Yi = squeeze(icYY(:,:,i)); Yi = reshape(Yi,[],1); cand = data_denoise_ICdd; index = ntime*(i-1)+1:ntime*i; candtmp = cand; candtmp(:,index) = []; cand0 = 1:nsen; cand0(i) = []; % tic; % [BB,totallos,bbb] = fusslasso(Yi,candtmp,l1,l2,l0,nstep,nsen,b0,DD); % toc; l00 = norm(candtmp)^2; tic; [BB,totallos,bbb] = fusslasso1(Yi,candtmp,l1,l2,l00,nstep,nsen,b0,DD); toc; aa = reshape(BB,ntime,nsen-1); ttt = zeros(nsen,ntime); ttt(cand0,:) = aa'; BBB(i,:,:) = ttt; pred = candtmp*BB; resid = Yi- pred; sumresid(iter,i) = (norm(resid,2))^2; resid = reshape(resid,ntime,nIC); covresid = cov(resid'); % covresid = diag(diag(cov(resid'))); % ttmp = diag(covresid); % ind = find(ttmp<1e-4); % aa = zeros(1,ntime); % aa(ind) = 1; % covresid = covresid + diag(aa); autocorr(:,:,i) = covresid; end % % test % residd = zeros(ntime,nIC,nsen); % figure('position',[40,100,1200,800]); % for i = 1:nsen % b0 = squeeze(BBB(i,:,:)); % b0 = reshape(b0',[],1); % tmp = zeros(nsen,ntime); % cand0 = 1:nsen; % Yi = squeeze(icY(:,:,i)); % Yi = reshape(Yi',[],1); % cand = data_denoise_ICddt; % candtmp = cand; % cand0 = 1:nsen; % pred = candtmp*b0; % resid = Yi- pred; % resid = reshape(resid,ntime,nIC); % residd(:,:,i) = resid; % % covresid = diag(diag(cov(resid'))); % covresid = cov(resid'); % % ttmp = diag(covresid); % % ind = find(ttmp<1e-4); % % aa = zeros(1,ntime); % % aa(ind) = 1; % % covresid = covresid + 1e-4*diag(aa); % autocorr(:,:,i) = covresid; % % subplot(3,4,i); % plot(icY(1,:,i)); % hold on; % plot(pred(1:ntime),'r'); % xlim([1 ntime]); % end % sumresid(iter,:) = sum(squeeze(sum(residd.^2,1)),1); aautocorr{iter} = autocorr; % figure('position',[40,100,1200,800]); % for i = 1:nsen % subplot(3,4,i); % imagesc(squeeze(autocorr(:,:,i))); % colorbar; % end BBBB{iter} = BBB; %% detect change point (temporal segment) diffB = BBB(:,:,2:end)-BBB(:,:,1:end-1); diffB = squeeze(mean(mean(abs(diffB),1))); ind = find(diffB>3*std(diffB)); if isempty(ind) failc(iter) = 1; else findc = ismember(ind,changepointset); indd = find(findc==0); ffalsec(iter) = length(indd); end identified_changept = [0,ind',ntime]; sss = diff(identified_changept); aaa = find(sss==1); identified_changept(aaa) = []; %% spatio segment num_time_segment = length(identified_changept)-1; false = zeros(num_time_segment,nsen); for s = 1:num_time_segment ss = identified_changept(s)+1; ee = identified_changept(s+1); for i = 1:nsen tmpB = squeeze(BBB(i,:,ss:ee)); tmpBsum = mean(tmpB,2); ind = find(abs(tmpBsum)>1e-2); dec = sum(ismember(ind,neigh{i})); if dec0); % end % % dd = zeros(nsen,ntime); % for i = 1:12 % dd(i,:) = diag(autocorr(:,:,i)); % end % k = 1; for i = 1:nsen figure('position',[40,100,1500,1000]); % subplot(5,7,i); % imagesc(squeeze(TrueBBB(i,:,:))); % colorbar; for k = 1 subplot(1,2,k); imagesc(squeeze(BBB1(i,:,:))); colorbar; end end for i = 1:nsen figure('position',[40,100,1200,800]); % subplot(2,2,1); % imagesc(squeeze(TrueBBB(i,:,:))); % colorbar; for k = 1 subplot(1,1,k); imagesc(squeeze(aautocorr{k}(:,:,i))); colorbar; end end % % aa = mvnrnd(zeros(1,ntime),Trueautocor,400); % imagesc(cov(aa)) % % aa = zeros(ntime,ntime,niter); % for i = 1:niter % aa(:,:,i) = aautocorr{i}(:,:,1); % end % % aa = reshape(aa,[],4); % figure('position',[40,100,1500,1000]); for i = 1:nsen subplot(4,6,i); for k = 1:10 plot(data_denoise_ICt0(k,:,i)); % plot(data_denoise_OCnew(k,:,i)); hold on; xlim([1 ntime]); end end % figure('position',[40,100,1500,1000]); % for i = 1:4 % subplot(2,2,i); % plot(TrueA(:,i)); % xlim([1 ntime]); % endst