text: from cvxpy import *
  import numpy as np 
  import time, collections, os, errno, sys
  import matplotlib
  matplotlib.use('Agg')
  import matplotlib.pyplot as plt
  from Visualization_function import visualize
  from solveCrossTime import *
  from scipy import stats
  
  from sklearn import mixture
  from sklearn import covariance
  import sklearn, random
  from sklearn.cluster import KMeans
  
  import pandas as pd
  pd.set_option('display.max_columns', 500)
  np.set_printoptions(formatter={'float': lambda x: "{0:0.4f}".format(x)})
  np.random.seed(102)
  #####################PARAMETERS TO PLAY WITH 
  window_size = 10
  maxIters = 100 ##number of Iterations of the smoothening + clustering algo
  beta = 5 ## Beta parameter
  lambda_parameter = 11e-2 ## Lambda regularization parameter
  number_of_clusters = 2
  threshold = 2e-5##Threshold for plots. Not used in TICC algorithm.
  write_out_file = False ##Only if True are any files outputted
  seg_len = 300##segment-length : used in confusion matrix computation
  
  ##INPUT file location
  input_file = "Synthetic Data Matrix rand_seed =[0,1] generated2.csv"
  
  ##Folder name to store all the OUPUTS
  prefix_string = "data_lambda=" + str(lambda_parameter)+"beta = "+str(beta) + "clusters=" +str(number_of_clusters)+"/"
  
  
  
  
  ########################################################
  
  ##parameters that are automatically set based upoon above
  num_blocks = window_size + 1
  switch_penalty = beta## smoothness penalty
  lam_sparse = lambda_parameter##sparsity parameter
  maxClusters = number_of_clusters+1## Number of clusters + 1
  write_out_file = False ##Only if True are any files outputted
  num_stacked = num_blocks - 1
  ##colors used in hexadecimal format
  hexadecimal_color_list = ["cc0000","0000ff","003300","33ff00","00ffcc","ffff00","ff9900","ff00ff","cccc66","666666","ffccff","660000","00ff00","ffffff","3399ff","006666","330000","ff0000","cc99ff","b0800f","3bd9eb","ef3e1b"]
  ##The basic folder to be created
  str_NULL = prefix_string
  
  
  print("lam_sparse", lam_sparse)
  print("switch_penalty", switch_penalty)
  print("num_cluster", maxClusters - 1)
  print("num stacked", num_stacked)
  
  Data = np.loadtxt(input_file, delimiter= ",")
  
  print("completed getting the data")
  Data_pre = Data
  UNNORMALIZED_Data = Data*1000
  (m,n) = Data.shape
  len_D_total = m
  size_blocks = n
  
  ##Add an optimization function
  # def optimize():
  
  def upper2Full(a, eps = 0):
      ind = (a<eps)&(a>-eps)
      a[ind] = 0
      n = int((-1  + np.sqrt(1+ 8*a.shape[0]))/2)  
      A = np.zeros([n,n])
      A[np.triu_indices(n)] = a 
      temp = A.diagonal()
      A = np.asarray((A + A.T) - np.diag(temp))             
      return A   
  
  
  def updateClusters(LLE_node_vals,switch_penalty = 1):
  	"""
  	Uses the Viterbi path dynamic programming algorithm
  	to compute the optimal cluster assigments
  
  	Takes in LLE_node_vals matrix and computes the path that minimizes
  	the total cost over the path
  	Note the LLE's are negative of the true LLE's actually!!!!!
  
  	Note: switch penalty > 0
  	"""
  	(T,num_clusters) = LLE_node_vals.shape
  	future_cost_vals = np.zeros(LLE_node_vals.shape)
  	##compute future costs
  	for i in range(T-2,-1,-1):
  		j = i+1
  		indicator = np.zeros(num_clusters)
  		future_costs = future_cost_vals[j,:]
  		lle_vals = LLE_node_vals[j,:]
  		for cluster in range(num_clusters):
  			total_vals = future_costs + lle_vals + switch_penalty
  			total_vals[cluster] -= switch_penalty
  			future_cost_vals[i,cluster] = np.min(total_vals)
  	##compute the best path
  	path = np.zeros(T)
  
  	##the first location
  	curr_location = np.argmin(future_cost_vals[0,:] + LLE_node_vals[0,:])
  	path[0] = curr_location
  
  	##compute the path
  	for i in range(T-1):
  		j = i+1
  		future_costs = future_cost_vals[j,:]
  		lle_vals = LLE_node_vals[j,:]
  		total_vals = future_costs + lle_vals + switch_penalty
  		total_vals[int(path[i])] -= switch_penalty
  
  		path[i+1] = np.argmin(total_vals)
  
  	##return the computed path
  	return path
  
  def find_matching(confusion_matrix):
  	"""
  	returns the perfect matching
  	"""
  	_,n = confusion_matrix.shape
  	path = []
  	for i in range(n):
  		max_val = -1e10
  		max_ind = -1
  		for j in range(n):
  			if j in path:
  				pass
  			else:
  				temp = confusion_matrix[i,j]
  				if temp > max_val:
  					max_val = temp
  					max_ind = j
  		path.append(max_ind)
  	return path
  
  def computeF1Score_delete(num_cluster,matching_algo,actual_clusters,threshold_algo,save_matrix = False):
  	"""
  	computes the F1 scores and returns a list of values
  	"""
  	F1_score = np.zeros(num_cluster)
  	for cluster in range(num_cluster):
  		matched_cluster = matching_algo[cluster]
  		true_matrix = actual_clusters[cluster]
  		estimated_matrix = threshold_algo[matched_cluster]
  		if save_matrix: np.savetxt("estimated_matrix_cluster=" + str(cluster)+".csv",estimated_matrix,delimiter = ",", fmt = "%1.4f")
  		TP = 0
  		TN = 0
  		FP = 0
  		FN = 0
  		for i in range(num_stacked*n):
  			for j in range(num_stacked*n):
  				if estimated_matrix[i,j] == 1 and true_matrix[i,j] != 0:
  					TP += 1.0
  				elif estimated_matrix[i,j] == 0 and true_matrix[i,j] == 0:
  					TN += 1.0
  				elif estimated_matrix[i,j] == 1 and true_matrix[i,j] == 0:
  					FP += 1.0
  				else:
  					FN += 1.0
  		precision = (TP)/(TP + FP)
  		recall = TP/(TP + FN)
  		f1 = (2*precision*recall)/(precision + recall)
  		F1_score[cluster] = f1
  	return F1_score
  
  def compute_confusion_matrix(num_clusters,clustered_points_algo, sorted_indices_algo):
  	"""
  	computes a confusion matrix and returns it
  	"""
  	seg_len = 200
  	true_confusion_matrix = np.zeros([num_clusters,num_clusters])
  	for point in range(len(clustered_points_algo)):
  		cluster = int(clustered_points_algo[point])
  
  
  		##CASE G: ABBACCCA
  		# num = (int(sorted_indices_algo[point]/seg_len) )
  		# if num in [0,3,7]:
  		# 	true_confusion_matrix[0,cluster] += 1
  		# elif num in[1,2]:
  		# 	true_confusion_matrix[1,cluster] += 1
  		# else:
  		# 	true_confusion_matrix[2,cluster] += 1
  
  		##CASE F: ABCBA
  		# num = (int(sorted_indices_algo[point]/seg_len))
  		# num = min(num, 4-num)
  		# true_confusion_matrix[num,cluster] += 1
  
  		#CASE E : ABCABC
  		num = (int(sorted_indices_algo[point]/seg_len) %num_clusters)
  		true_confusion_matrix[num,cluster] += 1
  
  		##CASE D : ABABABAB
  		# num = (int(sorted_indices_algo[point]/seg_len) %2)
  		# true_confusion_matrix[num,cluster] += 1
  
  		##CASE C: 
  		# num = (sorted_indices_algo[point]/seg_len)
  		# if num < 15:
  		# 	true_confusion_matrix[0,cluster] += 1
  		# elif num < 20:
  		# 	true_confusion_matrix[1,cluster] += 1
  		# else:
  		# 	true_confusion_matrix[0,cluster] += 1
  
  		##CASE B : 
  		# if num > 4:
  		# 	num = 9 - num
  		# true_confusion_matrix[num,cluster] += 1
  
  		##CASE A : ABA
  		# if sorted_indices_algo[point] < seg_len:
  		# 	true_confusion_matrix[0,cluster] += 1
  
  		# elif sorted_indices_algo[point] <3*seg_len:
  		# 	true_confusion_matrix[1,cluster] += 1
  		# else:
  		# 	true_confusion_matrix[0,cluster] += 1
  
  	return true_confusion_matrix
  
  def computeF1_macro(confusion_matrix,matching, num_clusters):
  	"""
  	computes the macro F1 score
  	confusion matrix : requres permutation
  	matching according to which matrix must be permuted
  	"""
  	##Permute the matrix columns
  	permuted_confusion_matrix = np.zeros([num_clusters,num_clusters])
  	for cluster in range(num_clusters):
  		matched_cluster = matching[cluster]
   		permuted_confusion_matrix[:,cluster] = confusion_matrix[:,matched_cluster]
   	##Compute the F1 score for every cluster
   	F1_score = 0
   	for cluster in range(num_clusters):
   		TP = permuted_confusion_matrix[cluster,cluster]
   		FP = np.sum(permuted_confusion_matrix[:,cluster]) - TP
   		FN = np.sum(permuted_confusion_matrix[cluster,:]) - TP
   		precision = TP/(TP + FP)
   		recall = TP/(TP + FN)
   		f1 = stats.hmean([precision,recall])
   		F1_score += f1
   	F1_score /= num_clusters
   	return F1_score
  
  def computeNetworkAccuracy(matching,train_cluster_inverse, num_clusters):
  	"""
  	Takes in the matching for the clusters
  	takes the computed clusters
  	computes the average F1 score over the network
  	"""
  	threshold = 1e-2
  	f1 = 0
  	for cluster in range(num_clusters):
  		true_cluster_cov = np.loadtxt("Inverse Covariance cluster ="+ str(cluster) +".csv", delimiter = ",")
  		matched_cluster = matching[cluster]
  		matched_cluster_cov = train_cluster_inverse[matched_cluster] 
  		(nrow,ncol) = true_cluster_cov.shape
  
  		out_true = np.zeros([nrow,ncol])
  		for i in range(nrow):
  			for j in range(ncol):
  				if np.abs(true_cluster_cov[i,j]) > threshold:
  					out_true[i,j] = 1
  		out_matched = np.zeros([nrow,ncol])
  		for i in range(nrow):
  			for j in range(ncol):
  				if np.abs(matched_cluster_cov[i,j]) > threshold:
  					out_matched[i,j] = 1
  		np.savetxt("Network_true_cluster=" +str(cluster) + ".csv",true_cluster_cov, delimiter = ",")
  		np.savetxt("Network_matched_cluster=" + str(matched_cluster)+".csv",matched_cluster_cov, delimiter = ",")
  
  
  		##compute the confusion matrix
  		confusion_matrix = np.zeros([2,2])
  		for i in range(nrow):
  			for j in range(ncol):
  				confusion_matrix[out_true[i,j],out_matched[i,j]] += 1
  		f1 += computeF1_macro(confusion_matrix, [0,1],2)
  	return f1/num_clusters
  
  ############
  
  if not os.path.exists(os.path.dirname(str_NULL)):
      try:
          os.makedirs(os.path.dirname(str_NULL))
      except OSError as exc: # Guard against race condition
          if exc.errno != errno.EEXIST:
              raise
  
  def hex_to_rgb(value):
  	"""
  	Return (red, green, blue) values
  	Input is hexadecimal color code: #rrggbb.
  	"""
  	lv = len(value)
  	out = tuple(int(value[i:i + lv // 3], 16) for i in range(0, lv, lv // 3))
  	out = tuple([x/256.0 for x in out])
  	return out
  
  color_list = []
  for hex_color in hexadecimal_color_list:
  	rgb_color = hex_to_rgb(hex_color)
  	color_list.append(rgb_color)
  colors = color_list
  
  
  train_cluster_inverse = {}
  log_det_values = {}
  computed_covariance = {}
  cluster_mean_info = {}
  cluster_mean_stacked_info = {}
  old_clustered_points = np.zeros(10)
  for iters in range(maxIters):
  	print("\n\n\nITERATION ###", iters)
  	num_clusters = maxClusters - 1
  
  	if iters == 0:
  		## Now splitting up stuff 
  		## split1 : Training and Test
  		## split2 : Training and Test - different clusters
  		training_percent = 0.90
  		training_idx = np.random.choice(m-num_blocks+1, size=int(m*training_percent),replace = False )
  		##Ensure that the first and the last few points are in
  		training_idx = list(training_idx)
  		if 0 not in training_idx:
  			training_idx.append(0)
  		if m - num_stacked  not in training_idx:
  			training_idx.append(m-num_stacked)
  		training_idx = np.array(training_idx)
  
  		sorted_training_idx = sorted(training_idx)
  		num_test_points = m - len(training_idx)
  		test_idx = []
  		##compute the test indices
  		for point in range(m-num_stacked+1):
  			if point not in sorted_training_idx:
  				test_idx.append(point)
  		sorted_test_idx = sorted(test_idx)
  		# np.savetxt("sorted_training.csv", sorted_training_idx, delimiter = ",")
  		# np.savetxt("sorted_test.csv", sorted_test_idx, delimiter = ",")
  
  		##Stack the complete data
  		complete_Data = np.zeros([m - num_stacked + 1, num_stacked*n])
  		len_data = m
  		for i in range(m - num_stacked + 1):
  			idx = i
  			for k in range(num_stacked):
  				if i+k < len_data:
  					idx_k = i + k
  					complete_Data[i][k*n:(k+1)*n] =  Data[idx_k][0:n]
  		# np.savetxt("Complete_Data_stacked_rand_seed="+str(0)+".csv", complete_Data,delimiter=",")
  
  		##Stack the training data
  		complete_D_train = np.zeros([len(training_idx), num_stacked*n])
  		len_training = len(training_idx)
  		for i in range(len(sorted_training_idx)):
  			idx = sorted_training_idx[i]
  			for k in range(num_stacked):
  				if i+k < len_training:
  					idx_k = sorted_training_idx[i+k]
  					complete_D_train[i][k*n:(k+1)*n] =  Data[idx_k][0:n]
  		# np.savetxt("Data_train_rand_seed="+str(0)+".csv", complete_D_train,delimiter=",")
  		##Stack the test  
  		complete_D_test = np.zeros([len(test_idx), num_stacked*n])
  		len_test = len(test_idx)
  
  		for i in range(len(sorted_test_idx)):
  			idx = sorted_test_idx[i]
  			idx_left = idx -1 
  			while idx_left not in sorted_training_idx:
  				idx_left -= 1
  			point_tr = sorted_training_idx.index(idx_left)
  			complete_D_test[i] = complete_D_train[point_tr]
  			complete_D_test[i][0:n] = Data[idx][0:n]
  		# np.savetxt("Data_test_rand_seed="+str(0)+".csv", complete_D_test,delimiter=",")
  		#####INITIALIZATION!!!
  		gmm = mixture.GaussianMixture(n_components=num_clusters, covariance_type="full")
  		gmm.fit(complete_D_train)
  		clustered_points = gmm.predict(complete_D_train)
  		clustered_points_test = gmm.predict(complete_D_test)
  		gmm_clustered_pts_test = gmm.predict(complete_D_test)
  		gmm_clustered_pts = clustered_points + 0
  
  		gmm_covariances = gmm.covariances_
  		gmm_means = gmm.means_
  
  		##USE K-means
  		kmeans = KMeans(n_clusters = num_clusters,random_state = 0).fit(complete_D_train)
  		clustered_points_kmeans = kmeans.labels_
  		clustered_points_test_kmeans = kmeans.predict(complete_D_test)
  
  		true_confusion_matrix_g = compute_confusion_matrix(num_clusters,gmm_clustered_pts,sorted_training_idx)
  		true_confusion_matrix_k = compute_confusion_matrix(num_clusters,clustered_points_kmeans,sorted_training_idx)
  
  	##Get the train and test points
  	train_clusters = collections.defaultdict(list)
  	test_clusters = collections.defaultdict(list)
  	len_train_clusters = collections.defaultdict(int)
  	len_test_clusters = collections.defaultdict(int)
  
  	counter = 0
  	for point in range(len(clustered_points)):
  		cluster = clustered_points[point]
  		train_clusters[cluster].append(point)
  		len_train_clusters[cluster] += 1
  		counter +=1 
  
  	for point in range(len(clustered_points_test)):
  		cluster = clustered_points_test[point]
  		test_clusters[cluster].append(point)
  		len_test_clusters[cluster] += 1
  		counter +=1 
  
  
  
  	##train_clusters holds the indices in complete_D_train 
  	##for each of the clusters
  	for cluster in range(num_clusters):
  		if len_train_clusters[cluster] != 0:
  			indices = train_clusters[cluster]
  			indices_test = test_clusters[cluster]
  
  
  
  			D_train = np.zeros([len_train_clusters[cluster],num_stacked*n])
  			for i in range(len_train_clusters[cluster]):
  				point = indices[i]
  				D_train[i,:] = complete_D_train[point,:]
  
  			D_test = np.zeros([len_test_clusters[cluster], num_stacked*n])
  			for i in range(len_test_clusters[cluster]):
  				point = indices_test[i]
  				D_test[i,:] = complete_D_test[point,:]
  			print("stacking Cluster #", cluster,"DONE!!!")
  			##Fit a model - OPTIMIZATION	
  			size_blocks = n
  			probSize = num_stacked * size_blocks
  			lamb = np.zeros((probSize,probSize)) + lam_sparse
  			S = np.cov(np.transpose(D_train) )
  
  			print("starting the OPTIMIZATION for cluster#", cluster)
  			#Set up the Toeplitz graphical lasso problem
  			gvx = TGraphVX()
  			theta = semidefinite(probSize,name='theta')
  			obj = -log_det(theta) + trace(S*theta)
  			gvx.AddNode(0, obj)
  			gvx.AddNode(1)
  			dummy = Variable(1)
  			gvx.AddEdge(0,1, Objective = lamb*dummy + num_stacked*dummy + size_blocks*dummy)
  			
  			##solve using customized ADMM solver
  			gvx.Solve(Verbose=False, MaxIters=1000, Rho = 1, EpsAbs = 1e-6, EpsRel = 1e-6)
  			print("\ncompleted solving the optimization problem for the cluster")
  
  			#THIS IS THE SOLUTION
  			val = gvx.GetNodeValue(0,'theta')
  			S_est = upper2Full(val, 0)
  			X2 = S_est
  			u, _ = np.linalg.eig(S_est)
  			cov_out = np.linalg.inv(X2)
  
  			inv_matrix = cov_out
  			##Store the log-det, covariance, inverse-covariance, cluster means, stacked means
  			log_det_values[num_clusters, cluster] = np.log(np.linalg.det(cov_out))
  			computed_covariance[num_clusters,cluster] = cov_out
  			cluster_mean_info[num_clusters,cluster] = np.mean(D_train, axis = 0)[(num_stacked-1)*n:num_stacked*n].reshape([1,n])
  			cluster_mean_stacked_info[num_clusters,cluster] = np.mean(D_train,axis=0)
  			train_cluster_inverse[cluster] = X2
  
  	cluster_norms = list(np.zeros(num_clusters))
  	print("printing the cluster len")
  	for cluster in range(num_clusters):
  		print("length of the cluster ", cluster,"------>", len_train_clusters[cluster])
  	##Computing the norms
  	if iters != 0:
  		for cluster in range(num_clusters):
  			cluster_norms[cluster] = (np.linalg.norm(old_computed_covariance[num_clusters,cluster]),cluster)
  		sorted_cluster_norms = sorted(cluster_norms,reverse = True)
  
  	##Add a point to the empty clusters 
  	##Assumption more non empty clusters than empty ones
  	counter = 0
  	for cluster in range(num_clusters):
  		if len_train_clusters[cluster] == 0:
  			##Add a point to the cluster
  			while len_train_clusters[sorted_cluster_norms[counter][1]] == 0:
  				# print "counter is:", counter
  				counter += 1
  				counter = counter % num_clusters
  				# print "counter is:", counter
  
  			cluster_selected = sorted_cluster_norms[counter][1]
  			# print "cluster that is zero is:", cluster, "selected cluster instead is:", cluster_selected
  			break_flag = False
  			while not break_flag:
  				point_num = random.randint(0,len(clustered_points))
  				if clustered_points[point_num] == cluster_selected:
  					clustered_points[point_num] = cluster
  					computed_covariance[num_clusters,cluster] = old_computed_covariance[num_clusters,cluster_selected]
  					cluster_mean_stacked_info[num_clusters,cluster] = complete_D_train[point_num,:]
  					cluster_mean_info[num_clusters,cluster] = complete_D_train[point,:][(num_stacked-1)*n:num_stacked*n]
  					break_flag = True
  			counter += 1
  
  	old_train_clusters = train_clusters
  	old_computed_covariance = computed_covariance
  	# print "UPDATED THE OLD COVARIANCE"
  
  
  
  	##Code -----------------------SMOOTHENING
  	##For each point compute the LLE 
  	print("beginning with the DP - smoothening ALGORITHM")
  
  	LLE_all_points_clusters = np.zeros([len(clustered_points),num_clusters])
  	for point in range(len(clustered_points)):
  		# print "Point #", point
  		if point + num_stacked-1 < complete_D_train.shape[0]:
  			for cluster in range(num_clusters):
  				# print "\nCLuster#", cluster
  				cluster_mean = cluster_mean_info[num_clusters,cluster] 
  				cluster_mean_stacked = cluster_mean_stacked_info[num_clusters,cluster] 
  
  				x = complete_D_train[point,:] - cluster_mean_stacked[0:(num_blocks-1)*n]
  				cov_matrix = computed_covariance[num_clusters,cluster][0:(num_blocks-1)*n,0:(num_blocks-1)*n]
  				inv_cov_matrix = np.linalg.inv(cov_matrix)
  				log_det_cov = np.log(np.linalg.det(cov_matrix))# log(det(sigma2|1))
  				lle = np.dot(   x.reshape([1,(num_blocks-1)*n]), np.dot(inv_cov_matrix,x.reshape([n*(num_blocks-1),1]))  ) + log_det_cov
  				LLE_all_points_clusters[point,cluster] = lle
  	
  	##Update cluster points - using dynamic programming smoothening
  	clustered_points = updateClusters(LLE_all_points_clusters,switch_penalty = switch_penalty)
  	print("\ncompleted smoothening algorithm")
  	print("\n\nprinting the length of points in each cluster")
  
  	for cluster in range(num_clusters):
  		print("length of cluster #", cluster, "-------->", sum([x== cluster for x in clustered_points]))
  	true_confusion_matrix = np.zeros([num_clusters,num_clusters])
  
  	##Save a figure of segmentation
  	plt.figure()
  	plt.plot(sorted_training_idx[0:len(clustered_points)],clustered_points,color = "r")#,marker = ".",s =100)
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("TRAINING_EM_lam_sparse="+str(lam_sparse) + "switch_penalty = " + str(switch_penalty) + ".jpg")
  	plt.close("all")
  
  	plt.figure()
  	plt.plot(sorted_training_idx[0:len(gmm_clustered_pts)],gmm_clustered_pts,color = "r")#,marker=".",s=100)
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("TRAINING_GMM_lam_sparse="+str(lam_sparse) + "switch_penalty =" + str(switch_penalty)+ ".jpg")
  	plt.close("all")
  	# print "Done writing the figure"
  
  	# print "Done writing the figure"
  
  	true_confusion_matrix = compute_confusion_matrix(num_clusters,clustered_points,sorted_training_idx)
  	# print "TRAINING TRUE confusion MATRIX:\n", true_confusion_matrix
  
  	####TEST SETS STUFF
  	### The closest point in training set is the cluster
  	### LLE + swtiching_penalty
  	clustered_test = np.zeros(len(clustered_points_test))
  	for point in range(len(clustered_points_test)):
  		idx = sorted_test_idx[point]
  		##Get the 2 closest points from training
  		idx1 = idx + 1
  		while (idx1 not in sorted_training_idx and idx1 < m):
  			idx1 += 1
  		idx2 = idx -1 
  		while (idx2 not in sorted_training_idx and idx2 > -1):
  			idx2 -= 1
  		if idx1 == m or idx2 == -1:
  			print("idx1 :", idx1 == m)
  			print("idx2 :", idx2 == -1)
  			print("point is:", point)
  			print("idx of the point is:", idx)
  			print("control should NOT reach here!!!!!!!")
  			break
  		vals = np.zeros(num_clusters)
  		right_clust = clustered_points[sorted_training_idx.index(idx1)]
  		left_clust = clustered_points[sorted_training_idx.index(idx2)]
  		point_tr = sorted_training_idx.index(idx2)
  		data_tt = complete_D_train[point_tr,:]
  		data_tt[0:n] = Data[idx,:]#complete_D_test[point,0:num_stacked] 
  
  		for cluster in range(num_clusters):
  			cluster_mean = cluster_mean_info[num_clusters,cluster] 
  			cluster_mean_stacked = cluster_mean_stacked_info[num_clusters,cluster] 
  
  			x = data_tt - cluster_mean_stacked[0:(num_blocks-1)*n]
  			cov_matrix = computed_covariance[num_clusters,cluster][0:(num_blocks-1)*n,0:(num_blocks-1)*n]
  			inv_cov_matrix = np.linalg.inv(cov_matrix)
  			log_det_cov = np.log(np.linalg.det(cov_matrix))# log(det(sigma2|1))
  			lle = np.dot(   x.reshape([1,(num_blocks-1)*n]), np.dot(inv_cov_matrix,x.reshape([n*(num_blocks-1),1]))  ) + log_det_cov
  			vals[cluster] = lle + switch_penalty*(cluster !=left_clust ) + switch_penalty*(cluster != right_clust )
  		out = np.argmin(vals)
  		clustered_test[point] = out
  
  	plt.figure()
  	plt.plot(sorted_test_idx[0:len(clustered_test)],clustered_test,color = "r")#,marker = ".", s= 100)
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("TEST_EM_lam_sparse="+str(lam_sparse) +"switch_penalty="+str(switch_penalty)+ ".jpg")
  	plt.close("all")
  	# print "done writing"
  
  
  	##GMM - TEST PREDICTIONS
  	clustered_test_gmm = np.zeros(len(clustered_points_test))
  	for point in range(len(clustered_points_test)):
  		idx = sorted_test_idx[point]
  		##Get the 2 closest points from training
  		idx1 = idx + 1
  		while (idx1 not in sorted_training_idx and idx1 < m):
  			idx1 += 1
  		idx2 = idx -1 
  		while (idx2 not in sorted_training_idx and idx2 > -1):
  			idx2 -= 1
  		if idx1 == m or idx2 == -1:
  			print("idx1 :", idx1 == m)
  			print("idx2 :", idx2 == -1)
  			print("point is:", point)
  			print("idx of the point is:", idx)
  
  			print("SOMETHING WRONG!!!!!!!")
  			break
  		vals = np.zeros(num_clusters)
  		right_clust = gmm_clustered_pts[sorted_training_idx.index(idx1)]
  		left_clust = gmm_clustered_pts[sorted_training_idx.index(idx2)]
  		point_tr = sorted_training_idx.index(idx2)
  		data_tt = complete_D_train[point_tr,:]
  		data_tt[0:n] = Data[idx,:]#complete_D_test[point,0:num_stacked] 
  
  		for cluster in range(num_clusters):
  			cluster_mean_stacked = gmm_means[cluster] 
  
  			x = data_tt - cluster_mean_stacked[0:(num_blocks-1)*n]
  			cov_matrix = gmm_covariances[cluster][0:(num_blocks-1)*n,0:(num_blocks-1)*n]
  			inv_cov_matrix = np.linalg.inv(cov_matrix)
  			log_det_cov = np.log(np.linalg.det(cov_matrix))
  			lle = np.dot(   x.reshape([1,(num_blocks-1)*n]), np.dot(inv_cov_matrix,x.reshape([n*(num_blocks-1),1]))  ) + log_det_cov
  			vals[cluster] = lle #+ switch_penalty*(cluster !=left_clust ) + switch_penalty*(cluster != right_clust )
  		out = np.argmin(vals)
  		clustered_test_gmm[point] = out
  
  	plt.figure()
  	plt.plot(sorted_test_idx[0:len(clustered_test_gmm)],clustered_test_gmm,color = "r")#,marker = ".", s= 100)
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("TEST_GMM_lam_sparse="+str(lam_sparse) + ".jpg")
  	plt.close("all")
  	# print "done writing"
  
  	plt.figure()
  	plt.plot(sorted_test_idx[0:len(clustered_points_test_kmeans)],clustered_points_test_kmeans,color = "r")#,marker = ".", s= 100)
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("TEST_Modified_KMEANS_NEW_lam_sparse="+str(lam_sparse) + ".jpg")
  	plt.close("all")
  	# print "done writing"
  
  	##Segment length
  	seg_len = 50
  	true_confusion_matrix_EM = compute_confusion_matrix(num_clusters,clustered_test,sorted_test_idx)
  	true_confusion_matrix_GMM = compute_confusion_matrix(num_clusters,gmm_clustered_pts_test,sorted_test_idx)
  	true_confusion_matrix_kmeans = compute_confusion_matrix(num_clusters,clustered_points_test_kmeans,sorted_test_idx)
  
  
  	true_answers = np.zeros(len(clustered_points))
  	for point in range(len(clustered_points)):
  		num = int(sorted_training_idx[point]/25.0)
  		if num <10 :
  			cluster = 0
  		elif num < 20:
  			cluster = 1
  		else:
  			cluster = 0
  		true_answers[point] = cluster
  
  	plt.figure()
  	plt.plot(sorted_training_idx[0:len(clustered_points)],true_answers,color = "k")
  	plt.ylim((-0.5,num_clusters + 0.5))
  	# plt.savefig("True Output modifed lam_sparse=" + str(lam_sparse)+ ".jpg")
  	plt.close("all")
  	binary_EM = (true_confusion_matrix_EM[0,0] + true_confusion_matrix_EM[1,1])/len(clustered_points_test)
  	binary_EM = np.max([binary_EM,1 -binary_EM])
  	# print "EM is -------->",binary_EM
  	binary_GMM = (true_confusion_matrix_GMM[0,0] + true_confusion_matrix_GMM[1,1])/len(clustered_points_test)
  	binary_GMM = np.max([binary_GMM,1-binary_GMM])
  
  	binary_Kmeans = (true_confusion_matrix_kmeans[0,0] + true_confusion_matrix_kmeans[1,1])/len(clustered_points_test)
  	binary_Kmeans = np.max([binary_Kmeans,1-binary_Kmeans])
  
  
  	##Create the F1 score from the graphs from k-means and GMM
  	##Get the train and test points
  	train_inverse_covariance_kmeans = {}
  	train_inverse_covariance_gmm = {}
  
  	counter = 0
  	for cluster in range(num_clusters):
  		##GMM
  		out = [(x == cluster) for x in gmm_clustered_pts]
  		len_cluster = sum(out)
  		D_train = np.zeros([len_cluster,num_stacked*n])
  		counter = 0
  		for point in range(len(gmm_clustered_pts)):
  			if gmm_clustered_pts[point] == cluster:
  				D_train[counter,:] = complete_D_train[point,:]
  				counter += 1
  
  		train_inverse_covariance_gmm[cluster] = np.linalg.inv(gmm_covariances[cluster])
  
  		##Kmeans
  		out = [(x == cluster) for x in clustered_points_kmeans]
  		len_cluster = sum(out)
  		D_train = np.zeros([len_cluster,num_stacked*n])
  		counter2 = 0
  		for point in range(len(clustered_points_kmeans)):
  			if clustered_points_kmeans[point] == cluster:
  				D_train[counter2,:] = complete_D_train[point,:]
  				counter2 += 1
  
  		train_inverse_covariance_kmeans[cluster] = X2
  
  	##Using inverses from kmeans, GMM and EM
  	threshold = 1e-5
  	##Thresholding
  	threshold_kmeans = {}
  	threshold_GMM = {}
  	threshold_EM = {}
  	threshold_actual = {}
  
  	##Change this to add thresholding function
  	##Kmeans - thresholding
  	for cluster in range(num_clusters):
  	    out = np.zeros(train_inverse_covariance_kmeans[0].shape, dtype = np.int)
  	    A = train_inverse_covariance_kmeans[cluster]
  	    for i in range(out.shape[0]):
  	        for j in range(out.shape[1]):
  				if (np.abs(A[i,j]) > threshold):
  					out[i,j] = 1
  		threshold_kmeans[cluster] = out
  
  	##GMM - thresholding
  	for cluster in range(num_clusters):
  		out = np.zeros(train_inverse_covariance_gmm[0].shape, dtype = np.int)
  		A = train_inverse_covariance_gmm[cluster]
  		for i in range(out.shape[0]):
  			for j in range(out.shape[1]):
  				if np.abs(A[i,j]) > threshold:
  					out[i,j] = 1
  		threshold_GMM[cluster] = out
  
      ## EM - thresholding
  	for cluster in range(num_clusters):
  		out = np.zeros(train_inverse_covariance_gmm[0].shape, dtype = np.int)
  		A = train_cluster_inverse[cluster]
  		for i in range(out.shape[0]):
  			for j in range(out.shape[1]):
  				if np.abs(A[i,j]) > threshold:
  					out[i,j] = 1
  		threshold_EM[cluster] = out
  
      ##compute the matching
      ##Assume its a 2x2 matrix?
  	actual_clusters = {}
  	for cluster in range(num_clusters):
  		print("getting the actual Inverse covariances")
  		actual_clusters[cluster] = np.loadtxt("Inverse Covariance cluster =" + str(cluster)+".csv", delimiter = ",")
  
      ##compute the appropriate matching
  	matching_Kmeans = find_matching(true_confusion_matrix_kmeans)
  	matching_GMM = find_matching(true_confusion_matrix_GMM)
  	matching_EM = find_matching(true_confusion_matrix_EM)
  
  	correct_EM = 0
  	correct_GMM = 0
  	correct_KMeans = 0
  	for cluster in range(num_clusters):
  		matched_cluster_EM = matching_EM[cluster]
  		matched_cluster_GMM = matching_GMM[cluster]
  		matched_cluster_Kmeans = matching_Kmeans[cluster]
  
  		correct_EM += true_confusion_matrix_EM[cluster,matched_cluster_EM]
  		correct_GMM += true_confusion_matrix_GMM[cluster,matched_cluster_GMM]
  		correct_KMeans += true_confusion_matrix_kmeans[cluster, matched_cluster_Kmeans]
  		# np.savetxt("computed estimated_matrix cluster =" + str(cluster) + ".csv", train_cluster_inverse[matched_cluster] , delimiter = ",", fmt = "%1.6f")
  	binary_EM = correct_EM/len(clustered_points_test)
  	binary_GMM = correct_GMM/len(clustered_points_test)
  	binary_Kmeans = correct_KMeans/len(clustered_points_test)
  
  
  	# print "\n\nKMEANS"
  	# print "true confusion_matrix", true_confusion_matrix_kmeans
  	# print "matching", matching_Kmeans
  
  	# print "\n\nGMM"
  	# print "true_confusion_matrix_GMM", true_confusion_matrix_GMM
  	# print "matching GMM", matching_GMM
  
  	# print "\n\nEM"
  	# print "true confusion_matrix", true_confusion_matrix_EM
  	# print "matching ", matching_EM
  
  
  	print("\n\n\n")
  	# print "The F1 scores is:", F1_EM,F1_GMM, F1_Kmeans
  	print("The binary accuracy at the end of iteration is:", binary_EM, binary_GMM, binary_Kmeans)
  
  	if np.array_equal(old_clustered_points,clustered_points):
  		print("\n\n\n\nCONVERGED!!! BREAKING EARLY!!!")
  		break
  	old_clustered_points = clustered_points
  
  ##Training confusion matrix
  train_confusion_matrix_EM = compute_confusion_matrix(num_clusters,clustered_points,sorted_training_idx)
  train_confusion_matrix_GMM = compute_confusion_matrix(num_clusters,gmm_clustered_pts,sorted_training_idx)
  train_confusion_matrix_kmeans = compute_confusion_matrix(num_clusters,clustered_points_kmeans,sorted_training_idx)
  test_confusion_matrix_EM = compute_confusion_matrix(num_clusters,clustered_test,sorted_test_idx)
  
  out = -1#computeNetworkAccuracy(matching_EM, train_cluster_inverse,num_clusters)
  
  # print "the NETWORK F1 score is:", out
  
  
  # print "train matrices are: EM\n", train_confusion_matrix_EM
  # print "\n GMM", train_confusion_matrix_GMM
  # print "\n kmeans", train_confusion_matrix_kmeans
  # print " \n"
  f1_EM_tr = computeF1_macro(train_confusion_matrix_EM,matching_EM,num_clusters)
  f1_GMM_tr = computeF1_macro(train_confusion_matrix_GMM,matching_GMM,num_clusters)
  f1_kmeans_tr = computeF1_macro(train_confusion_matrix_kmeans,matching_Kmeans,num_clusters)
  f1_EM_test = computeF1_macro(test_confusion_matrix_EM,matching_EM,num_clusters)
  # print "The TEST binary accuracy", binary_EM, binary_GMM, binary_Kmeans
  print("\n\n")
  print("FINAL TRAINING MACRO F1 score:", f1_EM_tr, f1_GMM_tr, f1_kmeans_tr)
  # print "TEST F1 score:", f1_EM_test
  correct_EM = 0
  correct_GMM = 0
  correct_KMeans = 0
  for cluster in range(num_clusters):
  	matched_cluster_EM = matching_EM[cluster]
  	matched_cluster_GMM = matching_GMM[cluster]
  	matched_cluster_Kmeans = matching_Kmeans[cluster]
  
  	correct_EM += train_confusion_matrix_EM[cluster,matched_cluster_EM]
  	correct_GMM += train_confusion_matrix_GMM[cluster,matched_cluster_GMM]
  	correct_KMeans += train_confusion_matrix_kmeans[cluster, matched_cluster_Kmeans]
  binary_EM = correct_EM/len(training_idx)
  binary_GMM = correct_GMM/len(training_idx)
  binary_Kmeans = correct_KMeans/len(training_idx)
  
  print("\n\n")
  print("\n\n\n")
  # print "The TRAINING binary accuracy", binary_EM, binary_GMM, binary_Kmeans
  # print "lam_sparse", lam_sparse
  # print "switch_penalty", switch_penalty
  # print "num_cluster", maxClusters - 1
  # print "num stacked", num_stacked
  
  
  
  
  
  
  
  
  