%% demo for SVM-MLRsub (Local and Global Probabilities)
%"M. Khodadadzadeh, J. Li, A. Plaza, H. Ghassemian, J. M. Bioucas-Dias and X. Li. 
%Spectral-Spatial Classification of Hyperspectral Data Using Local and Global Probabilities for Mixed Pixel Characterization. 
%IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 10, pp. 6298-6314, October 2014"

close all
clear all
clc

addpath .\LibSVM
addpath .\GraphCutMex


% load AVIRIS_IP
% indexes = train_test_random_new_i(trainall(2,:),50,[1 7 9 16]);
% train1 = trainall(:,indexes);
% 
% % the remaining used for test
% test1 = trainall;
% test1(:,indexes) = [];
% 
% ksub = 186;
% clear img;


load ROSIS_PU
indexes = train_test_random_new_i(trainall(2,:),50,[1 7 9 16]);
train1 = trainall(:,indexes);

% the original test set
test1 = GT_indexes;

ksub = 100;
clear img;

no_train = size(train1,2);

%% Probabilistic SVM Classification

disp (['SVM: ... '])
img = reshape(im',[no_lines,no_rows,size(im,1)]);

groundtruth = zeros(no_lines, no_rows);
groundtruth(trainall(1,:))=trainall(2,:);

training2D=zeros(no_lines, no_rows);
training2D(train1(1,:))=train1(2,:);

in_param.probability_estimates = 1;
%in_param.nfold = 5;  % cross validation for parameters
in_param.cost = 125;
in_param.gamma = 2^(-6);

[map_class outdata] = classify_svm(img,training2D,in_param);

% output probabilities
[pSVM,order,ordervalue] = ...
    aux_ordenar_v4(groundtruth,outdata.prob_estimates,no_lines, no_rows);
pSVM = pSVM';

[maxSVM,SVMmap] = max(pSVM);
[SVMOA,SVMk,SVMAA,SVMCA] =...
    calcError( test1(2,:)-1, SVMmap(test1(1,:))-1, 1: no_classes);
clear img training2D fixed_tr groundtruth im_test map_class

%%  Global Subspace Classifier
disp (['Global MLR: ... '])
train = im(:,train1(1,:));
[v,d]   = eig(train*train'/length(indexes));
d = diag(d);
%
dtrue = d(ksub);

g = [];
g_all= [];

% number of classes

for k_iter = 1:no_classes
    index_k = train1(2,:) == k_iter;
    train_k = train(:,index_k);
    n_k   = size(train_k,2);
    [v,d]   = eig(train_k*train_k'/n_k);
    d = diag(d);
    
    sub_sp = d<dtrue;
    sub(k_iter) = sum(sub_sp);
    tau = sub(k_iter);
    
    
    P = v(:,1:tau)*v(:,1:tau)'*train;
    gall = zeros(1,no_train);
    
    for num_k = 1:no_train
        gall(num_k) = sqrt(P(:,num_k)'*P(:,num_k));
    end
    g = [g; gall];
    
    ggall=zeros(1,no_samples);
    
    P_all = v(:,1:tau)*v(:,1:tau)'*im;
    
    for iter_all = 1:no_samples
        ggall(iter_all) = sqrt(P_all(:,iter_all)'*P_all(:,iter_all));
    end
    g_all = [g_all;ggall];
end


% regularization parameter
w = subspace_classifier(-g,train1(2,:),eps);

% compute the probablity
p_global= subspace_mlogistic(w,-g_all);

[maxp,classmap] = max(p_global);
[MLRgOA,MLRgkappa,MLRgAA,MLRgCA] =...
    calcError( test1(2,:)-1, classmap(test1(1,:))-1, 1: no_classes);
clear P_all

%% Class Combinations Map

M = 2;
prob_estimates = pSVM';

all_comb = zeros(no_samples,M);
for i=1:no_samples
    for j=1:M
        [m l]=max(prob_estimates(i,:));
        all_comb(i,j) = l;
        prob_estimates(i,l) = 0;
    end
end
all_comb = sort(all_comb')';
C = unique(all_comb, 'rows');

all_comb_prop = cell(size(C,1),M);
for i=1:size(C,1)
    all_comb_prop{i,1} = C(i,:);
    all_comb_prop{i,2} = [];
end
for i=1:no_samples
    for j=1:size(C,1)
        if all_comb(i,:) == all_comb_prop{j,1}
            all_comb_prop{j,2} = [all_comb_prop{j,2},i];
            break;
        end
    end
end

comb_map = [];
comb_map_lab = [];
for i=1:size(all_comb_prop,1)
    I = all_comb_prop{i,2};
    comb_map = [comb_map,I];
    comb_map_lab = [comb_map_lab,i*ones(1,length(I))];
end
comb_map = [comb_map;comb_map_lab];


%%  Local Subspace Classifier
disp (['Local MLR: ... '])
p_local = zeros(no_classes,no_samples);
for i=1:size(C,1)
    fprintf([num2str(i),'/',num2str(size(C,1)),' \n']);
    el = all_comb_prop{i,1}; % existing labels
    train_samples_el=[];train_el = [];
    for k=1:length(el)
        t=(train1(2,:)==el(k));
        train_el = [train_el,find(t)];
        indx=nonzeros(t.*train1(1,:));
        train_samples_el = [train_samples_el,[indx';k*ones(1,length(indx))]];
    end
    
    % second step: learn the MLRsub logistic regressors
    w = subspace_classifier(-g(all_comb_prop{i,1},train_el),train_samples_el(2,:),eps);
    
    % compute the probablity
    p1= subspace_mlogistic(w,-g_all(all_comb_prop{i,1},all_comb_prop{i,2}));
    
    p = zeros(no_classes,length(all_comb_prop{i,2}));
    for kk=1:M
        p(el(kk),:) = p1(kk,:);
    end
    
    p_local(:,all_comb_prop{i,2})= p;
    clear p p1
end

%% Decision Fusion
lambda = 0.6;
p_global_local = (lambda*p_global+(1-lambda)*p_local);
[maxp,MLRglmap] = max(p_global_local);
[MLRglOA,MLRglkappa,MLRglAA,MLRglCA] =...
    calcError( test1(2,:)-1, MLRglmap(test1(1,:))-1, 1: no_classes);

%% MRF-Based Spectral-Spatial Classification

mu = 2;

% SVM-MLL
Dc = reshape((log(pSVM+eps))',[no_lines, no_rows, no_classes]);
Sc = ones(no_classes) - eye(no_classes);
gch = GraphCut('open', -Dc, mu*Sc);
[gch SVMmllmap] = GraphCut('expand',gch);
gch = GraphCut('close', gch);

[SVMmllOA,SVMmllkappa,SVMmllAA,SVMmllCA] =...
    calcError( test1(2,:)-1, SVMmllmap(test1(1,:)), 1: no_classes);
clear Dc

% MLRsub-MLL
Dc = reshape((log(p_global+eps))',[no_lines, no_rows, no_classes]);
Sc = ones(no_classes) - eye(no_classes);
gch = GraphCut('open', -Dc, mu*Sc);
[gch MLRgmllmap] = GraphCut('expand',gch);
gch = GraphCut('close', gch);

[MLRgmllOA,MLRgmllkappa,MLRgmllAA,MLRgmllCA] =...
    calcError( test1(2,:)-1, MLRgmllmap(test1(1,:)), 1: no_classes);
clear Dc

% SVM-MLRsub-MLL
Dc = reshape((log(p_global_local+eps))',[no_lines, no_rows, no_classes]);
Sc = ones(no_classes) - eye(no_classes);
gch = GraphCut('open', -Dc, mu*Sc);
[gch MLRglmllmap] = GraphCut('expand',gch);
gch = GraphCut('close', gch);

[MLRglmllOA,MLRglmllkappa,MLRglmllAA,MLRglmllCA] = ...
    calcError( test1(2,:)-1, MLRglmllmap(test1(1,:)), 1: no_classes);
clear Dc

disp([SVMOA,MLRgOA,MLRglOA;
    SVMmllOA,MLRgmllOA,MLRglmllOA])