demo_multi_cluster_s.m

load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\conv_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\activity_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\dark_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\phonecharge_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\phonelock_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\audio_feature.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\pre_PHQ9.mat'
load 'C:\Users\asma\Desktop\Jins ALG\Depression_2\data\post_PHQ9.mat'

prefix='C:\Users\asma\Desktop\Jins ALG\Depression_2\';
addpath([prefix '\utility_code_1\clustering\ssvd-code']);
addpath([prefix '\utility_code_1\clustering\spectralCoClustering']);
addpath([prefix '\utility_code_1\clustering']);
addpath([prefix '\utility_code_1\short_functions']);
addpath([prefix '\utility_code_1\InfoTheory']);
addpath([prefix '\utility_code_1\kernel']);
addpath([prefix '\utility_code_1\evaluate']);
addpath([prefix '\utility_code_1\code']);
addpath([prefix '\code\exsiting_algorithm']);
addpath([prefix '\code']);
addpath([prefix '\proximal_2']);

%% read the raw data for activity, conversation, dark, audio, phonelock records and give the average on day values as view1. Each row of view1 is an instance from the population.
% avaerage one user/day
view1(:, 1) = average_c(activity_feature, 3); % for all 49 user get 1 value
view1(:, 2) = average_c(conv_feature, 3);
view1(:, 3) = average_c(dark_feature, 3);
view1(:, 4) = average_c(dark_feature, 4);
view1(:, 5) = average_c(conv_feature, 2);
view1(:, 6) = average_c(activity_feature, 4);
view1(:, 7) = average_c(activity_feature, 5);
view1(:, 8) = average_c(audio_feature, 3);
view1(:, 9) = average_c(audio_feature, 4);
view1(:, 10) = average_c(audio_feature, 5);
view1(:, 11) = average_c(phonelock_feature, 2);
view1(:, 12) = average_c(phonelock_feature, 3);


%% read the sine parameters which fit the time series data of activity, audio, conversation, phonelock as view2. 
view2(:, 1:4) = denoising_wl_sin(activity_feature, 3, [500 20 3 700]); % amp, phase, intercept and freq ///for all 49 user get 4 value
view2(:, 5:8) = denoising_wl_sin(audio_feature, 3, [500 20 3 700]);
view2(:, 9:12) = denoising_wl_sin(conv_feature, 3, [500 20 3 700]);
% view2(:, 13:16) = denoising_wl_sin(phonelock_feature, 3);
% view2 = [sin_act_5, sin_audio_5, sin_conv_5, sin_lock_5];
view1_data = view1;
view2_data = view2;

% normalize the data.
n = size(view1_data, 1);
d1 = size(view1_data, 2);
d2 = size(view2_data, 2);

M_1_norm = normc(view1_data);
M_2_norm = normc(view2_data);

M_1 = M_1_norm;
M_2 = M_2_norm;

M = cell(1, 1);
M{1} = M_1;
M{2} = M_2;

%% using cv to find the proper parameters lambda_z, lambda_2 then run 
% Multiview Biclustering method to get the 1st cluster. lambda_z is for the
% size of the cluster we want. lambda_2 is for the number of features to 
% differentiate the cluster from the rest of the population. 
lambda_z = 9;
lambda_2 = 5;
iSeedV1 = 12;
%% z1 is a binary vector as a idnetifier of our explored cluster. 1 means
% the ith instance on the ith row belong to the identified cluster. 0 means
% this instance does not belong to this cluster. V is a matrix shows the 
% significantly useful features that be detected to diffentiate the 
% identified cluster Clus1 from the rest by our method. 


% Depressed people identification, also find what features are important
% for our clustering method
% cluster people into group, and see how they are related to depression
[z1, U, V, obj] = proxi3_3(M, lambda_z, [lambda_2; lambda_2], iSeedV1);
rowClus1 = double(z1~=0);
Clus1 = rowClus1; %After finding one cluster, erase that from data and then move forward for next clustering

%% erase the instances in Clus1, then prepare for the next clustering process.
IND = find(z1~=0);
M2_1 = M_1;
M2_2 = M_2;
M2_1(IND, :)=[];
M2_2(IND, :)=[];
M2{1} = M2_1;
M2{2} = M2_2;

%% run multiview biclustering method again to detect the 2nd cluster Clus2.
lambda_z2 = 7;
iSeedV2 = 1;
n_1_pre = [];
n_2_pre = [];
n_3_pre = [];
n_1_post = [];
n_2_post = [];
n_3_post = [];
[z2, U2, V2, obj] = proxi3_3(M2, lambda_z2, [lambda_2; lambda_2], iSeedV2);
rowClus2 = double(z2~=0);
Clus2 = zeros(n, 1);
Clus2(IND) = 0;
IND2 = 1:1:49;
IND2(IND) = [];
TN = find(z2==0);
IND2(TN) = [];
Clus2(IND2) = 1;
Clus3 = ones(n, 1) - Clus1 - Clus2;

