%matplotlib inline
from pymatbridge import Matlab
mlab = Matlab()
mlab.start()
%load_ext pymatbridge

"""
Explanation: <div align="center">
  <h2> Méthodes quantitatives en neurosciences </h2>
</div>

<div align="center">
<b><i> Cours NSC-2006, année 2015</i></b><br>
<b>Laboratoire d'analyse de données multidimensionnelle-Réponses</b><br>
     *Pierre Bellec, Yassine Ben Haj Ali*
</div>

Objectif:
Ce laboratoire a pour but de vous initier à la manipulation
d’informations multidimensionnelles avec Matlab. Nous allons pour cela
analyser des données électrophysiologiques de décharges d’activité
neuronale. Nous allons effectuer différentes opérations visant à
visualiser, résumer et modéliser ces données. Les données
sont tirée de l’expérience de Georgopoulos1982 sur l’encodage
neuronale du mouvement du bras chez un macaque avec des implants
neuronaux. L’animal commence l’expérience en fixant un curseur au centre
d’une cible, ensuite il doit rejoindre des cibles périphériques
apparaissent dans une des 8 directions arrangé en cercle. Une fois 
la cible apparue, l’animal doit attendre ( 100-1500 ms) le signal de départ avant d’aller rejoindre la cible pour
une durée de 500ms, ensuite il retourne au point de départ (au centre).
Cette séquence de mouvement est appelée essai et dans cette expérience
il y en a 47. Le but de l’expérience de Georgopoulos et collègues est de
déterminer l’orientation spatiale préférentielle du neurone en question
dans la région MI, et qu’il est possible de prédire la direction du
mouvement à partir d’enregistrements physiologiques. Leurs résultats
indiquent qu’il y a bel et bien une préférence vers les angles de
mouvement entre 90 et 180 degrés. Durant ce travail nous allons reproduire certaines des analyses de données et
la visualisation des résultats de cette expérience.
Pour réaliser ce laboratoire, il est nécessaire de récupérer les
resources suivantes sur studium:


Chap17_Data.mat: le jeu de données tiré de Georgopoulos1982.


Les scripts diagramme_dispersion.m et diagramme_dispersion_essais.m pour la Section 1.


Les scripts histogramme_essai1.m, histogramme_essais.m pour la Section 2.


Notez que le laboratoire est noté. Il faudra remettre un rapport
détaillé incluant une réponse pour l’ensemble des questions numérotées
ci dessous. Chaque réponse fera typiquement quelques lignes, incluant du
code et une figure si demandé dans l’énoncé.
Ne pas tenir compte et ne pas executer cette partie du code:
End of explanation
"""


%%matlab
load('Chap17_Data')

"""
Explanation: Section 1 : Diagramme de dispersion
Nous allons commencer par effectuer un diagramme de dispersion
(scatter plot) de l’activation d’un neurone tout au long de la durée
d’un essai. Voir le script pour suivre ces étapes.
<ol start="1">
  <h4><li>Commençons par charger les données:</li></h4>
</ol>
End of explanation
"""


%%matlab
whos

"""
Explanation: <ol start="2">
  <h4><li>La commande `whos` nous permet de déterminer quelles variables sont disponibles dans l’espace de travail, ainsi que leur type.</li></h4>
</ol>
End of explanation
"""


%%matlab
fieldnames(spike)

"""
Explanation: <font color="red">Quelles variables sont présentes? Quel est le type de la variable spike? Quelle est sa taille?</font>

<font color="blue">Les variables sont direction, go, intruction, spike et unit. La variable spike est de type structure et de taille 1x47</font>
<ol start="3">
  <h4><li>La variable spike contient les temps des potentiels d’action détectés pour
    un neurone. Chaque entrée de la structure contient les données d’un
    essai différent. Il est possible de lister les champs de la
    structure avec la commande fieldnames:</li></h4>
</ol>
End of explanation
"""


%%matlab
size(spike(1).times)

"""
Explanation: Le champ "spikes(1).times" contient les temps de décharges de potentiels d’action pour
le premier essai. La commande permet de déterminer la taille de ce
vecteur, c’est à dire le nombre de décharges détectées:
End of explanation
"""


%%matlab
size(spike(2).times) %nb de décharges pour l'essai 2

"""
Explanation: <font color="red">Combien y-a-t-il eu de décharges pour l’essai 2? pour l’essai 10?</font>
End of explanation
"""


%%matlab
size(spike(10).times) %nb de décharges pour l'essai 10

