<jupyter_start><jupyter_text>Load the data from the csv, ignores lines with comments.The result is a simple python array<jupyter_code>import csv

data = []
with open('enzymes_5.csv', newline='') as csvfile:
    reader = csv.reader(csvfile, delimiter=',', quotechar='|')
    for row in reader:
        if len(row) > 0 and not row[0].startswith('#'):
            data.append([float(r) for r in row])
#data
!ls
def normalize_dataset(data):
    for column_id in range(len(data[0])):
        column_max = max([x[column_id] for x in data])
        for i, entry in enumerate([x[column_id] for x in data]):
             data[i][column_id] = entry/column_max
    return data
#data = normalize_dataset(data)
#for d in data:
#    print(d[0], d[1])
from scipy.stats import chisquare,chi2
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
from scipy.spatial import distance
import scipy as sp
import sys
from sklearn import mixture
#data = np.random.uniform(size=[50000, 10])
from sklearn.datasets.samples_generator import make_blobs

centers = [(0.2, 0.2), (0.2, 0.8), (0.5, 0.5)]
cluster_std = [0.05, 0.05, 0.05]

#data, y = make_blobs(n_samples=1000, cluster_std=cluster_std, centers=centers, n_features=3, random_state=1)

#plt.scatter(data[y == 0, 0], data[y == 0, 1], color="red", s=10, label="Cluster1")
#plt.scatter(data[y == 1, 0], data[y == 1, 1], color="blue", s=10, label="Cluster2")
#plt.scatter(data[y == 2, 0], data[y == 2, 1], color="green", s=10, label="Cluster3")
import csv

