source: trunk/src/MSGanalysiscrossvalidation.m @ 327

function [result] = MSGanalysiscrossvalidation(ij_intensity_concat, ij_cumul_concat, time_intensity_concat, time_cumul_concat, Longitudes, Latitudes, lonmin, lonmax, latmin, latmax)

%MSGANALYSISCROSSVALIDATION est une fonction permettant de faire la
%validation des relations entre les intensites issus de l'analyse des
%temperature de brillance par MSG et les precipitations d'apres EPSAT. Pour
%cela, des donnees aleatoires sont prises dans un domaine donnee et un
%diagramme de dispersion permet de trouver la formule reliant les deux
%grandeurs. Cette formule est ensuite appliquee a un autre echantillon
%aleatoire pour la valider.

%
%+
%
% =============================
% MSGanalysiscrossvalidation.m
% =============================
%
% .. function:: MSGanalysiscrossvalidation(ij_intensity_concat, ij_cumul_concat, time_intensity_concat, time_cumul_concat, Longitudes, Latitudes, lonmin, lonmax, latmin, latmax)
%
% DESCRIPTION
% ===========
%
% - Donnees d'entrees :
%
%   * ij_intensity_concat : Matrice comportant les donnees d'intensites et
%     dont les dimensions sont temps, longitudes, latitudes.
%   * ij_cumul_concat : Matrice contenant les donnees de precipitations
%     EPSAT cumulees et dont les dimensions sont temps, longitudes,
%     latitudes.
%   * time_intensity_concat : Variable temporelle des donnees d'intensites,
%     en jour julien.
%   * time_cumul_concat : Variable temporelle des donnees de
%     precipitations, en jour julien.
%   * Longitudes : Variable de longitudes des donnees de precipitations et
%     d'intensites.
%   * Latitudes : Variable de latitudes des donnees de precipitations et
%     d'intensites.
%   * lonmin : Borne inferieure en longitude des donnees selectionnees pour
%     le test et la validation.
%   * lonmax : Borne superieure en longitude des donnees selectionnees pour
%     le test et la validation.
%   * latmin : Borne inferieure en latitude des donnees selectionnees pour
%     le test et la validation.
%   * latmax : Borne superieure en longitude des donnees selectionnees pour
%     le test et la validation.
%
% Cette fonction permet de faire la validation des relations entre les 
% intensites issus de l'analyse des temperature de brillance par MSG et les 
% precipitations d'apres EPSAT. Pour cela, des donnees aleatoires sont 
% prises dans un domaine donnee et un diagramme de dispersion permet de 
% trouver la formule reliant les deux grandeurs. Cette formule est ensuite 
% appliquee a un autre echantillon aleatoire pour la valider.
%
% EXAMPLES
% ========
%
% Voir la 'demo'.
%
% SEE ALSO
% ========
%
% :func:`cumulEPSAT_concat`
% :func:`intensityMSG_concat`
% :func:`uniquerand`
% :func:`generalscatterplot`
%
% :func:`int2precip`
%
% :func:`extractedmatrixdata`
%
% TODO
% ====
%
% coding rules and describe demo
%
% fplod 20110804T121502Z aedon.locean-ipsl.upmc.fr (Darwin)
% 
% matlab 7.4 cratos 
%  Warning: Divide by zero.
%
% octave 3.4 cratos
%  Analyse du nuage 1 sur 1 vs  Analyse du nuage 1 sur n sous matlab 7.4
%  ...
%  invalid conversion from NaN to logical
%
%
% find a equivalent to matlab refline
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110804T121216Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add result
%   * octave and matlab 7.4 (zeus) do not like [argout1,~, argout3]=myfunc() so
%     replace ~ by temporary variable
%   * no call to refline when octave running
%
% - jaclod 2011-07-29
%
%   * Ajout de la documentation et de la demonstration.
%
% - jaclod 2011-07-28
%
%   * Creation.
%
%-
global application
%
result=-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extraction des donnees de la zone d'interet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[ij_intensity_concat, time_intensity_concat,  dim2extracted, dim3extracted] = extractedmatrixdata(ij_intensity_concat, time_intensity_concat, Longitudes, Latitudes, min(time_intensity_concat(:)), max(time_intensity_concat(:)), lonmin, lonmax, latmin, latmax);
clear dim2extracted
clear dim3extracted
[ij_cumul_concat, time_cumul_concat, dim2extracted, dim3extracted] = extractedmatrixdata(ij_cumul_concat, time_cumul_concat, Longitudes, Latitudes, min(time_intensity_concat(:)), max(time_intensity_concat(:)), lonmin, lonmax, latmin, latmax);
clear dim2extracted
clear dim3extracted

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Suppression de jours eventuels de donnees 'NaN'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Suppressions des donnees pour les jours ou nous avons que des 'NaN'.
% On cherche s'il y en a, quels sont les pas de temps concernes et leur
% nombre. On commence par les donnees d'intensites.
indnan = find(isnan(ij_intensity_concat(:,1,1)));
nbnan = size(indnan,1);