"""
Explanation: <font color="blue">Il y a 51 décharges pour l'essai 2.</font>
End of explanation
"""


 %%matlab
 t1 = spike(1).times; 
 t2 = spike(2).times; 

"""
Explanation: <font color="blue">Il y a 85 décharges pour l'essai 10.</font>
<ol start="4">
  <h4><li> La commande suivante va présenter l’ensemble des temps de décharge pour l’essai 1. </li></h4>
</ol>

&gt;&gt; spike(1).times
           -0.9893
           -0.9402
           -0.9158
        (...)

 <font color="red">Quelle est l’unité probable de ces temps? Pourquoi y-a-t-il des
valeurs négatives?</font>

<font color="blue">L'unité est en secondes et les valeurs négatives indiquent un temps de décharge qui précède le signal de départ le de l'expérience</font>
<ol start="5">
  <h4><li> On extrait les temps de décharges des deux premiers essais dans deux
    variables `t1` et `t2`: </li></h4>
</ol>
End of explanation
"""


%%matlab
figure 
hold on 

"""
Explanation: <ol start="6">
  <h4><li> On ouvre une nouvelle fenêtre, dédiée à la visualisation: </li></h4>
</ol>
End of explanation
"""


%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 1]) 
end

"""
Explanation: Attention à la deuxième instruction! Elle permet de dessiner
plusieurs objets sur une même figure, l’un à la suite de l’autre,
sans ré-initialiser la figure.

<ol start="7">
  <h4><li> Maintenant on va tracer la première ligne du diagramme. Notez que le
    nombre de décharges dans l’essai 1 est . On applique une boucle : </li></h4>
</ol>
End of explanation
"""


%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 1]) 
end
xlabel('Temp (sec)'); 
%Idem pour l’axe des y:
ylabel('Essai #')
%Enfin, on fixe les limites de l’axe des y
ylim([0 3])
% save the result
print('figure_dispersion.png','-dpng')

"""
Explanation: Notez l’utilisation de la commande line.

<ol start="8">
  <h4><li> On va maintenant ajouter un label sur l’axe des x et y . </li></h4>
</ol>
<ol start="9">
  <h4><li> Sauvegardez la figure dans un fichier png (utilisez la commande print), sous le nom figure_dispersion.png </li></h4>
</ol>
End of explanation
"""


%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 0.5]) 
end
xlabel('Temp (sec)')
%Idem pour l’axe des y:
ylabel('Essai #')
%Enfin, on fixe les limites de l’axe des y
ylim([0 3])

"""
Explanation: <font color="red"> <ol start="10">

<h4><li> Faire une nouvelle figure où chaque barre du diagramme a une hauteur
    de 0.5, plutôt que 1. Sauvegardez ce fichier sous le nom "figure_dispersion_lignes.png". </li></h4>
&lt;/ol&gt;&lt;/font&gt;
End of explanation
"""


%%matlab
% Charger les donnees
load('Chap17_Data')
% Preparer une figure
figure 
% permettre la superposition de plusieurs graphiques dans la meme figure
hold on 
% Donner un label à l'axe des x
xlabel('Temp (sec)'); 
% Donner un label à l'axe des y
ylabel('Essai #');
% Ajuster les limites de l'axe des y
ylim([0 length(spike)]);

for num_spike = 1:length(spike) %faire une boucle pour tout les essaies
    t = spike(num_spike).times; %definir la variable pour chaque essai
    for num_temps=1:length(t) %faire une boucle pour tous les points temps
        line([t(num_temps) t(num_temps)], [0+(num_spike-1) 1+(num_spike-1)]); %dessiner une line pour chaque point temps t1(i) avec longueur de 1
    end
end

"""
Explanation: <font color="red"> <ol start="11">

<h4><li> Vous allez compléter, à partir du fichier , les 4 lignes de code
    manquante à l’interieure de la boucle pour tracer tous les essais
    (47) dans une même figure. Le résultat ressemblerait à la figure
    suivante: </li></h4>
&lt;/ol&gt;&lt;/font&gt;
End of explanation
"""


%%matlab
clear         

"""
Explanation: Section 2 : Histogramme
Nous allons continuer l’exploration des données à travers un histogramme
qui résumerait la somme des activations dans un intervalle de temps
donné. Voir le script histograme_essai1.m pour reproduire les commandes suivantes:
<ol start="1">
  <h4><li>Commençons par nettoyer l’espace de travail: </li></h4>
</ol>
End of explanation
"""


