% function demo_LORSAL_MLL_AS()
% Author: Jose M. Bioucas-Dias
%         Jun Li, Dec. 2009
%
% Brief description 
%  
% This demo implements the algorithm introduced in [1] (see also [2,3]). 
% on the Indian Pines data set for 16 classes
%   
% In summary:
%
%   1- Based on a training set containing  spectral vectors and the
%      respective labels, learn the regressors of a multinomial 
%      logistic regression (MLR) using LORSAL (see appendix of [1] and [3])
%
%    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 the mutual information between the MLR regressors and the 
%       class labels, acively select new  samples.
%
%    4- Goto 1, until some stopping rule is meet
%
%
%
% -- Input parameters -----------------------------------------------------
%
% ---Algorithm parameters
%
% lambda - the Laplace parameter controlling the degree of sparsity of the 
%          regressor components.
%
% beta   -  LORSAL parameter setting the augmented Lagrance weight  and
%           algorithm convergency speed (see appendix of [1]; reasonable 
%           values: beta < lambda)
% mu      - the spatial prior regularization  parameter, controlling 
%           the degree of spatial smothness. Tune this parameter to 
%           obtain a good segmentation results (reasobnable values [1, 4]).
%
% -- 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
%
%
% k                 - number of classes (K in [1]) 
% n_train           - initial nuber of training sample (L_i in [1])
% train_perclass    - number of training samples per class. 
% n_newtrian        - number of labeled samples actively selected per 
%                     iteration (u in [1])
% tot_sim           - number of active selection iterations
% MMiter            - number of Monte Carlo runs to compute the overall
%                     accuracy
%
%
% NOTE: The final number of training samples is given by 
%       max{n_train, K*train_perclass} + n_newtrian*(tot_sim-1)
%
% --- output parameters ---------------------------------------------------
%
% OA_LORSAL_AL      - classification accuracies with the LORSAL-AL 
%                     (just classification)
% OA_LORSAL_AL_MLL  - segmentation accuracies with the LORSAL-AL-MLL
%      
%  
% --- Example1 (without active selaction) ---------------------------------
% 
% K=16;             (should remain fixed in this demo)
% n_train  2073;    (20% of the ground truth)
% n_newtrian = 0;
% tot_sim = 1;
% MMiter = 10;       (smaller values  ==> faster simulations and
%                     higher uncertainties in the overall accuracies (OAs))
%  
% --- Example2 (with active selaction) ------------------------------------
%
% n_train  500;    (5% of the ground truth)
% n_newtrian = 40;
% tot_sim = 11;
% MMiter = 10;        
% -------------------------------------------------------------------------
% 
% More details in
%
% [1] J. Li, J. Bioucas-Dias, and A. Plaza, "Hyperspectral segmentation  
%  with active learning," submitted to IEEE TGRS,  2010 (link).
%
%  [2] J. Li, J. Bioucas-Dias, and A. Plaza, "Semi-supervised Hyperspectral 
%  classification using active label selection," in SPIE Europe Remote 
%  Sensing, Vol.7477, 2009.
% 
%  [3] J. Bioucas-Dias and M. Figueiredo,  "Logistic regression via variable 
%  splitting and augmented lagrangian tools"  Tech. Rep., Insituto Superior 
%  Tecnico, TULisbon, 2009
%


clear all
close all
clc

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

% number of classes
k = 16;
% number of active selection iterations
tot_sim = 2;
% number of Monte Carlo runs
MMiter = 1;
% smoothness parameter
mu = 2;
% total initial trianing sample
n_train = 500;
% number of labeled samples actively selected per iteration
n_newtrian = 40;
% number of train samples per class
train_perclass = 5;

% ground truth image
load AVIRIS_Indiana_16class
% AVIRIS Indian Pines data set
load imgreal;
im = img;
sz = size(im);
clear img;

im = ToVector(im);
im = im';
im([104:108 150:163 220],:) =[];
trainall = trainall';

