% DEMO2_JSR - Model script for the launch of two methods
%
% Example on setting matrices in a cell array, setting advanced
% options, launching two methods one after another and
% retrieving and plotting information on a same plot
% for comparison.
%
% Explanations are in the comments of the script.
%
% type
%
% >> edit demo2_JSR
%
% to access this script without launching it.
%
% This demo uses jsr_norm_conitope which requires SeDuMi. One can verify
% that it is installed by watching the result of help sedumi for instance.
%
%
clc
disp(' ')
disp(' ')
disp('--------------------- DEMO2 for JSR Toolbox ---------------------')
disp(' ')
disp(' ')
disp('Explanations are in the comments of the script.')
disp(' ')
disp('Type')
disp(' ')
disp('>> edit demo2_JSR')
disp(' ')
disp('To read the file and comments in the editor.')
disp(' ')
disp('The execution of this example takes around 5min on a normal computer.')
disp('jsr_norm_conitope will successfully stop at iteration 13.')
disp(' ')
disp('...this gives the time to understand the script during computations.')
disp(' ')
disp(' ')
disp('This Toolbox is a beta version, please report any bug, comment or ')
disp('suggestion to jsr.louvain@gmail.com')
disp(' ')
disp(' ')
disp('Press any key to open it in your editor and start ')
disp('the computations...')
pause
s = which('sedumi.m');
if isempty(s)
error(['sedumi could not be found. Please install SeDuMi properly in order to run'...
'this demo (http://sedumi.ie.lehigh.edu/)']);
end
edit demo2_JSR
% HERE STARTS THE SCRIPT AND ITS EXPLANATION
%
% First, one takes a set of matrices into the workspace, for this demo one
% has A = rand(10,10,3) saved in demo2_JSR_rndmatrix.mat. One loads it :
load demo2_JSR_rndmatrix.mat
% In the same way one could have generated the tensor with a function,
% for instance :
%
% A = genMat;
% If one has a few matrices A1, A2, A3, A4. One can set them into
% a cell array this way :
%
% M{1} = A1;
% M{2} = A2;
% M{3} = A3;
% M{4} = A4;
%
% Or equivalently :
%
% M = [{A1}, {A2}, {A3}, {A4}];
%
% If A is a tensor (3D array) [nxnxm double], one can use tens2cell to make
% it a cell array. This is the case here, hence :
M = tens2cell(A);
% One will compare the evolution of the bounds on M with
% jsr_norm_conitope and jsr_prod_pruningAlgorithm (note that the latter
% algorithm requires the matrices to have non-negative entries).
%
% And also one will look at the evolution of the size of the population
% as both methods keep only a certain number of points.
%
% The set M or the matrices themselves might be large. So one asks the
% algorithms to save the results at each iteration.
% This allows to have some results even if one
% ends up stopping the algorithms before termination.
%
% For this purpose one can use the general (i.e. working for almost any method)
% saveinIt option.
%
% This will save the potential output at each
% iteration in a user-specified .mat file :
opts_conitope = jsrsettings('saveinIt','myfile_conitope_inIteration');
opts_pruningAlg = jsrsettings('saveinIt','myfile_pruningAlg_inIteration');
% Feeding these option structures will tell the algorithms to save at each
% iteration in myfile_conitope_inIteration.mat
% and myfile_pruningAlg_inIteration.mat respectively.
% One also asks the algorithms to save in a .mat file in the end.
% One can do this by using jsrsettings on the options again and enabling
% the general saveEnd option :
opts_conitope = jsrsettings(opts_conitope,'saveEnd','myfile_conitope_end');
opts_pruningAlg = jsrsettings(opts_pruningAlg,'saveEnd','myfile_pruningAlg_end');
% These files should have appeared in the current directory at the
% end of the execution.
% One can also tell the algorithms to generate plots themselves in the end.
% They can only generate separate plots in different figures. Nevertheless
% one will retrieve the output structures info that contains the data. From
% these one will make another figure with the two bounds.
%
% Finally, one wants the version of the "conitope" method in which the
% algorithm does not reinitialize the unit ball.
opts_conitope = jsrsettings(opts_conitope,'conitope.plotBounds',1,'conitope.plotpopHist',1,'conitope.reinitBall',0);
opts_pruningAlg = jsrsettings(opts_pruningAlg,'pruning.plotBounds',1,'pruning.plotpopHist',1);
% All the names of these option fields can be found in the help of each
% method and in help jsrsettings for the general options.
% Now one launches the methods one after another :
[bounds1, prodOpt, prodV, info_conitope] = jsr_norm_conitope(M,opts_conitope);
[bounds2, P, info_pruning] = jsr_prod_pruningAlgorithm(M,opts_pruningAlg);
% When those are done there should be four plots.
% The data is also contained in the respective info structures.
% In order to compare, on the same plot, the evolution of the bounds, one
% can do for instance :
% Retrieve information from the info structures
niter_con = info_conitope.niter;
Lb_con = info_conitope.allLb;
Ub_con = info_conitope.allUb;
niter_prng = info_pruning.niter;
Lb_prng = info_pruning.allLb;
Ub_prng = info_pruning.allUb;
% Generate the plot with title and legend, this should be the fifth figure:
figure
h1 = plot(0:niter_con, Lb_con, '-+g',1:niter_con, Ub_con,'*-r'); hold on
h2 = plot(1:niter_prng, Lb_prng,'-oc',1:niter_prng,Ub_prng,'-.b');
title('Evolution of the bounds found at each depth')
xlabel('depth')
legend([h1;h2],'Lb conitope','Ub conitope','Lb pruning','Ub pruning')
% Let us have a look at the optimal products given by each method.
%
% conitope gives the index of an optimal product. This means that its
% averaged spectral radius attains the lower bound. One can rebuild
% this product with buildProduct :
msg(-1,1,'\n Index of optimal product found by conitope : %s \n', num2str(prodOpt'))
msg(-1,1,'Corresponding matrix : opt1 = M{%d}*M{%d} = \n',prodOpt(2),prodOpt(1))
opt1 = buildProduct(M,prodOpt);
disp(opt1)
% Note the convention used for the order of product indices, it is further
% described in buildProducts
% pruning gives, in output P (a cell array), the indexes of two products
% attaining the lower and upper bound respectively.
msg(-1,1,'\n Index of product attaining lower bound found by pruningAlgorithm : %s \n', num2str(P{1}))
msg(-1,1,'Corresponding matrix : opt2 = M{%d}*M{%d} = \n',P{1}(2),P{1}(1))
opt2 = buildProduct(M,P{1});
disp(opt2)
% One can check that, as expected, the averaged spectral radii are the
% same :
%
msg(-1,1,'rho(opt1)^(1/2) = %.14g \n',rho(opt1)^(1/2))
msg(-1,1,'rho(opt2)^(1/2) = %.14g \n',rho(opt2)^(1/2))
% One can also look at the evolution of the bounds found with respect to
% computation time. The structures info contain the time (in sec) at the
% end of each iteration since starting of the algorithm.
%
% Let us retrieve this information from the structures and plot it
% on one graph.
time1 = (info_conitope.timeIter)/60; hold on
time2 = (info_pruning.timeIter)/60;
figure
hold on
h11 = plot([0 time1],info_conitope.allLb,'-+g',time1,info_conitope.allUb,'-*r','MarkerSize',10);
h22 = plot(time2,info_pruning.allLb,'-oc',time2,info_pruning.allUb,'-.b','MarkerSize',10);
xlabel('time from start in min')
title('Evolution of the bounds w.r.t. to time')
legend([h11;h22],'Lb conitope','Ub conitope','Lb pruning','Ub pruning')
% Note that conitope takes an iteration 0, this is why one has to plot for
% times [0 time1].
%
% One observes that the pruning algorithm for this set of matrices was
% much faster than conitope but also gave less tight bounds.
