source: trunk/private/script_correc.m @ 19

%SCRIPT_CORREC correction

%+
% module
% ======
%
% correction du TP
%
%
% DESCRIPTION
% ===========
%
% ++
%
% ``script_correc`` plot ans save several figures :
%
%  - ``script_correc_thetas.eps`` :
%
%   figure:: script_correc_thetas.eps
%   :scale: 50 %
%   :alt: script_correc_thetas.eps
%
%   This is the caption of the figure (a simple paragraph).
%
%   The legend consists of all elements after the caption.  In this
%   case, the legend consists of this paragraph and the following
%   table:

%
%  - ``script_correc_thetav.eps`` :
%
%  - ``script_correc_dphi.eps`` :
%
%  - ``script_correc_rho_670_.eps`` :
%
%  - ``script_correc_rho_765_.eps`` :
%
%  - ``script_correc_rho_865_.eps`` :
%
%  - ``script_correc_tau.eps`` :
%
%  - ``script_correc_diagramme_dispersion.eps`` :
%
%  - ``script_correc_histogramme_erreur.eps`` :
%
% If ``fid_log`` is not set, log informations are printed on stdout.
%
% If ``fid_log`` set previous launching ``script_correc`` to more than 1,
% log informations are saved in ``script_correc.log``.
%
% EXAMPLES
% ========
%
% To run and see log informations on screen :
%
% ::
%
% >> addpath('../toolbox')
% >> tp_mlp_aerosols_startup
% >> more off
% >> script_correc
%
% To run and save log informations in ``script_correc.log`` :
%
% ::
%
% >> addpath('../toolbox')
% >> tp_mlp_aerosols_startup
% >> fid_log = 6;
% >> more off
% >> script_correc
%
%
% SEE ALSO
% ========
%
% tp_mlp_aerosols_startup.m_
%
% .. _tp_mlp_aerosols_startup.m : tp_mlp_aerosols_startup.m.html
%
% intersect_octave.m_
%
% .. _intersect_octave.m : intersect_octave.m.html
%
% script_app.m_
%
% .. _script_app.m : script_app.m.html
%
% TODO
% ====
%
% improve description
%
% avoid overriding existing output files
%
% make it launchable from everywhere (path and IO directories)
%
% Running octave, ``more off`` is compulsary but I don't not yet why
%
% save output (stderr) to make comparisons possible
%
% If printer is set to png, one might see following warning using octave :
% ..
%
% gdImageStringFT: Could not find/open font while printing string ... with font Helvetica
%
% And no labels are visible on figure outputfiles.
%
% For now, find a portable solution (setting GDFONTPATH) is too tough, so
% printer is set to eps.
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% - fplod 2009-08-18T14:21:33Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * log information either on screen (stdout) or in ``script_correc.log``
%
% - fplod 2009-08-18T09:55:12Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * save figures to make comparisons possible
%   * replace error on suptitle using octave by a warning
%
% - fplod 2009-08-17T12:39:51Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add header
%   * add octave intersect_octave
%   * add octave suptitle_octave
%
%-



%Ajout du path
addpath ('../toolbox');


%% QUESTION 1
%NOTE : C'est la partie du TP à modifier car il y a trop de choses à coder.

if ~exist('LUT_NIR','var')
 load('../LUT/LUT_NIR.dat');
end

%indices pour chacune des longeurs d'onde
i670=find(LUT_NIR(:,6)==670);
i765=find(LUT_NIR(:,6)==765);
i865=find(LUT_NIR(:,6)==865);

%fusion 670nm avec 765nm
if (run_octave == 0)
 [parambase,I1,I2]=intersect(LUT_NIR(i670,[1:5]),LUT_NIR(i765,[1:5]),'rows');
else
 [parambase,I1,I2]=intersect_octave(LUT_NIR(i670,[1:5]),LUT_NIR(i765,[1:5]),'rows');
end
base=[parambase LUT_NIR(i670(I1),7) LUT_NIR(i765(I2),7)];

%fusion de base avec 670nm avec base
if (run_octave == 0)
 [parambase,I1,I2]=intersect(base(:,[1:5]),LUT_NIR(i865,[1:5]),'rows');
else
 [parambase,I1,I2]=intersect_octave(base(:,[1:5]),LUT_NIR(i865,[1:5]),'rows');
end

%base=[thetas,thetav,dphi,Rho_aer(670),Rho_aer(765),Rho_aer(865) tau]
base=[base(I1,[3:5 6:7]) LUT_NIR(i865(I2),7) base(I1,2)];
base_mix=mix_base(base); %%%%ALEA

clear i670 i765 i865 base parambase;

%% QUESTION 2
N=size(base_mix,1);
base_app=base_mix(1:round(0.02*N),:);
base_val=base_mix(round(0.8*N)+1:round(0.82*N),:);
base_test=base_mix(round(0.98*N)+1:end,:);

