%  function demo2_simulated_data()
%
%--------------Brief description------------------------------------------- 
% 
% This demo correponding to A.2 in [2]: Impact of incorporating unlabeled
% samples
%
%   
% In summary:
%
%   1- Based on a training set containing  spectral vectors and the
%      respective labels, learn the regressors of a multinomial 
%      logistic regression (MLR) GEM algorithm (see Algorithm 1 in [2])
%
%    2- Based on a multi-level logistic (MLL) prior, to model the 
%       contextual spatial information of the hyperspectral images, 
%       compute the MAP segmentation  via the \alpha-Expanasion 
%       Graph-cut algorithm.
%
%    3- Based on maximum entropy (or random selection), acively (or
%       randomly) select new unlabeled samples into the training set.
%
%    4- Goto 1, until some stopping rule is met
%
% -------------------------- Input parameters -----------------------------
%
% -------------------------Algorithm parameters----------------------------
%
% alpha  -  shape parameter of Gamma distribution, 
%           default setting: alpha = [1e-1 1e-3]
%
% beta   -  inverse scale parameter of Gamma distribution 
%           default setting: beta = [1e-1 1e-3]
%
% mu      - the spatial prior regularization  parameter, controlling 
%           the degree of spatial smothness. Tune this parameter to 
%           obtain a good segmentation results.
%
% --------------------hyperpsectral data parameters------------------------
%
% no_class          - number of classes 
% dim               - dimentionality of the hyperspectral data
% n                 - number of pixels in the image: n = dim*dim
% d                 - number of bands. 
% sigman            - noise standard deviation (std) of dataset
%
% ---------------------------- demo  parameters ---------------------------
%
% GC_toolbox        - this variable takes values in  {Bagon, Bioucas}
%                     default = Bagon.
%                     Graph cut toolbox to implement the alpha-Expansion
%                     algorithm. Both toolboxes are compiled for 
%                     Matlab R2009.
%                     If you face problems, compile the Bagon toolbox in 
%                     you computer by running compile_gc.m. See details in
%                     GCmex1.2
%
% hyperspectral_image - this variable takes values in 
%                       {new_image, existing_image}.   
%                       default = existing_image. 
%                      'new_image' generates a new class label image
%                       according to the MLL distribution.
%                       'existing_image' uses an existing binary image:
%                        demo1_classes.
%
% active_selection  - this variable takes values in 
%                     (maximum_entropy,random)
%                      default = maximum_entropy
%                      Active selection approach to select unlabeled
%                      samples. More details can be found in [2].
%
% train_perclass    - number of training samples per class. 
% size_train_set    - size of labeled samples
% no_u              - number o unlabeled samples selected per iteration
% tot_sim           - number of active selection iterations
% MMiter            - number of Monte Carlo runs to compute the overall
%                     accuracy
%
% --------------------------- output parameters ---------------------------
%
% ub            -  union bound
%
% xz_length     - size of training samples, including labeled and unlabeled
%                 samples
%  
% OA_class      - classification accuracies with the GEM algorithm 
%                 (Algorithm 1) in [2]
% OA_seg        - segmentation accuracies with the segmentation algorithm
%                 Algorithm 2 in [2]
%
%%                   NOTE:
%       1. The final number of labeled training samples is given by: 
%               max{size_train_set, no_class*train_perclass} 
%
%       2.  The final number of traing samples is give by:
%                    xz_length
%
%      3.  The number of unlabeled training samples for each iteration and
%           for each MC run is given by:
%                   xz_length - max{size_train_set, K*train_perclass}
%  
% ------------------------ Example1  --------------------------------------
%
% dim = 128;             
% n   = dim^2;           
% no_classes=2;          % binary problem
% d = 224;   
% train_perclass = 40;
% size_train_set = 400;   
% sigman  = 0.4;       
% mu=4;         
% no_u = 50;     
% tot_sim = 10;       % number of iterations used in the active selection
%                        procedure
% MMiter = 10;       (smaller values  ==> faster simulations and
%                     higher uncertainties in the overall accuracies (OAs))
%   
% -------------------------------------------------------------------------

% Author: Jun Li
%         Jose M. Bioucas-Dias, Sep. 2010
%
% More details in:
%
% [1] J. Li, J. Bioucas-Dias, and A. Plaza. Semi-supervised Hyperspectral 
%     classification.  in  First IEEE GRSS Workshop on Hyperspectral Image 
%     and Signal Processing-WHISPERS'2009, Grenoble, France,  2009.