% S'il y en a, on va les supprimer.
if ~isempty(indnan);
    
    % On boucle pour supprimer chaque pas de temps concernes en commencant
    % par les derniers (si on souhaite par exemple supprimer les pas de
    % temps 2 et 5 et que l'on supprime le 2 d'abord, le 5e pas de temps
    % sera ensuite decale au 4e. En commencant par la fin, cela ne se
    % produit pas).
    for nansuppr = 1 : nbnan;
       
        ij_intensity_concat(indnan(nbnan - nansuppr + 1),:,:) = [];
        
        time_intensity_concat(indnan(nbnan - nansuppr + 1)) = []; 
        
    end
    
end

clear('indnan','nbnan','nansuppr');

% De meme pour les donnees EPSAT.
indnan = find(isnan(ij_cumul_concat(:,1,1)));
nbnan = size(indnan,1);

if ~isempty(indnan);
    
    for nansuppr = 1 : nbnan;
       
        ij_cumul_concat(indnan(nbnan - nansuppr + 1),:,:) = [];
        
        time_cumul_concat(indnan(nbnan - nansuppr + 1)) = []; 
        
    end
    
end

clear('indnan','nbnan','nansuppr');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recuperation des donnees existantes pour MSG et EPSAT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% On passe d'informations sur la date et l'heure en des informations
% uniquement sur la date (jour julien entier).
time_intensity_concat = fix(time_intensity_concat);
time_cumul_concat = fix(time_cumul_concat);

% On ne garde que les jours commums.
time_concat = intersect(time_intensity_concat,time_cumul_concat);

% Nombre de jours en communs.
nbjour = size(time_concat,2);

% Initialisation des matrices contenant pour MSG et EPSAT, les indices
% temporelles des donnees des jours en commun a recuperer.
indrecupMSG = zeros(1,nbjour);
indrecupEPSAT = zeros(1,nbjour);

for jour = 1 : nbjour;
   
    % On cherche quel est l'indice du jour dans les donnees temporelles
    % pour chacun des cas.
    indrecupMSG(1,jour) = find(time_intensity_concat == time_concat(jour));
    indrecupEPSAT(1,jour) = find(time_cumul_concat == time_concat(jour));
    
end

% On recupere seulement ces donnees.
ij_intensity_concat = ij_intensity_concat(indrecupMSG,:,:);
ij_cumul_concat = ij_cumul_concat(indrecupEPSAT,:,:);

clear('time_intensity_concat','time_cumul_concat','nbjour','jour','indrecupMSG','indrecupEPSAT');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fabrication des jeux de donnees
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Selection de 10 000 points uniques aleatoirements pour les donnees de
% test et de validation.
nbpoints = size(ij_intensity_concat,1)*size(ij_intensity_concat,2)*size(ij_intensity_concat,3);
[indtest] = uniquerand(10000,1,nbpoints);
[indvalid] = uniquerand(10000,1,nbpoints);

% Recuperation des donnees correspondant aux points selectionnes.
datatestMSG = ij_intensity_concat(indtest);
datatestEPSAT = ij_cumul_concat(indtest);
datavalidMSG = ij_intensity_concat(indvalid);
datavalidEPSAT = ij_cumul_concat(indvalid);

% Calcul du nombre de donnees se recouvrant.
allind = cat(2,indtest,indvalid);
recover = 20000 - size(unique(allind),2);