%%matlab
load('Chap17_Data')

"""
Explanation: <ol start="2">
  <h4><li> Chargeons de nouveau les données: </li></h4>
</ol>
End of explanation
"""


%%matlab
centres = [-0.95:0.1:0.95];

"""
Explanation: <ol start="3">
  <h4><li>Définissons les bords et le pas des catégories de l’histogramme </li></h4>
</ol>
End of explanation
"""


%%matlab
histo = zeros(1,length(centres));

"""
Explanation: <ol start="4">
  <h4><li>Initialiser une matrice de zéros dont la longueur est égale au
    nombre d’intervalles: </li></h4>
</ol>
End of explanation
"""


%%matlab
histo = hist(spike(1).times,centres); 

"""
Explanation: <ol start="5">
  <h4><li>Récupérez le nombre de décharges par intervalle de temps dans
    l’essai numéro 1, à l’aide avec la fonction . </li></h4>
</ol>
End of explanation
"""


%%matlab
whos histo % elle est de taille 1x20

%%matlab
min(histo)
max(histo)
mean(histo)

"""
Explanation: <font color="red">Examinez le contenu de la variable histo. Quelle est sa taille ? Son minimum, maximum, sa moyenne (voir les fonctions matlab min, max, mean ).</font>
End of explanation
"""


%%matlab
bar(centres,histo);
 %Ajuster les limites de l’axe des x
xlim([-1.1 1]);
xlabel('Temps (sec)');  %Donner un label à l’axe des x
ylabel('# essai');%Donner un label à l’axe des y

"""
Explanation: <ol start="6">
  <h4><li>Dessinez l’histogramme avec la fonction `bar`. </li></h4>
</ol>
End of explanation
"""


%%matlab
%Charger les donnees
load('Chap17_Data')

% Definir les centres des intervalles pour l'histogramme
centres = [-0.95:0.1:0.95];

% Initialiser une matrice de zéros histo dont la longueur est égale au nombre d'intervalles:
histo = zeros(length(centres),1);

% Faire une boucle à travers tous les essais et recuperer le nombre de decharges par intervalle avec la fonction histc 
for jj = 1:47
    histo=histo+histc(spike(jj).times,centres);
end

% Dessiner l'histograme avec la fonction bar
bar(centres,histo);

%Ajuster les limites de l'axe des x
xlim([-1.1 1]);

%Donner un label à l’axe des x
xlabel('Temps (sec)');  

%Donner un label à l’axe des y
ylabel('# essai');

"""
Explanation: <font color="red"> <ol start="10">

<h4><li>Reprenez le code du fichier et remplissez la boucle afin de réaliser
    un histogramme pour l’ensemble des essais </li></h4>
</ol></font>
End of explanation
"""


%%matlab
clear

"""
Explanation: Section 3 : Régression
Nous allons maintenant implémenter une régression à l’aide de Matlab.
<ol start="1">
  <h4><li>Commençons par nettoyer l’espace de travail: </li></h4>
</ol>
End of explanation
"""


%%matlab
x = [ 165 165 157 170 175 165 182 178 ]';
y = [  47  56  49  60  82  52  78  90 ]';          

"""
Explanation: <ol start="2">
  <h4><li>Maintenant nous allons récuperer les données de “notes” du cours. </li></h4>
</ol>
End of explanation
"""


%%matlab
whos

"""
Explanation: <font color="red"> Quels est la taille et le contenu des vecteurs x et y? A quoi sert l'opération ' ? </font>
End of explanation
"""


%%matlab
ftheta = inline('theta(1)+theta(2)*x','theta','x');

"""
Explanation: <font color="blue"> Les variables x et y sont tout deux des vecteurs de 1 colonne et 8 lignes. L'operation ' transposer un vecteur ligne en colonne ou le contraire.</font>
<ol start="3">
  <h4><li>On définit une nouvelle fonction en ligne: </li></h4>
</ol>
End of explanation
"""


%%matlab
whos ftheta

"""
Explanation: <font color="red"> Quel est le type de la variable ftheta?</font>
End of explanation
"""


%%matlab
theta_chap = nlinfit(x, y, ftheta, [1 1] );

"""
Explanation: <font color="blue">ftheta est une fonction en ligne (inline).</font>
<ol start="4">
  <h4><li>Estimez les coefficients de régression à l’aide de la fonction : </li></h4>
</ol>
End of explanation
"""


