ema

Complexity

Total Complexity

39

Size/Duplication

Total Lines	152
Duplicated Lines	0 %

Importance

Changes

0

14 Methods

Rating	Name	Duplication	Size	Complexity
A	get_centers()	0	2	1
A	__update_mean()	0	7	2
A	process()	0	14	4
B	__init__()	0	17	5
B	get_clusters()	0	14	5
A	__get_random_means()	0	12	3
A	__expectation_step()	0	8	4
A	__maximization_step()	0	8	2
A	__probabilities()	0	7	2
A	get_covariances()	0	2	1
A	__log_likelihood()	0	11	3
A	__update_covariance()	0	8	2
A	__get_random_covariances()	0	8	2
A	__get_stop_flag()	0	7	3

"""!


@brief Cluster analysis algorithm: Expectation-Maximization Algorithm (EMA).

@details Implementation based on article:

         - 


@authors Andrei Novikov (pyclustering@yandex.ru)

@date 2014-2017

@copyright GNU Public License


@cond GNU_PUBLIC_LICENSE

    PyClustering is free software: you can redistribute it and/or modify

    it under the terms of the GNU General Public License as published by

    the Free Software Foundation, either version 3 of the License, or

    (at your option) any later version.

    

    PyClustering is distributed in the hope that it will be useful,

    but WITHOUT ANY WARRANTY; without even the implied warranty of

    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

    GNU General Public License for more details.

    

    You should have received a copy of the GNU General Public License

    along with this program.  If not, see <http://www.gnu.org/licenses/>.

@endcond


"""



import numpy;

The import numpy could not be resolved.

from pyclustering.utils import pi;


import matplotlib.pyplot as plt;

The import matplotlib.pyplot could not be resolved.

Unused matplotlib.pyplot imported as plt

from _operator import index

Unused index imported from _operator

def gaussian(data, mean = None, covariance = None):

    dimension = len(data[0]);

 

    if (mean is None):

        mean = numpy.mean(data);

     

    if (covariance is None):

        covariance = numpy.cov(data, rowvar = False);

     

    inv_variance = numpy.linalg.inv(covariance);

    right_const = 1.0 / ( (pi * 2.0) ** (dimension / 2.0) * numpy.linalg.norm(covariance) ** 0.5 );

     

    result = [];

     

    for point in data:

        mean_delta = point - mean;

        point_gaussian = right_const * numpy.exp( -0.5 * mean_delta.dot(inv_variance).dot(numpy.transpose(mean_delta)) );

        result.append(point_gaussian);

     

    return result;



class ema:

    def __init__(self, data, amount_clusters, means = None, variances = None):

        self.__data = numpy.array(data);

        self.__amount_clusters = amount_clusters;

        

        self.__means = means;

        if (means is None):

            self.__means = self.__get_random_means(data, amount_clusters);


        self.__variances = variances;

        if (variances is None):

            self.__variances = self.__get_random_covariances(data, amount_clusters);

        

        self.__rc = [ [0.0] * len(self.__data) for _ in range(amount_clusters) ];

        self.__pic = [1.0] * amount_clusters;

        self.__clusters = [];

        self.__gaussians = [ [] for _ in range(amount_clusters) ];

        self.__stop = False;



    def process(self):

        self.__clusters = None;

        

        previous_likelihood = -10000500;

        current_likelihood = -10000000;

        

        while((self.__stop is False) and (abs(numpy.min(previous_likelihood) - numpy.min(current_likelihood)) > 0.00001) and (current_likelihood < 0.0)):

            self.__expectation_step();

            self.__maximization_step();

            

            previous_likelihood = current_likelihood;

            current_likelihood = self.__log_likelihood();

            self.__stop = self.__get_stop_flag();

            print(previous_likelihood, current_likelihood);



    def get_clusters(self):

        if (self.__clusters is not None):

            return self.__clusters;

        

        self.__clusters= [];

        for index_cluster in range(self.__amount_clusters):

            cluster = [];

            for index_point in range(len(self.__data)):

                if (self.__rc[index_cluster][index_point] >= 0.5):

                    cluster.append(index_point);

            

            self.__clusters.append(cluster);

        

        return self.__clusters;



    def get_centers(self):

        return self.__means;



    def get_covariances(self):

        return self.__variances;



    def __log_likelihood(self):

        likelihood = 0.0;

        

        for index_point in range(len(self.__data)):

            particle = 0.0;

            for index_cluster in range(self.__amount_clusters):

                particle += self.__pic[index_cluster] * self.__gaussians[index_cluster][index_point];

            

            likelihood += numpy.log(particle);

        

        return likelihood;



    def __probabilities(self, index_cluster, index_point):

        divider = 0.0;

        for i in range(self.__amount_clusters):

            divider += self.__pic[i] * self.__gaussians[i][index_point];

        

        rc = self.__pic[index_cluster] * self.__gaussians[index_cluster][index_point] / divider;

        return rc;



    def __expectation_step(self):

        for index in range(self.__amount_clusters):

