%% Ce script fait appel a des fonctions/instructions incluses dans 
% le package SAISIR
% Vous pouvez telecharger le package SAISIR ‡ l'adresse
% suivante :
% http://www.chimiometrie.fr/saisirdownload.html
%
%% Representation des donnees
figure; % spectre Raman
plot(Bacillus_subtilis_II.d(1,:),'r');
xlabel('Raman shift (cm-1)'); 
ylabel('Signalintensity');
title('Raman spectrum'); 
subtitle('microorganism: Bacillus subtilis')

figure; % spectre Raman
plot(Aeromonas_bestiarum.d(1,:),'b'); title('Raman spectrum');
subtitle('microorganism: Aeromonas bestiarum')
xlabel('Raman shift (cm-1)'); ylabel('Signal intensity');

%% mise en forme de la matrice de donnees n = 86; p = 5641
mydata = appendrow1(Bacillus_subtilis_II, Aeromonas_bestiarum);
% appendrow1 est une fonction du package SAISIR telechargeable 
% sur www.chimiometrie.fr

%% Analyse en Composantes Principales avec le package SAISIR
mypca = pca(mydata); % la structure "mypca" contient l'ensemble des 
                     % resultats
figure; % Creation d'une fenetre graphique pour permettre la creation 
        % d'une figure
        % on regarde les SCORES, c'est a dire la position relative des 
        % echantillons dans l'espace
        % des composantes principales
colored_map1(mypca.score, 1,2, 3,6); % 1,2 signifie qu'on souhaite voir 
% le plan forme par PC1 et PC2
% version de l'affichage sans les commandes du package SAISIR, code Matlab
% classique plot : fonction Matlab permettant l'affichage de courbes, 
% de points de toutes valeurs numeriques
plot(mypca.score.d(1:43,1),mypca.score.d(1:43,2),'bo'); hold on;
plot(mypca.score.d(44:86,1),mypca.score.d(44:86,2),'rs');
xlabel('Composante principale n°1 (PC1)'); % on identifie l'axe horizontal
ylabel('Composante principale n°2 (PC2)'); % on identifie l'axe vertical
 
figure;
plot(mypca.eigenvec.d(:,1),'b'); % Affichage des loadings sur PC1
figure;
plot(mypca.eigenvec.d(:,2),'r'); % Affichage des loadings sur PC2

%% Version de l'ACP avec le code Matlab classique
[n,p]=size(mydata.d);       % on calcul la taille de l'ensemble de donnees
average.d=mean(mydata.d);   % on fait la moyenne par colonne
centred = mydata.d-ones(n,1)*average.d; % on centre les donnees
[U S V] = svd(centred); T = U*S;        % ACP avec l'algorithme SVD 
                                        % (Singular Value Decomposition)
figure;                                 % + calcul des SCORES par le 
                                        % produit U*S

% Affichage des SCORES
plot(T(1:43,1),T(1:43,2),'bo'); hold on; % hold on permet de superposer 
                                         % les 2 series de points
plot(T(44:86,1),T(44:86,2),'rs');
% on constate des differences avec les resultats precedents dans les
% Echelles, pourquoi ? il y a un coefficient de ponderation
plot(T(1:43,1)/sqrt(n),T(1:43,2)/sqrt(n),'bo'); hold on; % hold on permet 
                                % de superposer les 2 series de points
plot(T(44:86,1)/sqrt(n),-T(44:86,2)/sqrt(n),'rs');