% start Monte Carlo 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,n_train); 
    train1 = trainall(:,indexes);
    test1 = trainall;
    test1(:,indexes) = [];
    
    % start active section iterations
    for n_sim = 1:tot_sim
        fprintf('Active selection iteration  %d \n', n_sim);
        train = im(:,train1(1,:));
        test = im(:,test1(1,:));
        y=train1(2,:);
        % reindex to natural indeces
        x = train;
        % space dimension
        sigma = 0.8;        
        [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(:));
        % build design matrix (kernel) 
        K=exp(-dist/2/scale/sigma^2);
        % set first line to one 
        K = [ones(1,n); K];
        
        % sparcity parameter
        lambda = 0.001;
        beta = 0.0001;
        % learn the regressors
        [w,L] = LORSAL(K,y,lambda,beta);
 
        p = splitimage2(im,x,w,scale,sigma);
        
if strcmp(GC_toolbox,'Bagon')

        %% graphcut by Matlab wrapper;
        %% http://www.wisdom.weizmann.ac.il/~bagon/matlab.html  

        Dc = reshape(log(p+eps)',[sz(1) sz(2) k]);
        Sc = ones(k) - eye(k);

        gch = GraphCut('open', -Dc, mu*Sc);
        [gch seg] = GraphCut('expand',gch);
        gch = GraphCut('close', gch);
       
        
else   % GC_toolbox == Bioucas
      
        %% graphcut by another toolbox; Jose Bioucas Dias        
        % energy for each class
        e=-log(p+eps);% log(p12);
        for i=1:k
            eLabel(:,:,i)=reshape(e(i,:),sz(1),sz(2));
        end

        % cliques definition
        cliques = [ 0  1;
            1  0];
        [cliquesm,cliquesn] = size(cliques);
        [em,en,nrLabels] = size(eLabel);
        labels = [1:nrLabels];
        labels = labels-1;
        % Expantion Algorithm
        beta =2;
        initLabelling = ones(sz(1),sz(2));
        seg = alfaExpansion(eLabel,labels,initLabelling,beta,cliques);
        
end
        [maxp,class] = max(p);
        
        % LORSAL classification accuracy 
        OA_LORSAL_AL(iter,n_sim)     =  sum(class(test1(1,:))==test1(2,:))/length(test1(2,:))
        
%         figure
%         imagesc(reshape(class,sz(1),sz(2))');
        % LORSAL + MLL segmentation accuracy
        OA_LORSAL_AL_MLL(iter,n_sim) =  sum(seg(test1(1,:))+1==test1(2,:))/length(test1(2,:))
        
%         figure
%         imagesc(seg')
 
        % active selection procedure. Samples are actively selected from
        % the test set.
        pall = p(:,test1(1,:));
        [minp,classp] = max(pall);     
        [sortmminpp, indexsortminpp] = sort(minp);
        xp = indexsortminpp(1:n_newtrian);
        trainnew = test1(:,xp);
        train1 = [train1,trainnew];
        test1(:,xp) = [];        
    end
end

%-------------------------------------------------------------------------%
%  evaluation of the algorithm performance
%-------------------------------------------------------------------------%
%
% compute the mean OAs over MMiter MC runs
mean_classification = mean(OA_LORSAL_AL)
mean_segmentation = mean(OA_LORSAL_AL_MLL)


% the training set size
U   = n_newtrian*(tot_sim-1);
L_i = max(n_train, k*train_perclass);
L   = U+L_i;
fprintf('The number of initial training samples:  L_i = %d \n',L_i );
fprintf('The number of active selected training samples:  U = %d \n', U);
fprintf('The final number of  training samples: L = L_i + U =  %d \n', L);



% ground truth image
figure
true_im = zeros(sz(1),sz(2));
true_im(trainall(1,:)) = trainall(2,:);
imagesc(true_im')
title('\bf{Ground truth map}','Fontsize', 14)

% classification map
figure
class_im = zeros(sz(1),sz(2));
class_im(trainall(1,:)) = class(trainall(1,:));
imagesc(class_im')
title('\bf{Classification map obtained with LORSAL-AL algorithm}','Fontsize', 12)
xlabel('\bf{L_i= ,U= ,L= }', 'Fontsize', 14)

x = [L_i U L];
y = [L_i U L];
fid = fopen('exp.txt', 'wt');
fprintf(fid, '%6.2f %12.8f\n', y);
% fclose(fid)

% segmentation map
figure
seg_im = zeros(sz(1),sz(2));
seg_im(trainall(1,:)) = seg(trainall(1,:))+1;
imagesc(seg_im')
title('\bf{Segmentation map obtained with LORSAL-AL-MLL algorithm}','Fontsize', 12)
xlabel('\bf{L_i=}', 'Fontsize', 14)
% %-----------------------------------------------%