%
% [2] J. Li, J. Bioucas-Dias, and A. Plaza. Semi-supervised Hyperspectral 
%    image segmentation using multinomial logistic regression with active 
%    learning.  IEEE Transactions on Geoscience and Remote Sensing, 2009 
%    (accepted).
%


% RBF with with all bands
clear all
close all
clc

%%
%--------------------------------------------------------------------------
%                          data parameter
%--------------------------------------------------------------------------

% set GC toolbox
GC_toolbox = 'Bagon';
%GC_toolbox = 'Bioucas';

% set hyperspectral image
hyperspectral_image = 'existing_image';
% hyperspectral_image = 'new_image';

% set active selection approach
active_selection = 'maximum_entropy';
%active_selection = 'random';

dim=128;        % dimentionality
no_class = 10;  % number of classes
sigman = 0.4;  % noise standard deviation (std)
n_bands = 224;  % number of bands

%--------------------------------------------------------------------------
%                   generate an MLL image
%--------------------------------------------------------------------------
if strcmp(hyperspectral_image,'new_image')
% generate a new label image
 classes = mll(dim,mu,no_class,2);
else
% use an existing label image
 load SimulationA2_class_demo1_class_A
end

image = A(:,classes(:)) + sigman*randn(n_bands,dim*dim);

% ground truth image
truth_image = classes;

train_perclass = 40;  % number of labeled samples per class
size_train_set = 400;  % number of labeled samples 
no_u = 50;      % number of unlabeled samples selected per iteration
mu=2;                   % MLL soomthness parameter 
tot_sim = 10;  % number of iterations used in the active selection procedure
MMiter  = 1;  % number of simulations
xz_length=zeros(MMiter,tot_sim);  % number of training samples

% compute the union bound
 for no_i = 1:no_class
     for no_j = 1:no_class
         dis_ij(no_i,no_j) = norm(A(:,no_i)-A(:,no_j),'fro');
     end
 end
 dis_ij = sort(dis_ij(:));
 dis_ij0 = dis_ij == 0;
 dis_ij(dis_ij0) = [];
 mindis = min(dis_ij);
 ub = (no_class-1)/2*erfc(sqrt(mindis)/(2*sigman));

trainall(:,1) = 1:(dim*dim);
trainall(:,2) = classes(:);

%%
%--------------------------------------------------------------------------
%                         Start MC runs
%--------------------------------------------------------------------------
fprintf('\nStarting Monte Carlo runs \n\n');