index is re-defining a name which is already available in the outer-scope (previously defined on line 34).

            self.__gaussians[index] = gaussian(self.__data, self.__means[index], self.__variances[index]);

        

        for index_cluster in range(self.__amount_clusters):

            for index_point in range(len(self.__data)):

                self.__rc[index_cluster][index_point] = self.__probabilities(index_cluster, index_point);

            print(self.__rc[index_cluster]);



    def __maximization_step(self):

        for index_cluster in range(self.__amount_clusters):

            mc = numpy.sum(self.__rc[index_cluster]);

            

            self.__pic[index_cluster] = mc / len(self.__data);

            self.__means[index_cluster] = self.__update_mean(index_cluster, mc);

            

            self.__variances[index_cluster] = self.__update_covariance(index_cluster, mc);



    def __get_stop_flag(self):

        for covariance in self.__variances:

            print(covariance[0])

            if (min(covariance[0]) == 0):

                return True;

        

        return False;



    def __update_covariance(self, index_cluster, mc):

        covariance = 0.0;

        for index_point in range(len(self.__data)):

            deviation = numpy.array( [ self.__data[index_point] - self.__means[index_cluster] ]);

            covariance += self.__rc[index_cluster][index_point] * deviation.T.dot(deviation);

        

        covariance = covariance / mc;

        return covariance;



    def __update_mean(self, index_cluster, mc):

        mean = 0.0;

        for index_point in range(len(self.__data)):

            mean += self.__rc[index_cluster][index_point] * self.__data[index_point];

        

        mean = mean / mc;

        return mean;



    def __get_random_covariances(self, data, amount):

        covariances = [];

        data_covariance = numpy.cov(data, rowvar = False);

        for _ in range(amount):

            random_appendix = numpy.min(data_covariance) * 0.2 * numpy.random.random();

            covariances.append(data_covariance + random_appendix);

         

        return covariances;



    def __get_random_means(self, data, amount):

        means = [];

        mean_indexes = [];

        for _ in range(amount):

            random_index = numpy.random.randint(0, len(data));

            while(random_index in mean_indexes):

                mean_indexes.append(random_index);

                random_index = numpy.random.randint(0, len(data));

            

            means.append(numpy.array(data[random_index]));

        

        return means;




# from pyclustering.samples.definitions import SIMPLE_SAMPLES, FCPS_SAMPLES;

# from pyclustering.utils import read_sample;

#   

# # sample = read_sample(SIMPLE_SAMPLES.SAMPLE_SIMPLE9);

# # ema_instance = ema(sample, 2);

# 

# sample = read_sample(SIMPLE_SAMPLES.SAMPLE_SIMPLE2);

# ema_instance = ema(sample, 3);

# 

# ema_instance.process();

# clusters = ema_instance.get_clusters();

#   

# print(clusters);