% view 1 is ok
% view 2, wavelet filtering, f transformation -> raw data...can improve ft
% method. using ft to find most important point. Also we can add useful
% descrption - more features
%------------------------------------
phqmat = -ones(60,5);
    for i = 1: 49
        te = audio_feature{1,i};
        id(i) = te(1,1);
    end
    for i = 1 :46
        phqmat(pre_PHQ9(i, 1)+1, 1) = pre_PHQ9(i, 2);
    end
    for i = 1 :38
        phqmat(post_PHQ9(i, 1)+1, 2) = post_PHQ9(i, 2);
    end
    for i = 1 :49
        phqmat(id(i)+1, 3) = Clus1(i);
    end
    for i = 1 :49
        phqmat(id(i)+1, 4) = Clus2(i);
    end
    for i = 1 :49
        phqmat(id(i)+1, 5) = Clus3(i);
    end
    sum_1_pre = 0;
    t =0;
    for i = 1 : 60
        if(phqmat(i,3)==1 && phqmat(i, 1)>0)
           t = t+1;
           sum_1_pre = sum_1_pre + phqmat(i, 1);
           n_1_pre(t) = phqmat(i, 1);
        end
    end
    sum_1_pre = sum_1_pre/t;

    t=0;
    sum_1_post = 0;
    for i = 1 : 60
        if(phqmat(i,3)==1 && phqmat(i, 2)>0)
            t = t+1;
            sum_1_post = sum_1_post + phqmat(i, 2);
            n_1_post(t) = phqmat(i, 2);
        end
    end
    sum_1_post = sum_1_post/t;
    
    t=0;
    sum_2_pre = 0;
    for i = 1 : 60
        if(phqmat(i,4)==1 && phqmat(i, 1)>0)
        t = t+1;
        sum_2_pre = sum_2_pre + phqmat(i, 1);
        n_2_pre(t) = phqmat(i, 1);
        end
    end
    sum_2_pre = sum_2_pre/t;
    
    t=0;
    sum_2_post = 0;
    for i = 1 : 60
        if(phqmat(i,4)==1 && phqmat(i, 2)>0)
            t = t+1;
            sum_2_post = sum_2_post + phqmat(i, 2);
            n_2_post(t) = phqmat(i, 2);
        end
    end
sum_2_post = sum_2_post/t;

    
    t=0;
    sum_3_pre = 0;
    for i = 1 : 60
        if(phqmat(i,5)==1 && phqmat(i, 1)>0)
        t = t+1;
        sum_3_pre = sum_3_pre + phqmat(i, 1);
        n_3_pre(t) = phqmat(i, 1);
        end
    end
    sum_3_pre = sum_3_pre/t;
    
    t=0;
    sum_3_post = 0;
    for i = 1 : 60
        if(phqmat(i,5)==1 && phqmat(i, 2)>0)
            t = t+1;
            sum_3_post = sum_3_post + phqmat(i, 2);
            n_3_post(t) = phqmat(i, 2);
        end
    end
sum_3_post = sum_3_post/t;

var_1_pre = std(n_1_pre);
var_2_pre = std(n_2_pre);
var_3_pre = std(n_3_pre);
var_1_post = std(n_1_post);
var_2_post = std(n_2_post);
var_3_post = std(n_3_post);
var_pre = std(pre_PHQ9(:, 2));
var_post = std(post_PHQ9(:, 2));
avar_1_pre = var_1_pre/var_pre ;
avar_2_pre = var_2_pre/var_pre;
avar_3_pre = var_3_pre/var_pre;
avar_1_post = var_1_post/var_post;
avar_2_post = var_2_post/var_post;
avar_3_post = var_3_post/var_post;

sum_pre = sum(pre_PHQ9(:, 2))/46;
sum_post = sum(post_PHQ9(:, 2))/38;
aver_1_pre = sum_1_pre / sum_pre;%sum_1_pre = sum(where clus=1 & score >0 /total such instances)
aver_2_pre = sum_2_pre / sum_pre;
aver_3_pre = sum_3_pre / sum_pre;
aver_1_post = sum_1_post / sum_post;
aver_2_post = sum_2_post / sum_post;
aver_3_post = sum_3_post / sum_post;

vari1 = zeros(3, 14);
vari1(1,1) = avar_1_pre/ sum_pre; 
vari1(1,2) = avar_1_post/ sum_post;
vari1(2,1) = avar_2_pre/ sum_pre; 
vari1(2,2) = avar_2_post/ sum_post;
vari1(3,1) = avar_3_pre/ sum_pre; 
vari1(3,2) = avar_3_post/ sum_post;

vari2 = zeros(3, 18);
vari2(1,1) = avar_1_pre/ sum_pre; 
vari2(1,2) = avar_1_post/ sum_post;
vari2(2,1) = avar_2_pre/ sum_pre; 
vari2(2,2) = avar_2_post/ sum_post;
vari1(3,1) = avar_3_pre/ sum_pre; 
vari1(3,2) = avar_3_post/ sum_post;