for iter = 1:MMiter
    
    fprintf('MC run %d \n', iter);
    % randomly select the initial training set from the ground truth
    % image
    indexes=train_test_random_new(trainall(:,2)',train_perclass,size_train_set);
    train = image(:,indexes);
    test = trainall;
    test(indexes,:) = [];
    y1 = test(:,2);
    indx1 =  1:(dim*dim);
    indx1(indexes) = [];
    indexesall = indx1;
    y=classes(indexes);
    
    % unlabeled set
    xu = [];
    
    % reindex to natural indeces
    x = [train xu];
   
    xtest = x;
    %----------------------------------------------------------------------
    %            start active section iterations
    %----------------------------------------------------------------------
    for n_sim = 1:tot_sim    
        fprintf('Active selection iteration  %d \n', n_sim);
        x = xtest;
        % space dimension
        [d,n] =size(x);
        % build |x_i-x_j| matrix 
        nx = sum(x.^2);
        [X,Y] = meshgrid(nx);
        dist=X+Y-2*x'*x;

        scale = mean(dist(:));
        clear X Y nx
        
        sigma = 0.6;
        % build design matrix (kernel) 
        K=exp(-dist/2/scale/sigma^2);
        % set first line to one 
        K = [ones(1,n); K];

        lambda=0.004;
        tic
        [w,L1] = semi_segmentation_graph(K,y);
        t1 = toc;
    
        p = splitimage(image,xtest,w,scale,sigma);
        
        [maxp,class] = max(p);
        class_map = reshape(class,dim,dim);
        % classification overall accuracy
        OA_class(iter,n_sim) =  sum(class(indx1)==y1')/length(y1);
        
        if strcmp(GC_toolbox,'Bagon')
        %------------------------------------------------------------------
        %              graphcut by Matlab wrapper; 
        %        http://www.wisdom.weizmann.ac.il/~bagon/matlab.html          
        %------------------------------------------------------------------
        
        Dc = reshape((log(p+eps))',[dim dim no_class]);%(log(p+eps)-H)
        % energy for each class0
        Sc = ones(no_class) - eye(no_class);
        % spatialy varying part

        gch = GraphCut('open', -Dc, mu*Sc);
        [gch seg] = GraphCut('expand',gch);
        gch = GraphCut('close', gch);
        
        else
        
        %% 
        %------------------------------------------------------------------
        %          graphcut by another toolbox; Jose Bioucas Dias        
        %------------------------------------------------------------------
        % energy for each class
        e=-log(p+eps);% log(p12);
        for i=1:no_class
            eLabel(:,:,i)=reshape(e(i,:),dim,dim);
        end
 
        % % % % % %
        cliques = [ 0  1;
            1  0];
        % % % % % % %
   
        [cliquesm,cliquesn] = size(cliques);
        [em,en,nrLabels] = size(eLabel);
        labels = [1:nrLabels];
        labels = labels-1;
        % Expantion Algorithm
        initLabelling = ones(dim,dim);
        seg = alfaExpansion(eLabel,labels,initLabelling,mu,cliques);
        end
        
        %  segmentation overall accuracy
        OA_seg(iter,n_sim) =   sum(seg(indx1) + 1 ==y1')/length(y1);

        %------------------------------------------------------------------
        %           active selection:  maximum entropy
        %------------------------------------------------------------------
        if strcmp(active_selection,'maximum_entropy')
        
        pindexall = p(:,indexesall);
        [maxpall,classall] = max(pindexall);
        [sortmaxpp, indexsortmaxpp] = sort(maxpall);
        xp = indexesall(indexsortmaxpp((1:no_u)));
        % eliminate the [xp] unlabeld samples from the background unlabeled
        % set, see line 8 in Algorithm 2 of reference [2]
        indexesall(indexsortmaxpp((1:no_u))) = [];
        
        else
        
        %------------------------------------------------------------------
        %              active selection: random
        %------------------------------------------------------------------

        per_index1 = randperm(length(indexesall));
        xp = indexesall(per_index1(1:no_u));
        % eliminate the [xp] unlabeld samples from the background unlabeled
        % set, see line 8 in Algorithm 2 of reference [2]
        indexesall(per_index1(1:no_u)) = [];
        end

        xu = image(:,xp);
        z = [xtest xu];
        
        [d1,n2] =size(z);
        % number of training samples
        xz_length(iter,n_sim) = n2;
    
        xtest = z;
        
    end
end

%-------------------------------------------------------------------------%
%           evaluation of the algorithm performance
%-------------------------------------------------------------------------%
%
% compute the mean OAs over MMiter MC runs
mean_classification = mean(OA_class,1);
mean_segmentation = mean(OA_seg,1);

L_total = mean(xz_length,1);

% the training set size
L   = max(size_train_set, no_class*train_perclass) ;
U   = L_total(end) - L;
fprintf('The number of labeled training samples:  L = %d \n',L );
fprintf('The number of unlabeled training samples:  U = %d \n', U);
fprintf('The final number of  training samples: L_total = L + U =  %d \n', L_total(end));
fprintf('The classification results: OA =  %2.2f (U = 0), OA =  %2.2f (U = %d) \n',...
        mean_classification(1),mean_classification(end),U);
fprintf('The segmentation results: OA =  %2.2f (U = 0), OA =  %2.2f (U = %d) \n',...
        mean_segmentation(1),mean_segmentation(end),U);


% ground truth image
figure(1), imagesc(truth_image)
title('\bf{Ground truth map}','Fontsize', 14)

% classification map
figure(2), imagesc(class_map)
st = sprintf('classification map: OA = %2.2f, (L = %d, U = %d)',...
    mean_classification(end),L, U);
title(st,'Fontsize', 12)

% segmentation map
figure(3), imagesc(seg)
st = sprintf('segmentation map: OA = %2.2f, (L = %d, U = %d)', ...
     mean_segmentation(end),L, U);
title(st,'Fontsize', 12)