%%matlab
theta_chap

"""
Explanation: <font color="red"> Quelles sont les valeurs de theta_chap? A quoi sert l’argument [1 1]? Essayez de reproduire l’estimation avec d’autres valeurs pour cet argument, est-ce que cela affecte le résultat?</font>

<font color="blue">Le paramètre theta_chap vaut:</font>
End of explanation
"""


%%matlab
figure
plot(x,y,'b.');
hold on
plot(x,ftheta(theta_chap,x),'r');
            

"""
Explanation: <font color="blue">L'argument [1 1] est une valeur initiale de la méthode qui cherche la valeur theta_chap. En répétant l'expérience pour plusieurs valeurs ([2 2], [30 30], [-30 -30]) on voit que la résultat theta_chap ne semble pas dépendre ici de ce paramètre. </font>
<ol start="5">
  <h4><li>Maintenant représenter le résultat de la régression. </li></h4>
</ol>
End of explanation
"""


%%matlab
figure
plot(x,y,'b.');
hold on
plot(x,ftheta(theta_chap,x),'r');
ylim([40 95])
xlabel('taille')
ylabel('poids')
print('regression_notes.png','-dpng')

"""
Explanation: <font color="red"><ol start="6">

<h4><li>Utilisez la fonction `ylim` pour changer les limites de l’axe y de 40 à 95.
    Ajouter le label `taille` sur l’axe des x avec la commande `xlabel`, et le label `poids` sur
    l’axe des y avec la commande `ylabel`. Faites une sauvegarde de cette image,
    dans un fichier `regression_notes.png`. </li></h4>
&lt;/ol&gt;&lt;/font&gt;
End of explanation
"""


%%matlab
clear
x = 0:0.1:30;
y = cos(x) + randn(1,301); 

"""
Explanation: <ol start="7">
  <h4><li>Maintenant nous allons ajuster une courbe plus complexe, un cosinus.
    On commence par simuler des données: </li></h4>
</ol>
End of explanation
"""


%%matlab
size(x)
size(y)

"""
Explanation: <font color="red">Quelle est la taille de x? La taille de y?</font>
End of explanation
"""


%%matlab
help randn

"""
Explanation: <font color="blue">Les variables x et y sont des vecteurs lignes de longueur 301.</font>

<font color="red"> A quoi sert la fonction randn (utilisez la commande help ).</font>
End of explanation
"""


%%matlab
figure
plot(x,y,'.')
print('donnees_cosinus.png','-dpng')

"""
Explanation: <font color="blue">La commande randn permet de simuler du bruit suivant une distribution normale (Gaussienne) de moyenne nulle et de variance 1.</font> 

<font color="red"> Générer un graphe de la relation entre x et y, et sauvegardez cette image dans un fichier donnees_cosinus.png.</font>
End of explanation
"""


%%matlab
ftheta = inline('theta(1)+theta(2)*cos(x-theta(3))','theta','x');

"""
Explanation: <ol start="8">
  <h4><li>On va maintenant définir une fonction de trois paramètres: </li></h4>
</ol>
End of explanation
"""


%%matlab
ftheta([0 1 1],0)

"""
Explanation: <font color="red">Quelle est la valeur de la fonction ftheta pour theta=[0 1 1] et x=0 ?</font>
End of explanation
"""


%%matlab
theta_chap = nlinfit(x, y, ftheta, [0 1 1] )    

"""
Explanation: <ol start="9">
  <h4><li>Estimez les coefficients de régression à l’aide de la fonction : </li></h4>
</ol>


<font color="red">Quelles sont les valeurs de theta_chap ?</font>
End of explanation
"""


%%matlab
figure
plot(x,y,'b');
hold on
plot(x,ftheta(theta_chap,x),'r');
            

"""
Explanation: <font color="red">A quoi sert l’argument [0 1 1] ? Essayez de reproduire l’estimation avec d’autres valeurs pour cet argument, est-ce que cela affecte le résultat?</font>

<font color="blue">L'argument [0 1 1] jour le même rôle que le paramètre [1 1] à la question 4. En utilisant une valeur différente du paramètre (par exemple [0 1 10]) on trouve une valeur différente pour theta_chap. C'est parce que la fonction ftheta est périodique, et il existe donc une infinité de valeurs d'entrée qui donnent la même sortie.</font>
<font color="red"><ol start="10">

<h4><li> Maintenant représenter le résultat de la régression. </li></h4>
</ol></font>
End of explanation
"""


