.. _chap-tpComposantesDiscriminantes:

############################################################
TP++ - Analyse factorielle discriminante - Étape descriptive
############################################################
 
(les TP++ sont des compléments à l’enseignement de base et ne font pas partie des séances de travaux pratiques avec encadrement)

.. only:: html

    .. container:: notebook

        .. image:: _static/zeppelin_classic_logo.png
            :class: svg-inline

        `Cahier Zeppelin <jupyter/tpComposantesPrincipales.json>`_

Références externes utiles :

  * `Documentation Spark <https://spark.apache.org/docs/latest/>`_
  * `Documentation API Spark en Scala <https://spark.apache.org/docs/latest/api/scala/index.html#org.apache.spark.package>`_
  * `Documentation Scala <http://docs.scala-lang.org>`_

**L'objectif** de cette séance de TP est de montrer une mise en œuvre avec Spark de l'étape descriptive de l'analyse factorielle discriminante (AFD) ou *Linear Discriminant Analysis* (LDA). Cette méthode d'analyse n'est pas (encore) proposée par Spark dans MLLib.

L'analyse factorielle discriminante traite des observations caractérisées par des variables quantitatives plus une variable nominale, de « classe », et cherche à trouver dans l'espace initial constitué par les variables quantitatives un sous-espace linéaire dans lequel la séparation entre classes est optimale. Une présentation rapide de l'AFD est dans `cette section de RCP216 <http://cedric.cnam.fr/vertigo/Cours/RCP216/coursReductionVolume.html#analyse-factorielle-discriminante>`_ et une présentation plus détaillée dans `cette section de RCP208 <http://cedric.cnam.fr/vertigo/Cours/ml/coursMethodesFactorielles.html#analyse-factorielle-discriminante-afd>`_.
    
L'approche employée ici pour réaliser l'AFD consiste à

- effectuer une ACP normée des observations,

- projeter les observations et les centres des classes sur les axes principaux correspondant à des valeurs propres non nulles,

- appliquer aux observations une transformation qui égalise les variances des projections sur les axes principaux (*whitening*),

- calculer les centres de gravité des classes,

- appliquer une nouvelle ACP aux seuls centres de gravité, les axes principaux résultants étant ainsi les axes discriminants de l'analyse discriminante des observations,

- projeter les observations sur les trois premiers axes discriminants.

La mise en œuvre de l'AFD est illustrée ici sur les données « textures », voir par exemple `cette présentation <http://cedric.cnam.fr/vertigo/Cours/ml/coursMethodesFactorielles.html#exemple-donnees-textures>`_. Ces données consistent en 5500 observations caractérisées par 40 variables numériques de départ et appartenant à 11 classes.

L'approche proposée ici peut être adaptée à la visualisation d'autres données vectorielles qui sont caractérisées aussi par des étiquettes de **classes** (adaptées aux modèles décisionnels de classement) ou des étiquettes de **groupes** (résultant de l'application de modèles descriptifs issus de la classification automatique, par ex. *k-means*).

Pour la lecture des données, il faut entrer les commandes suivantes dans une fenêtre terminal :

.. code-block:: bash

    %%bash
    mkdir -p tpafd/data
    cd tpafd/data
    wget http://cedric.cnam.fr/vertigo/Cours/RCP216/docs/texture.csv

Ensuite, si vous utilisez Zeppelin, il suffit d'exécuter la cellule ci-dessous. Sinon, dans une autre fenêtre terminal, vous devez vous positionner dans le répertoire ``~/tpafd``, lancer ``spark-shell`` et entrer dans ``spark-shell`` les commandes qui suivent (noter que la séquence d'opérations indiquée ci-dessous peut être améliorée) :

.. only:: jupyter

    .. code-block:: scala

        %spark.dep
        z.load("org.vegas-viz:vegas_2.11:0.3.11")
        z.load("org.vegas-viz:vegas-spark_2.11:0.3.11")


