<commit_before>#!/usr/bin/env python
########################################################################################################################
#
# Implementation of Principal Component Analysis using scatter matrix inspired by <NAME>
# http://sebastianraschka.com/Articles/2014_pca_step_by_step.html#sc_matrix
#
#
# Step 1: Take the whole dataset consisting of J N-dimensional flattened images [NxJ], N=nb of vox
#
# Step 2: Compute the mean image: PSI (called mean_data_vect in the code)
#
# Step 3: Compute the covariance matrix of the dataset
#
# Step 4: Compute the eigenvectors and corresponding eigenvalues (from the covariance matrix)
#
# Step 5: Sort the eigenvectors by decreasing eigenvalues and choose k in order to keep V of them such as:
#       sum(kept eigenvalues)/sum(all eigenvalues) < k (kappa in asman article)
#           This gives W, a NxV matrix with N=nb_vox and V=nb of kept eigenvectors
#
# Step 6: Transform the input image onto the new subspace, this can be done by:
#       y = W.T*(x - mean)  where:
#           x is the input N*1 flatened image
#           .T: transpose operator
#           W [NxV]
#           y is x projected in the PCA space
#
#   TODO: add datashape arg in __init__ in order to remove if split etc...
#
# ----------------------------------------------------------------------------------------------------------------------
# Copyright (c) 2014 Polytechnique Montreal <www.neuro.polymtl.ca>
# Authors: <NAME>
# Modified: 2014-12-05
#
# About the license: see the file LICENSE.TXT
########################################################################################################################
import matplotlib.pyplot as plt
import numpy as np
import sct_utils as sct
from math import sqrt
import os


class PCA:

    def __init__(self, dataset, mean_vect=None, eig_pairs=None, k=0.80):
        # STEP 1
        self.dataset = dataset  # This should be a J*N dimensional matrix of J N-dimensional flattened images
        self.N = dataset.shape[0]  # The number of rows is the dimension of flattened images
        self.J = dataset.shape[1]  # The number of columns is the number of images
        # STEP 2
        if mean_vect is not None:
            self.mean_data_vect = mean_vect
        else:
            self.mean_data_vect = self.mean()


        n = int(sqrt(self.N))
        self.mean_image = self.mean_data_vect.reshape(n,n)

        fig=plt.figure()
        plt.imshow(self.mean_image.astype(np.float))
        plt.plot()

        if eig_pairs == None:
            # STEP 3
            self.covariance_matrix = self.covariance_matrix()
            # STEP 4 eigpairs consist of a list of tuple (eigenvalue, eigenvector) already sorted by decreasing eigenvalues
            self.eig_pairs = self.sorted_eig()
        else:
             # STEP 3
            self.covariance_matrix = None
            # STEP 4 eigpairs consist of a list of tuple (eigenvalue, eigenvector) already sorted by decreasing eigenvalues
            self.eig_pairs = eig_pairs
        # STEP 5
        self.k = k
        self.W, self.kept = self.generate_W()
        print '\n\n-------> IN PCA : '
        print '\n-> W:', self.W
        print '\n-> kept:', self.kept
        # omega is a matrix of k rows and J columns, each columns correspond to a vector projection of an image from
        # the dataset
        self.omega = self.project_dataset()

    # STEP 2
    def mean(self):
        mean_im = []
        for row in self.dataset:
            m = sum(row)/self.J
            mean_im.append(m)
        mean_im = np.array([mean_im]).T
        return mean_im

    # STEP 3
    def covariance_matrix(self):
        covariance_matrix = np.zeros((self.N, self.N))
        for j in range(0, self.J):
            covariance_matrix += float(1)/self.J*(self.dataset[:, j].reshape(self.N, 1) - self.mean_data_vect)\
                .dot((self.dataset[:, j].reshape(self.N, 1) - self.mean_data_vect).T)
        return covariance_matrix

    # STEP 4
    def sorted_eig(self):
        eigenvalues, eigenvectors = np.linalg.eig(self.covariance_matrix)

        eigenvectors = eigenvectors.astype(np.float)

        # Create a list of (eigenvalue, eigenvector) tuple
        eig_pairs = [(np.abs(eigenvalues[i]), eigenvectors[:, i]) for i in range(len(eigenvalues))
                     if np.abs(eigenvalues[i]) > 0.0000001]
        # Sort the (eigenvalue, eigenvector) tuples from high to low

        def getKey(item):
            return item[0] #sorting by eigenvalues
        eig_pairs = sorted(eig_pairs, key=getKey, reverse=True)
        return eig_pairs

    # STEP 5
    def generate_W(self):
        eigenvalues_kept = []
        s = sum([eig[0] for eig in self.eig_pairs])
        print '\n####################################\n ---> sum of eigenvalues : ', s
        first = 1
        for eig in self.eig_pairs:
            if first:
                W = np.asarray(eig[1]).reshape(self.N, 1)
                eigenvalues_kept.append(eig[0])
                first = 0
            else:
                if (sum(eigenvalues_kept) + eig[0])/s <= self.k:
                    eigenvalues_kept.append(eig[0])
                    W = np.hstack((W, np.asarray(eig[1]).reshape(self.N, 1)))
                else:
                    break
        kept = len(eigenvalues_kept)
        print 'kept eigenvalues (PCA space dimension)  : ', kept
        return W, kept

    # STEP 6
    # Project the all dataset in order to get images "close" to the target & to plot the dataset
    def project_dataset(self):
        omega = []
        for column in self.dataset.T:
            omega.append(self.project_array(column).reshape(self.kept,))
        return np.asarray(omega).T

    def project(self, target):
        # check if the target is a volume or a flat image
        coord_projected_slices = []
        # loop across all the slice
        for slice in target.data:<commit_msg>NEW: msct_pca : possibility to ignore domes modes of the pca<commit_after>#!/usr/bin/env python