Clusmat = [Clus1, Clus2, Clus3];
View1_perc = Clusmat' * view1;
View2_perc = Clusmat' * view2;
View1_perc(1, :) = View1_perc(1, :) / sum(Clus1);
View1_perc(2, :) = View1_perc(2, :) / sum(Clus2);
View1_perc(3, :) = View1_perc(3, :) / sum(Clus3);
View2_perc(1, :) = View2_perc(1, :) / sum(Clus1);
View2_perc(2, :) = View2_perc(2, :) / sum(Clus2);
View2_perc(3, :) = View2_perc(3, :) / sum(Clus3);
Clusall = double(Clus1|Clus2|Clus3);
averView1 = (Clusall' * view1) / sum(Clusall);
averView2 = (Clusall' * view2) /  sum(Clusall);
PHQ = [aver_1_pre,aver_1_post;aver_2_pre,aver_2_post;aver_3_pre,aver_3_post];
for i = 1 : 3
View1_perc(i, :) = View1_perc(i, :) ./ averView1;
View2_perc(i, :) = View2_perc(i, :) ./ averView2;
end
TT = [PHQ, View1_perc];
figure; bar(TT,1);

set(gca,'XTickLabel',{'Cluster1','Cluster2','Cluster3'})
xlabel('View-Meanvalue','FontSize',16),ylabel('Relative distance to the average','FontSize',16)
legend('PRE-PHQ','POST-PHQ','Act_1', 'Conv_t', 'Dark_n', 'Dark_t', 'Conv_n', 'Act_2', 'Act_3', 'Aud_1', 'Aud_2','Aud_3','Loc_n','Loc_t',-1)
axis([0,4 ,0.6, 1.8 ]);
hold on;
numgroups = size(View1_perc, 1); 
numbars = size(View1_perc, 2)+2; 
groupwidth = min(0.8, numbars/(numbars+1.5));
    for u = 1:numbars
      % Based on barweb.m by Bolu Ajiboye from MATLAB File Exchange
      x = (1:numgroups) - groupwidth/2 + (2*u-1) * groupwidth / (2*numbars);  % Aligning error bar with individual bar
      errorbar(x, TT(:, u), vari1(:,u), 'k', 'linestyle', 'none');
    end

TT = [PHQ, View2_perc];
figure; bar(TT,1);
set(gca,'XTickLabel',{'Cluster1','Cluster2','Cluster3'})
xlabel('View-Process','FontSize',16),ylabel('Relative distance to the average','FontSize',16)
legend('PRE-PHQ','POST-PHQ','act_1', 'act_2','act_3','act_4','audio_1', 'audio_2','audio_3','audio_4','conv_1','conv_2','conv_3','conv_4', 'lock_1','lock_2','lock_3','lock_4',-1)
axis([0,4 ,0,3 ]);        
hold on;
numgroups = size(View2_perc, 1); 
numbars = size(View2_perc, 2)+2; 
groupwidth = min(0.8, numbars/(numbars+1.5));
    for u = 1:numbars
      % Based on barweb.m by Bolu Ajiboye from MATLAB File Exchange
      x = (1:numgroups) - groupwidth/2 + (2*u-1) * groupwidth / (2*numbars);  % Aligning error bar with individual bar
      errorbar(x, TT(:, u), vari2(:,u), 'k', 'linestyle', 'none');
    end

% 
% wewant = zeros(nt, 1); 
% for j= 1: nt
%        if
% end
% 
% 
% Clus1 = rowClus1;
% Clus2 = ones(n, 1) - Clus1;
% Clusmat = [Clus1, Clus2];
% View1_perc = Clusmat' * view1;
% View2_perc = Clusmat' * view2;
% % 
% View1_perc(1, :) = View1_perc(1, :) / sum(Clus1); 
% View1_perc(2, :) = View1_perc(2, :) / sum(Clus2);
% 
% Clusall = double(Clus1|Clus2);
% averView1 = (Clusall' * view1) / sum(Clusall);
% averView2 = (Clusall' * view2) /  sum(Clusall);
% 
% for i = 1 : 2
%     View1_perc(i, :) = View1_perc(i, :) ./ averView1; 
%     View2_perc(i, :) = View2_perc(i, :) ./ averView2; 
% end
% figure; bar(View1_perc,1);
% set(gca,'XTickLabel',{'Cluster1','Cluster2'}) 
% xlabel('Location1','FontSize',16),ylabel('Relative distance to the average','FontSize',16)  
% legend('Height', 'Branch', 'Nod', 'Pod', 'Seed', 'Yield', 'Weight', 'Pro', 'Oil','10','11','12',-1)
% axis([0,4 ,0.6,2.2 ]);
% 
% figure; bar(View2_perc,1);
% set(gca,'XTickLabel',{'Cluster1','Cluster2','Cluster3'}) 
% xlabel('Location2','FontSize',16),ylabel('Relative distance to the average','FontSize',16)  
% legend('Height', 'Branch', 'Nod', 'Pod', 'Seed', 'Yield', 'Weight', 'Pro', 'Oil','10','11','12','13','14','15','16',-1)
% axis([0,4 ,0.5,1.4 ]);