% Calcul du pourcentage de donnees se retrouvant dans les deux jeux.
pcrecover = recover / 10000 * 100;

% Affichaque d'un message d'information pour connaitre le pourcentage de
% donnees se recouvrant sur les deux jeux de donnees.
disp([num2str(recover) ' points se retrouvent dans le test et la validation, soit ' num2str(pcrecover) ' pour cent des donnees dans chaque diagramme.']);

clear('nbpoints','indtest','indvalid','allind','recover');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Ouverture d'une figure blanche en plein ecran.
h=figure;
set(h,'color','w','units','normalized','position',[0 0 1 1]);
    
%%%%%%%%%%%%%%%%%%
% Donnees de test
%%%%%%%%%%%%%%%%%%

% On trace le diagramme de dispersion tout en recuperant les informations
% sur les coefficients de regression lineaire et de correlation.
subplot(1,2,1)
[coeffs, r2] = generalscatterplot(datatestEPSAT, datatestMSG, 'Cumuls EPSAT (mm)', 'intensites MSG');
clear r2
set(gca,'FontSize',15,'FontAngle','italic','FontName','Times New Roman');

% Ajout de la formule de regression.
formule = (['Regression : y = ' num2str(coeffs(1)) ' * x + ' num2str(coeffs(2))]);
text(0.5,0.97,formule,...
    'HorizontalAlignment','center','verticalAlignment','middle',...
    'Units','Normalized','FontSize',13,'FontAngle','italic',...
    'FontName','Times New Roman');  

%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%
% Donnees de validation
%%%%%%%%%%%%%%%%%%%%%%%%

% Passage des donnees d'intensites en donnees de precipitation par la
% formule trouvee avec le test.
datavalidMSG = coeffs(1) * datavalidMSG + coeffs(2);

% On trace le diagramme de dispersion entre les donnees de precipitations
% d'EPSAT et celles obtenues d'apres la formule reliant les intensites aux
% precipitations.
subplot(1,2,2)
[coeffs,r2] = generalscatterplot(datavalidEPSAT, datavalidMSG, 'Cumuls EPSAT (mm)', 'cumuls par l''analyse des Tb MSG (mm)');
clear coeffs
clear r2
set(gca,'FontSize',15,'FontAngle','italic','FontName','Times New Roman');

% Ajout d'une ligne de reference (y=x).
switch application
  case {'matlab'}
     lineref = refline(1,0);
     set(lineref,'color','k','linewidth',2);
  case {'octave'}
     disp(['www : no refline available using ', application]);
  otherwise
     error('eee : unknown application running');
end

% Ajout d'un commentaire pour la legende.
comm = ('Regression (rouge), Reference y=x (noir)');
text(0.5,0.97,comm,...
    'HorizontalAlignment','center','verticalAlignment','middle',...
    'Units','Normalized','FontSize',13,'FontAngle','italic',...
    'FontName','Times New Roman');

%%%%%%%%%%%%%%%%%%%%%%%%

% Ajout d'un titre general.
Titlegen = (['Validation croisee de la relation precipitations-intensites sur une zone de ' num2str(lonmin) ' a ' num2str(lonmax) ' en longitudes et de ' num2str(latmin) ' a ' num2str(latmax) ' en latitudes (' num2str(pcrecover) ' pour cent de recouvrement)']);
text(-0.2,-0.10,Titlegen,...
    'HorizontalAlignment','center','verticalAlignment','middle',...
    'Units','Normalized','FontSize',15,'FontAngle','italic',...
    'FontName','Times New Roman');  

result=0;

end

%!demo
%! clear all
%! close all
%! varamma_startup
%! more off
%! [ij_intensity_concat, time_intensity_concat, Longitudes, Latitudes] = intensityMSG_concat('extracted', 0, 0, 0, 0, 200, 200, 2006, 08, 01, 08, 31);
%! [ij_cumul_concat, time_cumul_concat, Longitudes, Latitudes] = cumulEPSAT_concat('extracted', 0, 0, 0, 0, 200, 200, 2006, 08, 01, 08, 31);
%! result=MSGanalysiscrossvalidation(ij_intensity_concat, ij_cumul_concat, time_intensity_concat, time_cumul_concat, Longitudes, Latitudes, -20, -10, 10, 20);