########################################################################################################################
#
# Implementation of Principal Component Analysis using scatter matrix inspired by <NAME>
# http://sebastianraschka.com/Articles/2014_pca_step_by_step.html#sc_matrix
#
#
# Step 1: Take the whole dataset consisting of J N-dimensional flattened images [NxJ], N=nb of vox
#
# Step 2: Compute the mean image: PSI (called mean_data_vect in the code)
#
# Step 3: Compute the covariance matrix of the dataset
#
# Step 4: Compute the eigenvectors and corresponding eigenvalues (from the covariance matrix)
#
# Step 5: Sort the eigenvectors by decreasing eigenvalues and choose k in order to keep V of them such as:
#       sum(kept eigenvalues)/sum(all eigenvalues) < k (kappa in asman article)
#           This gives W, a NxV matrix with N=nb_vox and V=nb of kept eigenvectors
#
# Step 6: Transform the input image onto the new subspace, this can be done by:
#       y = W.T*(x - mean)  where:
#           x is the input N*1 flatened image
#           .T: transpose operator
#           W [NxV]
#           y is x projected in the PCA space
#
#   TODO: add datashape arg in __init__ in order to remove if split etc...
#
# ----------------------------------------------------------------------------------------------------------------------
# Copyright (c) 2014 Polytechnique Montreal <www.neuro.polymtl.ca>
# Authors: <NAME>
# Modified: 2014-12-05
#
# About the license: see the file LICENSE.TXT
########################################################################################################################
import matplotlib.pyplot as plt
import numpy as np
import sct_utils as sct
from math import sqrt
import os