with open('synthetic.csv', mode='w') as employee_file:
    employee_writer = csv.writer(employee_file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
    for l in data:
        employee_writer.writerow(l)<jupyter_output><empty_output><jupyter_text>Class to store data about bins.Makes it easier to keep track of marked bins and merging them.Holds the index of the bin in the initial array, mapping to the dataset.Holds an array of all subsequent bins that should be merged to itself.<jupyter_code># TODO create datapoint class

class Bin:
    def __init__(self, index, interval, dim):
        self.interval = interval
        self.marked = False
        self.support = 0
        self.index = index
        self.merge_with = []
        self.assigned_points = []
        self.dimension = dim
        
    def add_point(self, point):
        self.support += 1
        self.assigned_points.append(point)
        
    def get_width(self):
        return self.interval.length
    
    def __str__(self):
        return 'Interval: {0}; Marked: {1}; Support: {2}; Index: {3}; \
            Merge-with {4}; # of assigned points {5}; Dimension: {6}'.format(
            self.interval,
            self.marked,
            self.support,
            self.index,
            self.merge_with,
            len(self.assigned_points),
            self.dimension
        )
    
class Interval:
    def __init__(self, start, end):
        self.start = start
        self.end = end
        self.length = end - start
    
    def __str__(self):
        return 'Start: {0}; End: {1}; Length {2}'.format(self.start, self.end, self.length)<jupyter_output><empty_output><jupyter_text>Calculates the interval, splits the data into bins.TODO: Works on one dimension only, should be generalised to an arbitrary number of dimensions according to the dataset<jupyter_code>def split_into_bins(data):
    bins = []
    for col_idx in range(len(data[0])):
        column_bins = []
        column_data = [x[col_idx] for x in data]
        col_min = min(column_data)
        col_max = max(column_data)

        interval_length = col_max - col_min
        # divide interval 1 + log2(n) bins
        nr_of_bins = math.floor(1 + math.log2(len(column_data)))
        
        column_bins = []
        b = interval_length / nr_of_bins
        for i in range(nr_of_bins):
            # adds a bin with a start and end interval
            bin_interval = Interval(col_min + b * i, col_min + b * (i + 1))
            
            column_bins.append(Bin(i, bin_interval, col_idx))

        # compute support for each bin
        for i, datapoint in enumerate(column_data):
            bin_index = int((datapoint - col_min) / b)
            bin_index = bin_index if bin_index < nr_of_bins else bin_index - 1
            column_bins[bin_index].add_point(data[i])

        bins.append(column_bins)
    return bins, nr_of_bins

bins, nr_of_bins = split_into_bins(data)
#plt.scatter([x[0] for x in data], [y[1] for y in data])
#plt.scatter([b.interval[0] for b in bins[0]], [0.5 for _ in range(len(bins[0]))], c='green', s=100)
#plt.scatter([b.interval[1] for b in bins[0]], [0.5 for _ in range(len(bins[0]))], c='red')
#plt.show()<jupyter_output><empty_output><jupyter_text>"On the attributes deemed non-uniform, the bin with the largest support is marked. The remaining un-marked bins are tested again using the Chi-square test for uniform distribution. If the Chi-square test indicates that the un-marked bins “look” uniform, then we stop. Otherwise, the bin with the second-largest support is marked. Then, we repeat testing the remaining un-marked bins for the uniform distribution and marking bins in decreasing order of support, until the current set of un-marked bins satisfies the Chi-square test for uniform distribution."<jupyter_code># find the bin with the highest support and mark it
def mark_highest_support(column_bins):
    max_support = 0
    max_index = 0
    for _bin in column_bins:
        if _bin.marked:
            continue
        if _bin.support > max_support:
            max_support = _bin.support
            max_index = _bin.index
    #print(max_support, max_index)
    column_bins[max_index].marked = True

# perform chisquared for the support 
def mark_bins(column_bins, alpha=0.001, stat=1):
    while (stat > alpha):
        # support list of all *unmarked* bins
        support_list = [column_bins[i].support for i in range(nr_of_bins) if not column_bins[i].marked]
        #print(support_list)
        # if there are no unmarked bins, end the process
        if len(support_list) == 0: 
            break
        (stat, p) = chisquare(support_list)
        #print('stat', stat)
        #print('p', p)
        if (stat > alpha):
            mark_highest_support(column_bins)
            

    #print(list(column_bins[i].support for i in range(nr_of_bins) if not column_bins[i].marked))
    
for column_bins in bins:
    mark_bins(column_bins)<jupyter_output>F:\ProgramData\Anaconda3\lib\site-packages\scipy\stats\stats.py:5745: RuntimeWarning: invalid value encountered in true_divide
  terms = (f_obs - f_exp)**2 / f_exp
F:\ProgramData\Anaconda3\lib\site-packages\scipy\stats\_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in greater
  return (a < x) & (x < b)
F:\ProgramData\Anaconda3\lib\site-packages\scipy\stats\_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in less
  return (a < x) & (x < b)
F:\ProgramData\Anaconda3\lib\site-packages\scipy\stats\_distn_infrastructure.py:1912: RuntimeWarning: invalid value encountered in less_equal
  cond2 = cond0 & (x <= _a)<jupyter_text>Find the adjacent bins that have the same marked status. Add all following bins with the same status to `merge_with` array in the Bin objects<jupyter_code># merge bins
def mark_merge_bins(column_bins):
    for i, _bin1 in enumerate(column_bins):
        for j, _bin2 in enumerate(column_bins[i+1:]):
            if _bin1.marked == _bin2.marked:
                _bin1.merge_with.append(_bin2.index)
            else:
                break
                
for column_bins in bins:
    mark_merge_bins(column_bins)<jupyter_output><empty_output><jupyter_text>Merges the bins according to the merge list. New bins in a new array.<jupyter_code># merge each bin in the list by extending the interval and combining support
def merge_bin(all_column_bins, column_bin):
    for bin_index in column_bin.merge_with:
            column_bin.interval.end = all_column_bins[bin_index].interval.end
            column_bin.support += all_column_bins[bin_index].support
            column_bin.assigned_points.extend(all_column_bins[bin_index].assigned_points)

# merge bins of a single column
def merge_column_bins(column_bins):
    i = 0
    new_bins = []
    while i < len(column_bins):
        # if bin has no following bins to merge with, keep it as is, and go to the next one
        if len(column_bins[i].merge_with) == 0:
            new_bins.append(column_bins[i])
            i += 1
            continue
            
        merge_bin(column_bins, column_bins[i])
        
        new_bins.append(column_bins[i])
        
        # skip all of the bins that were included in the current one
        i = max(column_bins[i].merge_with) + 1
        
    return new_bins

def merge_all_bins(bins):
    new_bins = []
    for column_bins in bins:
        new_bins.append(merge_column_bins(column_bins))
    return new_bins

new_bins = merge_all_bins(bins)
dim1 = 17
dim2 = 22

dims = [i for i, x in enumerate(new_bins) if len(x) > 3]

for i in range(0, len(dims)-1, 2):
    dim1 = dims[i]
    dim2 = dims[i + 1]
    plt.scatter([x[dim1] for x in data], [y[dim2] for y in data])
    plt.scatter([b.interval.start for b in new_bins[dim1]], [0.5 for _ in range(len(new_bins[dim1]))], c='green', s=100)
    plt.scatter([b.interval.end for b in new_bins[dim2]], [0.5 for _ in range(len(new_bins[dim2]))], c='red')

    plt.scatter([0.5 for _ in range(len(new_bins[dim1]))], [b.interval.start for b in new_bins[dim1]], c='yellow', s=100)
    plt.scatter([0.5 for _ in range(len(new_bins[dim2]))], [b.interval.end for b in new_bins[dim2]], c='purple')

    plt.show()

sum(x.support for x in new_bins[0]) # confirm all points are kept
from scipy.stats import poisson

def create_new_candidate(candidate, dim_bin, reevaluated_points):
    current_bins_list = []
    current_bins_list.extend(candidate.bins)
    current_bins_list.append(dim_bin)
    return PSignature(current_bins_list, reevaluated_points)

def generate_candidate_list(candidate_list, current_dim_bins, threshold, current_dim):
    new_candidates = []
    for candidate in candidate_list:
        for dim_bin in current_dim_bins:
            if dim_bin.marked:
                expected_sup = candidate.get_support() * dim_bin.get_width()            
                reevaluated_points = candidate.reevaluate_assigned_points(dim_bin, current_dim)
                r_support = len(reevaluated_points)
                if r_support == 0:
                    continue
                print('R support {0}, expected support {1}'.format(r_support, expected_sup))
                print('Poisson distribution:', poisson.pmf(r_support, expected_sup), r_support, expected_sup)

                if poisson.pmf(r_support, expected_sup) < threshold:
                    new_candidate = create_new_candidate(candidate, dim_bin, reevaluated_points)
                    new_candidates.append(new_candidate)
                    print("Length of new candidates after poisson", len(new_candidates))
    return new_candidates
class PSignature:
    def __init__(self, bins, assigned_points=[]):
        self.bins = bins
        self.assigned_points = assigned_points
        
    def get_support(self):
        return len(self.assigned_points)
    
    def add_bin(self, _bin):
        self.bins.append(_bin)
        self.assigned_points.append(_bin.assigned_support)
        self.reevaluate_assigned_points(current_dim)
    
    def reevaluate_assigned_points(self, _bin, current_dim):
        evaluated_points = []
        current_interval = _bin.interval
        for point in self.assigned_points:
            if point[current_dim] > current_interval.start and point[current_dim] <= current_interval.end:
                evaluated_points.append(point)
        return evaluated_points
    
    def get_means(self):
        return np.average(np.array(self.assigned_points), axis = 0)
candidate_list = []
for _bin in new_bins[0]:
    if _bin.marked:
        candidate_list.append(PSignature([_bin], assigned_points=_bin.assigned_points))

poisson_threshold = 1e-4
for dim in range(1, len(data[0])):    
    candidate_list = generate_candidate_list(candidate_list, new_bins[dim], poisson_threshold, dim)
for c in candidate_list:
    print("---")
    for b in c.bins:
        print(b.index, b.dimension, len(b.assigned_points), b.marked)
        print(b.interval)
    print(len(c.assigned_points))
for bin in new_bins[1]:
    print(bin.marked)
for candidate in candidate_list:
    for _bin in candidate.bins:
        print('start {0}, end {1}, length {2}'.format(_bin.interval.start, _bin.interval.end, _bin.interval.length))
class DataPoint:
    def __init__(self, coords):
        self.coords = coords
        self.assigned_clusters = []
        
    def __eq__(self, other):
        if len(other.coords) != len(self.coords):
            return False
        for i, x in enumerate(self.coords):
            if x - other.coords[i] > 1e-9:
                return False
        return True
for c in candidate_list:
    print(c.get_means())

inv_cov_cluster_dict = dict()

for i,can in enumerate(candidate_list):   
    cov = np.cov(np.array(can.assigned_points).T)
    inv_covmat= np.linalg.inv(cov)
    inv_cov_cluster_dict[i] = inv_covmat
#fuzzy membership matrix
#initialize matrix with datapoints in one column and found cluster (e.g. 1,2,3) in other column
#initialize clusterpoints with a 1 at the matrix intersection
matrix = np.zeros(dtype='float', shape=(len(data), len(candidate_list)))
dps = []

print(matrix.shape)

cov_dat = np.cov(np.array(data).T)
inv_covmat_dat= np.linalg.inv(cov_dat)

for i, point in enumerate(data):
    if i % 100 == 0:
        print(i)
    data_point = DataPoint(point)
    for j, candidate in enumerate(candidate_list):
        candidate_data_points = [DataPoint(p) for p in candidate.assigned_points]
        if data_point in candidate_data_points:
            matrix[i][j] = 1
            data_point.assigned_clusters.append(j)
    fraction = 1 if len(data_point.assigned_clusters) == 0 else 1 / len(data_point.assigned_clusters) 
    for r in range(len(candidate_list)):
        if matrix[i][r] == 1:
            matrix[i][r] = fraction
    #"""
    if len(data_point.assigned_clusters) == 0:
        closest = sys.maxsize
        closest_candidate_idx = 0
        for idx, c in enumerate(candidate_list):
            mh_distance = distance.mahalanobis(data_point.coords, c.get_means(), inv_cov_cluster_dict[idx])
            if mh_distance < closest:
                closest = mh_distance
                closest_candidate_idx = idx
        data_point.assigned_clusters.append(closest_candidate_idx)
        matrix[i][closest_candidate_idx] = 1
    #"""
    dps.append(data_point)
                
#compute mean of support set of cluster

#compute the shortest mahalanobis distance(scipy.spatial.distance.mahalanobis) 
# of unassigned points to cluster core and assign

# EM -> probably need to implement ourself
means_before = np.array([c.get_means() for c in candidate_list])

gmm = mixture.BayesianGaussianMixture(n_components=len(candidate_list), covariance_type='full').fit(matrix)

#gmm2 = mixture.BayesianGaussianMixture(n_components=len(candidate_list), covariance_type='full').fit(data)

result = gmm.predict(matrix)

#result2 = gmm2.predict(data)
clustered_points = list()
projected_cluster_dict = dict()
for c in range(len(candidate_list)):
    projected_cluster_dict[c] = []

for assigned_cluster, p in list(zip(result, data)):
    clustered_points.append((assigned_cluster,p))
    if assigned_cluster in projected_cluster_dict:
        projected_cluster_dict[assigned_cluster].append(p)

means_after_bgm = list()        
for pj in projected_cluster_dict.keys():
    means_after_bgm.append(np.mean(np.array(projected_cluster_dict[pj]), axis = 0))


    
amount = 0    
for pj in projected_cluster_dict.keys():
    amount += len(projected_cluster_dict[pj])

print(amount)
    
plt.scatter([x[0] for x in data], [y[1] for y in data], c=result, s=20)
plt.scatter([x[0] for x in means_after_bgm], [y[1] for y in means_after_bgm], c="green")
plt.scatter([x[0] for x in means_before], [y[1] for y in means_before], c="red")

plt.show()
print(means_after_bgm)
inv_cov_dict = dict()

for key in projected_cluster_dict.keys():   
    cov = np.cov(np.array(projected_cluster_dict[key]).T)
    inv_covmat= np.linalg.inv(cov)
    inv_cov_dict[key] = inv_covmat
    
degree_of_freedom = (len(data[0])-1) *(len(data[1])-1)
#degree_of_freedom = 10
chi_crit = chi2.ppf(0.001, df=degree_of_freedom)

noise_cluster_idx  = len(candidate_list)
print(chi_crit)
for i, c in enumerate(clustered_points):
    cluster_mean = means_after_bgm[c[0]] 
    md = distance.mahalanobis(c[1],cluster_mean,inv_cov_dict[c[0]])
    if md > chi_crit:
        clustered_points[i] = (noise_cluster_idx,c[1])
        print(md, clustered_points[i], chi_crit)
plt.scatter([x[1][0] for x in clustered_points], [y[1][1] for y in clustered_points], c = [z[0] for z in clustered_points])
plt.scatter([x[0] for x in means_after_bgm], [y[1] for y in means_after_bgm], c="green")
plt.show()
print(means_after_bgm)
print(gmm.means_)<jupyter_output>[[3.99514325e-01 6.00485675e-01]
 [9.99899443e-01 1.00557106e-04]]
StarcoderdataJupyterStructuredCleanDedup