%% QUESTION 3 : rien a faire

%% QUESTION 4
moy=mean(base_app);
et=std(base_app);

base_app_n=cenred(base_app,moy,et);
base_val_n=cenred(base_val,moy,et);
base_test_n=cenred(base_test,moy,et);


stitle={'Thetas','Thetav','dphi','Rho(670)','Rho(765)','Rho(865)','Tau'};
for j=1:size(base_mix,2)
figure(j)
clf

subplot(3,2,1)
hist(base_app(:,j),20);
title('App denorm');

subplot(3,2,2)
hist(base_app_n(:,j),20);
title('App norm');

subplot(3,2,3)
hist(base_val(:,j),20);
title('Val denorm');

subplot(3,2,4)
hist(base_val_n(:,j),20);
title('Val norm');

subplot(3,2,5)
hist(base_test(:,j),20);
title('Test denorm');

subplot(3,2,6)
hist(base_test_n(:,j),20);
title('Test norm');

if (run_octave == 0)
 suptitle(stitle{j});
else
 warning('www : octave can not write supertitle');
end

% save figure(j)
printer='eps';
print_printer=['-d', printer];
% figure name is based on stitle{j}.
% transformation is lower and replace ( and ) by _
if (run_octave == 0)
 figname=regexprep(lower(stitle{j}), '(', '_');
else
 figname=regexprep(lower(stitle{j}), '\(', '_');
end
if (run_octave == 0)
 figname=regexprep(figname, ')', '_');
else
 figname=regexprep(figname, '\)', '_');
end
fullfilename=['./', mfilename, '_', figname , '.', printer];
clear figname;
clear printer;
print(print_printer,fullfilename);
clear print_printer;
clear fullfilename;

end

%QUESTION 5
script_app

% définition du fichier de sortie "Affichage"
% if not set, fid_log is set to 1 (stdout)
if ~ exist('fid_log','var')
 fid_log = 1;
end

% ouverture du fichier de sortie "Affichage" si différent de stdout
fid_log_is_stdout = isequal(fid_log,1);

if fid_log_is_stdout == 0
 fullfilename = ['./', mfilename, '.log'];
 fid_log = fopen(fullfilename,'w');
 clear fullfilename;
end

%Affichage du tableau des performances
fprintf(fid_log,'\n\t\tInit1\t\tInit2\n');
for n=1:length(archi)
 fprintf(fid_log,'%d neurones\t%3.2e\t%3.2e\n',archi(n),perf(n,1),perf(n,2));
end


%On determine l'architecture optimale
mini=min(perf(:));
[n,i]=find(perf==mini);
fprintf(fid_log,'\n Architecture optimale : %d neurones (RMS=%3.2e)\n',archi(n),perf(n,i));

%FACULTATIF : On prolonge l'apprentissage avec une base plus grosse
base_app=base_mix(1:round(0.1*N),:);
  %Apprentissage du rseau
    [nets(n,i).W1 nets(n,i).W2]=MLPfit(base_val_n(:,1:end-1),base_val_n(:,end),base_app_n(:,1:end-1), ...
                                       base_app_n(:,end),nets(n,i).W1,nets(n,i).W2,'tah','lin',[1500,1e-6],[100,1]);


%QUESTION 6

    %Restitution du reseau
  tau_nn=MLPval(base_val_n(:,1:end-1),nets(n,i).W1,nets(n,i).W2);

  %Performance du reseau
  delta=decenred(tau_nn,moy(end),et(end)) - decenred(base_val_n(:,end),moy(end),et(end));
  perf_opti=sqrt((1./length(delta))*sum(delta.^2));

  fprintf(fid_log,'\nRMS finale : %3.2e\n',perf_opti);

  % fermeture du ficher de sortie "Affichage" si différent de stdout
  if fid_log_is_stdout == 0
   fclose(fid_log);
  end
  clear fid_log_is_stdout;

  figure
  plot(decenred(base_val_n(:,end),moy(end),et(end)),decenred(tau_nn,moy(end),et(end)),'.');
  hold on
  plot([0 0.5],[0 0.5],'-k');
  title('Diagramme de dispersion')
  xlabel('TAU(vrai)')
  ylabel('TAU(estime)');

  % save figure
  printer='eps';
  print_printer=['-d', printer];
  fullfilename=['./', mfilename, '_diagramme_dispersion', '.', printer];
  clear printer;
  print(print_printer,fullfilename);
  clear print_printer;
  clear fullfilename;

  figure
  hist(delta,50);
  xlabel('TAU(estime)-TAU(vrai)');
  ylabel('Nombre');

  % save figure
  printer='eps';
  print_printer=['-d', printer];
  fullfilename=['./', mfilename, '_histogramme_erreur', '.', printer];
  clear printer;
  print(print_printer,fullfilename);
  clear print_printer;
  clear fullfilename;