.. code-block:: scala

    // Stage 1:
    // - load the textures dataset
    // - then find the centers of each of the 11 classes
    
    import org.apache.spark.sql.types._
    
    def genTexturesFields(from: Int, to: Int) = for (i <- from until to)
                          yield StructField("val"+i.toString, DoubleType, true)
    
    val texturesFields = genTexturesFields(0, 40)
    
    val texturesSchema = StructType(texturesFields).add("label", DoubleType, true)
    
    // read the textures dataset in a DataFrame
    val texturesDF = spark.read.format("csv").option("sep"," ").schema(texturesSchema)
                          .load("tpafd/data/texture.csv").cache()
    
    // check reading (there are 500 samples for each of the 11 classes)
    texturesDF.groupBy("label").count().show()
    
    // define the aggregation operation by a Map
    val colMap = {for (i <- 0 until 40) yield ("val"+i,"avg")}.toMap
    
    // find the center of each of the 11 classes
    val cgclassesDF = texturesDF.groupBy("label").agg(colMap)
    
    // verification
    cgclassesDF.groupBy("label").count().show
    
    // Stage 2:
    // - assemble the 1D features to prepare scaling and PCA,
    // - then fit scaler, transform data, fit PCA
    // - then project both data and class centers on the PCs
    
    import org.apache.spark.ml.feature.VectorAssembler
    
    // define assembler for data
    val colsEntree = {for (i <- 0 until 40) yield "val"+i.toString}.toArray
    
    val assembleur = new VectorAssembler().setInputCols(colsEntree).setOutputCol("features")
    
    // assemble data
    val texturesDFA = assembleur.transform(texturesDF).select("label","features")
    
    // verification
    texturesDFA.groupBy("label").count().show
    
    // define assembler for class centers (1D features have different names)
    val colsEntreeCG = {for (i <- 0 until 40) yield "avg(val"+i+")".toString}.toArray
    
    val assembleurCG = new VectorAssembler().setInputCols(colsEntreeCG).setOutputCol("features")
    
    //assemble class centers
    val cgclassesDFA = assembleurCG.transform(cgclassesDF).select("label","features")
    
    import org.apache.spark.ml.feature.StandardScaler
    import org.apache.spark.ml.feature.PCA
    
    // define scaler
    val scaler = new StandardScaler().setInputCol("features")
                                     .setOutputCol("scaledFeatures")
                                     .setWithStd(true)
                                     .setWithMean(true)    
    // define PCA
    val pca = new PCA().setInputCol("scaledFeatures").setOutputCol("pcaFeatures").setK(40)
    
    // fit scaler and inspect
    val scalerModel = scaler.fit(texturesDFA)

    scalerModel.mean
    scalerModel.std
    
    // apply scaling to the data
    val scaledTexturesDFA = scalerModel.transform(texturesDFA).select("label","scaledFeatures")
    scaledTexturesDFA.printSchema()
    
    // fit PCA and inspect
    val pcaModel = pca.fit(scaledTexturesDFA)
    
    pcaModel.explainedVariance
    
    // Stage 3:
    // - whiten the class centers,
    // - then fit a PCA to the centers
    //   (the resulting PCs are the discriminant components of the data)
    
    import org.apache.spark.ml.feature.ElementwiseProduct
    import org.apache.spark.ml.linalg.Vectors
    import scala.math.sqrt
    
    // prepare vector for PCA whitening
    val wv = Vectors.dense(pcaModel.explainedVariance.toArray.map(x => 1/(sqrt(x))))
    
    // preparer whitening transformer
    val whitener = new ElementwiseProduct().setScalingVec(wv)
                                           .setInputCol("pcaFeatures")
                                           .setOutputCol("whitenedFeatures")
    
    // apply scaling, PCA projection and then whitening to class centers
    val cgclassesDFAwhitened = whitener.transform(pcaModel.transform(scalerModel
                                                          .transform(cgclassesDFA)))
    
    // define new PCA on whitened class centers
    //  the resulting PCs are the discriminant components of the initial data
    val pcaCG = new PCA().setInputCol("whitenedFeatures").setOutputCol("ldaFeatures").setK(3)
    
    // fit new PCA and inspect
    val ldaModel = pcaCG.fit(cgclassesDFAwhitened)
    ldaModel.explainedVariance
    
    // Stage 4:
    // - apply PCA1 projection, whitening and then PCA2 projection to the initial data
    // - save resulting data to prepare visualization
    
    // apply PCA1 projection, whitening and then PCA2 projection to the initial data
    val ldaTexturesDFA = ldaModel.transform(whitener.transform(pcaModel
                                 .transform(scaledTexturesDFA)))
                                 .select("label","ldaFeatures")


Nous pouvons ensuite réaliser une visualisation avec Vegas :

.. code-block:: scala

    import org.apache.spark.ml.linalg.Vector
    import org.apache.spark.sql.functions.udf
    // Extrait le premier scalaire d'un vecteur et le place dans une colonne
    val first = udf((v: Vector) => v.toArray(0))
    // Extrait le deuxième scalaire d'un vecteur et le place dans une colonne
    val second = udf((v: Vector) => v.toArray(1))

    val coordinates = ldaTexturesDFA.withColumn("x", first($"ldaFeatures")).withColumn("y", second($"ldaFeatures"))
    coordinates.printSchema()
  
    implicit val render = vegas.render.ShowHTML(s => print("%html " + s))

    import vegas._
    import vegas.data.External._
    import vegas.sparkExt._
    Vegas("AFD 2D").withDataFrame(coordinates).mark(Point).encodeX("x", Quant).encodeY("y", Quant).encodeColor("label", Nom).show

Les données appartenant à une même classe sont représentées avec une même couleur. Vu le nombre élevé de classes (11), certaines nuances sont proches entre classes différentes. Vous pouvez trouver une projection des mêmes données sur les 2 premiers axes discriminants dans la Fig. 13 de `cette section de RCP216 <http://cedric.cnam.fr/vertigo/Cours/RCP216/coursReductionVolume.html#analyse-factorielle-discriminante>`_.