class PCA:

    def __init__(self, dataset, mean_vect=None, eig_pairs=None, k=0.80):
        # STEP 1
        self.dataset = dataset  # This should be a J*N dimensional matrix of J N-dimensional flattened images
        self.N = dataset.shape[0]  # The number of rows is the dimension of flattened images
        self.J = dataset.shape[1]  # The number of columns is the number of images
        # STEP 2
        if mean_vect is not None:
            self.mean_data_vect = mean_vect
        else:
            self.mean_data_vect = self.mean()


        n = int(sqrt(self.N))
        self.mean_image = self.mean_data_vect.reshape(n,n)

        fig=plt.figure()
        plt.imshow(self.mean_image.astype(np.float))
        plt.plot()

        if eig_pairs == None:
            # STEP 3
            self.covariance_matrix = self.covariance_matrix()
            # STEP 4 eigpairs consist of a list of tuple (eigenvalue, eigenvector) already sorted by decreasing eigenvalues
            self.eig_pairs = self.sorted_eig()
        else:
             # STEP 3
            self.covariance_matrix = None
            # STEP 4 eigpairs consist of a list of tuple (eigenvalue, eigenvector) already sorted by decreasing eigenvalues
            self.eig_pairs = eig_pairs
        # STEP 5
        self.k = k
        self.W, self.kept = self.generate_W(modes_to_ignore=0)
        print '\n\n-------> IN PCA : '
        print '\n-> W:', self.W
        print '\n-> kept:', self.kept
        # omega is a matrix of k rows and J columns, each columns correspond to a vector projection of an image from
        # the dataset
        self.omega = self.project_dataset()

    # STEP 2
    def mean(self):
        mean_im = []
        for row in self.dataset:
            m = sum(row)/self.J
            mean_im.append(m)
        mean_im = np.array([mean_im]).T
        return mean_im

    # STEP 3
    def covariance_matrix(self):
        covariance_matrix = np.zeros((self.N, self.N))
        for j in range(0, self.J):
            covariance_matrix += float(1)/self.J*(self.dataset[:, j].reshape(self.N, 1) - self.mean_data_vect)\
                .dot((self.dataset[:, j].reshape(self.N, 1) - self.mean_data_vect).T)
        return covariance_matrix

    # STEP 4
    def sorted_eig(self):
        eigenvalues, eigenvectors = np.linalg.eig(self.covariance_matrix)

        eigenvectors = eigenvectors.astype(np.float)

        # Create a list of (eigenvalue, eigenvector) tuple
        eig_pairs = [(np.abs(eigenvalues[i]), eigenvectors[:, i]) for i in range(len(eigenvalues))
                     if np.abs(eigenvalues[i]) > 0.0000001]
        # Sort the (eigenvalue, eigenvector) tuples from high to low

        def getKey(item):
            return item[0] #sorting by eigenvalues
        eig_pairs = sorted(eig_pairs, key=getKey, reverse=True)
        return eig_pairs

    # STEP 5
    def generate_W(self, modes_to_ignore=0):
        eigenvalues_kept = []
        s = sum([eig[0] for eig in self.eig_pairs])
        print '\n####################################\n ---> sum of eigenvalues : ', s
        first = 1
        start=modes_to_ignore
        for eig in self.eig_pairs[start:]:
            if first:
                W = np.asarray(eig[1]).reshape(self.N, 1)
                eigenvalues_kept.append(eig[0])
                first = 0
            else:
                if (sum(eigenvalues_kept) + eig[0])/s <= self.k:
                    eigenvalues_kept.append(eig[0])
                    W = np.hstack((W, np.asarray(eig[1]).reshape(self.N, 1)))
                else:
                    break
        kept = len(eigenvalues_kept)
        print 'kept eigenvalues (PCA space dimension)  : ', kept
        return W, kept

    # STEP 6
    # Project the all dataset in order to get images "close" to the target & to plot the dataset
    def project_dataset(self):
        omega = []
        for column in self.dataset.T:
            omega.append(self.project_array(column).reshape(self.kept,))
        return np.asarray(omega).T

    def project(self, target):
        # check if the target is a volume or a flat image
        coord_projected_slices = []
        # loop across all the slice
        for slice in target.data:
{'fraction_non_alphanumeric': 0.1553033464074105, 'fraction_numerical': 0.010681799215555371, 'mean_word_length': 2.2275787772690547, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3581995', 'n_tokens_mistral': 3738, 'n_tokens_neox': 3479, 'n_words': 1100}