text:
''' sparseae.py

	Implementation of Sparse Autoencoder

	Author: Cem Karaoguz
	Date: 13.03.2015
	Version: 1.0
'''

import sys
import numpy as np
import pylab as pl
from scipy.optimize import minimize
import scipy.io
import scipy.linalg

from UFL.common import DataInputOutput, DataNormalization, AuxFunctions, Visualization
from UFL.PCA import PCA

ACTIVATION_FUNCTION_SIGMOID = 0;
ACTIVATION_FUNCTION_IDENTITY = 1;
SPARSEAE_ACTIVATION_FUNCTIONS = [ACTIVATION_FUNCTION_SIGMOID, ACTIVATION_FUNCTION_IDENTITY]

class SparseAutoencoder:
	'''
	Sparse Autoencoder

	Implemented as 3 layered neural network
	'''

	def __init__(self,
	             dimLayers,
				 lambda_w=0.0001,
				 beta=3,
				 sparsityParam=0.01,
				 actFunctions=[ACTIVATION_FUNCTION_SIGMOID, ACTIVATION_FUNCTION_SIGMOID],
				 debug=0):
		'''
		Initialization function of the Sparse Autoencoder class

		Arguments
		dimLayers		: Dimensions of the layers, must be in the form [Input layer dim., hidden layer dim., output layer dim.] where input dim. = output dim.
		lambda_w		: weight decay parameter, default is 0.0001
		beta			: weight of sparsity penalty term, default is 3
		sparsityParam	: weight of the sparsity in the cost function, default is 0.01
		actFunctions	: activation function, default is [ACTIVATION_FUNCTION_SIGMOID, ACTIVATION_FUNCTION_SIGMOID]
		debug			: Debugging flag
		'''
		self.isInitialized = False;

		self.debug = debug;
		self.dimLayers = dimLayers;
		self.lambda_w = lambda_w;
		self.beta = beta;
		self.sparsityParam = sparsityParam;
		self.actFunctions = actFunctions;

		assert self.dimLayers[0]>0, 'ERROR:SparseAutoencoder:init: Input size must be >0'
		assert self.dimLayers[1]>0, 'ERROR:SparseAutoencoder:init: Feature size must be >0'
		assert self.dimLayers[0]==self.dimLayers[-1], 'ERROR:SparseAutoencoder:init: Input and output size should be the same'

		self.nLayers = len(self.dimLayers);

		assert self.nLayers==3, 'ERROR:SparseAutoencoder:init: Autoencoder should be 3 layered network'
		assert len(self.actFunctions)==(self.nLayers-1), 'ERROR:SparseAutoencoder:init: Please give activation functions for hidden and output layers'
		for f in self.actFunctions:
			assert f in SPARSEAE_ACTIVATION_FUNCTIONS, 'ERROR:SparseAutoencoder:init: Activation function not recognized'

		weights = [];
		biases = [];
		self.weightPrototypes = [];
		self.biasPrototypes = [];
		self.sizeParams = 0;
		# Choose weights uniformly from the interval [-r, r]
		r  = np.sqrt(6) / np.sqrt(self.dimLayers[0]+self.dimLayers[1]+1);
		for i in range(self.nLayers-1):
			weights.append( np.random.rand(self.dimLayers[i+1], self.dimLayers[i]) * 2 * r - r );
			biases.append( np.zeros((self.dimLayers[i+1], 1)) );
			# Set network topology
			self.weightPrototypes.append((self.dimLayers[i+1], self.dimLayers[i]));
			self.biasPrototypes.append((self.dimLayers[i+1], 1));
			# Total number of parameters
			self.sizeParams = self.sizeParams + (self.dimLayers[i+1] * self.dimLayers[i]) + (self.dimLayers[i+1]);

		self.params = self.rollParameters(weights, biases);

		if debug:
			print 'DEBUG:SparseAutoencoder:init: initialized for inputSize: ', self.dimLayers[0];
			print 'DEBUG:SparseAutoencoder:init: initialized for featureSize: ', self.dimLayers[1];
			print 'DEBUG:SparseAutoencoder:init: initialized for lambda_w: ', self.lambda_w;
			print 'DEBUG:SparseAutoencoder:init: initialized for beta: ', self.beta;
			print 'DEBUG:SparseAutoencoder:init: initialized for sparsityParam: ', self.sparsityParam;
			print

		self.isInitialized = True;

	def rollParameters(self, weights, biases):
		'''
		Converts the parameters in matrix form into vector

		Arguments
		weights	: list of weight matrices of each layer
		biases	: list of bias vectors of each layer

		Returns
		params	: parameter vector
		'''
		assert AuxFunctions.checkNetworkParameters(weights, self.weightPrototypes), 'ERROR:SparseAutoencoder:rollParameters: weight dimension does not match the network topology' ;
		assert AuxFunctions.checkNetworkParameters(biases, self.biasPrototypes), 'ERROR:SparseAutoencoder:rollParameters: bias dimension does not match the network topology';

		params = np.array([]);
		for i in range(len(weights)):
			params = np.hstack((params, weights[i].flatten(), biases[i].flatten()))

		return params

	def unrollParameters(self, params):
		'''
		Converts the vectorized parameters into list of matrices

		Arguments
		params	: parameter vector

		Returns
		weights	: list of weight matrices of each layer
		biases	: list of bias vectors of each layer
		'''
		assert len(params)==self.sizeParams, 'ERROR:SparseAutoencoder:unrollParameters: Parameter size mismatch'

		weights = [];
		biases = [];
		read_start = 0;
		read_end = 0;

		for i in range(self.nLayers - 1):
			# set the end index for read
			read_end = read_start + self.dimLayers[i+1]*self.dimLayers[i];
			# read the weights for the current layer
			w = params[read_start:read_end];
			# reshape the weights
			w = np.reshape(w, (self.dimLayers[i+1], self.dimLayers[i]))
			weights.append( w );
			# set the start index for the next read
			read_start = read_end;
			# set the end index for the next read
			read_end = read_start + self.dimLayers[i+1];
			# read the bias terms
			b = params[read_start:read_end];
			# reshape and store the bias
			biases.append( np.reshape(b, (self.dimLayers[i+1], 1)) )
			# set the start index for the next read
			read_start = read_end;

		return weights, biases;

	def getWeights(self):
		'''
		Returns the list of weight matrices of the network
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getWeights: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);

		return weights

	def getBiases(self):
		'''
		Returns the list of bias vectors of the network
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getBiases: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);

		return biases

	def getParameters(self):
		'''
		Returns the model parameters of the network
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getParameters: The instance is not properly initialized'

		return  self.params

	def setWeights(self, weights_new):
		'''
		Updates the weights of the model parameters of the network

		Arguments
		weights_new	: New weights to set
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setWeights: The instance is not properly initialized'
		assert AuxFunctions.checkNetworkParameters(weights_new, self.weightPrototypes), 'ERROR:SparseAutoencoder:setWeights: weight dimension does not match the network topology' ;

		[weights, biases] = self.unrollParameters(self.params);
		weights = weights_new;
		self.params = self.rollParameters(weights, biases);

	def setBiases(self, biases_new):
		'''
		Updates the biases of the model parameters of the network

		Arguments
		biases_new	: New biases to set
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setBiases: The instance is not properly initialized'
		assert AuxFunctions.checkNetworkParameters(biases_new, self.biasPrototypes), 'ERROR:SparseAutoencoder:setBiases: bias dimension does not match the network topology' ;

		[weights, biases] = self.unrollParameters(self.params);
		biases = biases_new;
		self.params = self.rollParameters(weights, biases);

	def setParameters(self, weights, biases):
		'''
		Updates the internal Sparse Autoencoder parameters with the given ones

		Arguments
		weights	: list of weights to set for the the network
		biases	: list of biases to set for the network
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setParams: The instance is not properly initialized'

		self.params = self.rollParameters(weights, biases);

	def getParameterSize(self):
		'''
		Returns the size of the model parameters i.e. (weights and biases) of the network
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getParameterSize: The instance is not properly initialized'

		return self.sizeParams

	def doForwardPropagation(self, X, weights, biases):
		'''
		Computes the forward propagation of the input in the network:

		Z{l+1} = W{l}*H{l} + B{l}
		H{l+1} = f(Z{l+1})

		where
		{l} and {l+1} denote layers,
		B is the bias matrix, columnwise repetition of the bias vector with the number of samples,
		Z is the output matrix of neurons before the activation function is applied,
		f(.) is the activation function
		H is the output matrix of neurons after the activation function is applied (h{1}=X),

		Arguments
		X			: data matrix in the form [input dim., number of samples]
		weights		: list of weight matrices of each layer
		biases		: list of bias vectors of each layer

		Returns
		outputs		: list of output matrices (z) of each layer (output of neuron before activation function)
		activities	: list of activation matrices (h) of each layer (output of neuron after activation function)
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:doForwardPropagation: The instance is not properly initialized'
		assert AuxFunctions.checkNetworkParameters(weights, self.weightPrototypes), 'ERROR:SparseAutoencoder:doForwardPropagation: weight dimension does not match the network topology' ;
		assert AuxFunctions.checkNetworkParameters(biases, self.biasPrototypes), 'ERROR:SparseAutoencoder:doForwardPropagation: bias dimension does not match the network topology';

		outputs = [];
		activities = [];
		for layer in range(self.nLayers-1):

			if layer==0:
				x = X;
			else:
				x = activities[layer-1];

			z = np.dot(weights[layer], x) + np.repeat(biases[layer], x.shape[1], 1);

			if self.actFunctions[layer]==ACTIVATION_FUNCTION_SIGMOID:
				h = AuxFunctions.sigmoid(z);
			elif self.actFunctions[layer]==ACTIVATION_FUNCTION_IDENTITY:
				h = 1.0 * z;
			else:
				# Should not be here
				print 'ERROR:SparseAutoencoder:doForwardPropagation: Activation function not recognized'
				sys.exit()

			outputs.append(z);
			activities.append(h);

		return [outputs, activities];

	def computeCost(self, theta, X):
		'''
		Computes the value of the Sparse Autoencoder objective function for given
		features (theta), data matrix (X) and corresponding labels (y):

		f = 1/2 * sum((X - H)^2) + beta * sum_j(KL(rho||rho_hat_j))

		where H is the activity matrix with columns corresponding to the samples and
		rows corresponding to the output layer units,
		rho is the target average activity value,
		rho_hat_j is the average activity of unit j in the hidden layer:

		rho_hat_j = 1/nSamples * sum_i(h_j(x_i)),		i=1...nSamples

		h_j(x_i) is the activity of hidden unit j given input x_i,
		KL(.) is Kullback-Leibler divergence function:

		KL(rho||rho_hat) = rho*log(rho/rho_hat) + (1-rho)*log((1-rho)/(1-rho_hat)).

		KL function penalizes any average values that diverges from rho	hence, sparsity
		in the hidden units is encouraged.

		Arguments
		theta	: function parameters in the form (feature dim * input dim, )
		X		: data matrix in the form [input dim, number of samples]
		y		: labels in the form [1, number of samples].

		Returns
		f		: computed cost (floating point number)
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:computeCost: The instance is not properly initialized'

		f = 0;

		nSamples = X.shape[1];

		[weights, biases] = self.unrollParameters(theta);

		[outputs, activities] = self.doForwardPropagation(X, weights, biases);

		cost_fidelty = (1.0/nSamples) * np.sum( 0.5 * (X - activities[-1])**2 );

		W_sum = 0;
		for l in range(len(weights)):
			W_sum = W_sum + np.sum((weights[l])**2);

		cost_regularization = (self.lambda_w/2) * W_sum;

		a2_mean = np.mean(activities[0], 1);
		cost_sparsity = self.beta * np.sum(self.sparsityParam * np.log(self.sparsityParam/a2_mean) + (1 - self.sparsityParam) * np.log( (1 - self.sparsityParam)/(1 - a2_mean) ));

		f = cost_fidelty + cost_regularization + cost_sparsity;

		return f

	def computeGradient(self, theta, X):
		'''
		Computes gradients of the Sparse Autoencoder objective function for given parameters,
		data and corresponding labels using error back propagation:

		E_{3} = (H_{3} - X)
		E_{2} = (W_{2}'*E_{3}  + beta*(-rho/rho_hat + (1-rho)/(1-rho_hat)) ) * df(Z_{2})/dz

		where rho is the target average activity value,
		rho_hat is the average activity vector of the units in the hidden layer,
		df(z)/dz is the derivatives of acitvation function f(.) at points Z which is

		df(Z)/dz = f(Z)*(1-f(Z)) 	for sigmoid activation function and
		df(Z)/dz = Z 				for identity activation function.

		Gradients are computed via:

		dJ(W,b;X,y)/dW_{l} = E_{l+1} * H_{l}'
		dJ(W,b;X,y)/db_{l} = sum(E_{l+1})

		where sum(.) is taken columnwise i.e. over samples

		Arguments
		theta	: function parameters in the form (feature dim * input dim, )
		X		: data matrix in the form [input dim, number of samples]
		y		: labels in the form [1, number of samples].
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:computeGradient: The instance is not properly initialized'

		nSamples = X.shape[1];

		[weights, biases] = self.unrollParameters(theta);

		[outputs, activities] = self.doForwardPropagation(X, weights, biases);

		# Output layer
		if self.actFunctions[-1]==ACTIVATION_FUNCTION_SIGMOID:
			delta3 = 1.0/nSamples * (activities[-1] - X) * ( activities[-1] * (1 - activities[-1]) );
		elif self.actFunctions[-1]==ACTIVATION_FUNCTION_IDENTITY:
			delta3 = 1.0/nSamples * (activities[-1] - X);
		else:
			# Should not be here
			print 'ERROR:SparseAutoencoder:doForwardPropagation: Activation function not recognized'
			sys.exit()

		# Hidden layer
		if self.actFunctions[-2]==ACTIVATION_FUNCTION_SIGMOID:
			aux1 = np.dot(np.transpose(weights[1]), delta3);
			a2_mean = np.mean(activities[0], 1);
			aux2 = self.beta/(nSamples*1.0) * ( (1-self.sparsityParam)/(1-a2_mean) - (self.sparsityParam/a2_mean) );
			aux3 = np.repeat(np.reshape(aux2, [len(aux2), 1]), nSamples, 1 );
			delta2 = (aux1 + aux3) * ( activities[-2] * (1 - activities[-2]) );
		else:
			# Should not be here
			print 'ERROR:SparseAutoencoder:doForwardPropagation: Activation function not recognized'
			sys.exit()

		deltas = [delta2, delta3];
		gradients_W = [];
		gradients_b = [];
		for layer in range(len(weights)):
			if layer==0:
				x_in = X;
			else:
				x_in = activities[layer-1];

			gradients_W.append( np.dot(deltas[layer], np.transpose(x_in)) + (self.lambda_w * weights[layer]) );
			gradients_b.append( np.sum(deltas[layer], 1) );

		return self.rollParameters(gradients_W, gradients_b);

	def testGradient(self, X):
		'''
		Tests the analytical gradient computation by comparing it with the numerical gradients

		Arguments
		X		: data matrix the form [input dim., number of samples]

		Returns
		result	: 0 if passed, -1 if failed
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:testGradient: The instance is not properly initialized'
		assert X.shape[0]==self.dimLayers[0], 'ERROR:SparseAutoencoder:testGradient: Dimensions of given data do not match with the number of parameters'

		if self.debug: print 'DEBUG:SparseAutoencoder:testGradient: Testing gradient computation...'

		result = 0;

		grad = self.computeGradient(self.params, X);

		numGrad = AuxFunctions.computeNumericalGradient( func=self.computeCost, params=self.params, args=((X,)) );

		errorGrad = np.sqrt(np.sum((grad - numGrad)**2));

		if errorGrad<1e-4:
			if self.debug:
				print 'DEBUG:SparseAutoencoder:testGradient:Gradient error: ', errorGrad
				print 'DEBUG:SparseAutoencoder:testGradient:Gradient check PASSED!'
				print

			result = 0;
		else:
			if self.debug:
				print 'DEBUG:SparseAutoencoder:testGradient:Gradient error: ', errorGrad
				print 'DEBUG:SparseAutoencoder:testGradient:Gradient check FAILED!'
				print

			result = -1;

		return result

	def optimizeParameters(self, X):
		'''
		Optimizes the parameters of the Sparse Autoencoder model

		Arguments
		X		: data in the form [input dim., number of samples]

		Returns
		result	: result of the optimization (success or failure)
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:optimizeParameters: The instance is not properly initialized'
		assert X.shape[0]==self.dimLayers[0], 'ERROR:SparseAutoencoder:optimizeParameters: Dimensions of given data do not match with the number of parameters'

		if self.debug: print "DEBUG:SparseAutoencoder:optimizeParameters: Optimizing parameters..."

		# Set optimization options
		method = 'L-BFGS-B'
		options = {};
		options['maxiter'] = 400;

		if self.debug:
			options['disp'] = True;

		# Optimize the cost function
		result = minimize(fun=self.computeCost, jac=self.computeGradient, x0=self.params, args=(X,), method=method, options=options)

		# Set the new values
		self.params = result.x;

		if self.debug: print 'DEBUG:SparseAutoencoder:optimizeParameters: Optimization result: ', result.message

		return result.success;

	def predict(self, X):
		'''
		Applies the Sparse Autoencoder model to the given data

		Arguments
		X		: data in the form [input dim., number of samples]
		W		: weights of the first layer (optional)
		b 		: biases of the first layer (optional)

		Returns
		output of the network in the form [output dim., number of samples]
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:predict: The instance is not properly initialized'
		assert X.shape[0]==self.dimLayers[0], 'ERROR:SparseAutoencoder:predict: Dimensions of given data do not match with the number of parameters'

		[weights, biases] = self.unrollParameters(self.params);

		[outputs, activities] = self.doForwardPropagation(X, weights, biases);

		return activities[0];

	def doForwardPropagation_asLayer(self, X, theta):
		'''
		Wrapper function for doForwardPropagation for cases where Spaerse Autoencoder is
		used as a layer of a deep network.

		Arguments
		X			: data matrix in the form [input dim., number of samples]
		theta		: model parameters for the first layer, must be packed as [weights+biases]

		Returns
		activation	: activation if the first layer
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:doForwardPropagation_asLayer: The instance is not properly initialized'
		assert np.size(theta)==self.weightPrototypes[0][0]*self.weightPrototypes[0][1]+self.biasPrototypes[0][0]*self.biasPrototypes[0][1], 'ERROR:SparseAutoencoder:doForwardPropagation_asLayer: Dimensions of the given model parameter do not match the internal structure'

		theta_local = self.params;
		theta_local[0:np.size(theta)] = theta;
		[weights, biases] = self.unrollParameters(theta_local);

		[outputs, activities] = self.doForwardPropagation(X, weights, biases);

		return activities[0];

	def doBackPropagateError_asLayer(self, error, theta, layer_in, layer_out):
		'''
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:doBackPropagateError_asLayer: The instance is not properly initialized'
		assert np.size(theta)==self.weightPrototypes[0][0]*self.weightPrototypes[0][1]+self.biasPrototypes[0][0]*self.biasPrototypes[0][1], 'ERROR:SparseAutoencoder:doBackPropagateError_asLayer: Dimensions of the given model parameter do not match the internal structure'

		weights = np.reshape(theta[0:self.weightPrototypes[0][0]*self.weightPrototypes[0][1]], self.weightPrototypes[0]);

		delta = error * (layer_out * (1 - layer_out))

		gradients_W = np.dot(delta, np.transpose(layer_in));
		gradients_b = np.sum(delta, 1);

		grad = np.hstack((gradients_W.flatten(), gradients_b.flatten()));

		error_prop = np.dot(np.transpose(weights), delta);

		return grad, error_prop;

	def getParameterSize_asLayer(self):
		'''
		Wrapper function for getParameterSize for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getParameterSize: The instance is not properly initialized'

		return self.weightPrototypes[0][0] * self.weightPrototypes[0][1] + self.biasPrototypes[0][0] * self.biasPrototypes[0][1]

	def getWeights_asLayer(self):
		'''
		Wrapper function for getWeights for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getWeights_asLayer: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);

		return weights[0]

	def getBiases_asLayer(self, layer=0):
		'''
		Wrapper function for getBiases for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getBiases_asLayer: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);

		return biases[0]

	def getParameters_asLayer(self):
		'''
		Wrapper function for getParameters for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:getParameters_asLayer: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);

		return  np.hstack((weights[0].flatten(), biases[0].flatten()))

	def setWeights_asLayer(self, W):
		'''
		Wrapper function for setWeights for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setWeights_asLayer: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);
		weights[0] = W;

		self.params = self.rollParameters(weights, biases);

	def setBiases_asLayer(self, b):
		'''
		Wrapper function for setBiases for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setBiases_asLayer: The instance is not properly initialized'

		[weights, biases] = self.unrollParameters(self.params);
		biases[0] = b;
		self.params = self.rollParameters(weights, biases);

	def setParameters_asLayer(self, theta):
		'''
		Wrapper function for setParameters for cases where Spaerse Autoencoder is
		used as a layer of a deep network. Discards the second layer.
		'''
		assert self.isInitialized, 'ERROR:SparseAutoencoder:setParameters_asLayer: The instance is not properly initialized'
		assert np.size(theta)==np.size(self.params[0:(self.weightPrototypes[0][0]*self.weightPrototypes[0][1]+self.biasPrototypes[0][0]*self.biasPrototypes[0][1])])

		self.params[0:(self.weightPrototypes[0][0]*self.weightPrototypes[0][1]+self.biasPrototypes[0][0]*self.biasPrototypes[0][1])] = theta;

if __name__ == '__main__':

	# Examples:
	# 1) Sparse Autoencoder example on natural images
	# 2) Linear decoder example on color images
	example = 2;

	# Test gradient computation?
	doTest = True;

	if (doTest):
		debug 			= 2;
		numPatches 		= 10;
		patchWidth 		= 2;
		patchHeight 	= 2;
		imChannels		= 1;
		inputDim 		= patchWidth * patchHeight * imChannels;
		numFeatures 	= 5;

		patches = np.random.rand(inputDim, numPatches);

		dimLayers = [inputDim, numFeatures, inputDim];

		SAE_Test = SparseAutoencoder(dimLayers=dimLayers, debug=debug);

		print 'Checking gradient...'

		SAE_Test.testGradient(patches);

	if(example==1):

		print "-------------------"
		print "Autoencoder Example"
		print "-------------------"

		if 1:
		  filename_data = '/home/cem/develop/UFL/data/stlSampledPatches.mat';
		else:
		  filename_data = 'C://develop//python//UFL//data//stlSampledPatches.mat';

		debug 			= 1;
		numPatches 		= 10000;
		patchWidth 		= 8;
		patchHeight 	= 8;
		imChannels		= 1;
		inputDim 		= patchWidth * patchHeight * imChannels;
		numFeatures 	= 25;
		lambda_w 		= 0.0001;     # weight decay parameter
		beta 			= 3;          # weight of sparsity penalty term
		sparsityParam 	= 0.01;
		actFunctions	= [ACTIVATION_FUNCTION_SIGMOID, ACTIVATION_FUNCTION_SIGMOID]

		# Read data from file
		data = scipy.io.loadmat(filename_data);

		patches = DataInputOutput.samplePatches(data['IMAGES'], patchWidth, patchHeight, numPatches);

		# Normalize data
		patches = DataNormalization.normZeroToOne(patches);

		if debug>1:
			Visualization.displayNetwork(patches[:, 0:100]);

		if debug:
			print 'Number of samples: ', patches.shape[1];

		dimLayers = [inputDim, numFeatures, inputDim];

		SAE = SparseAutoencoder(dimLayers=dimLayers, lambda_w=lambda_w, beta=beta, sparsityParam=sparsityParam, actFunctions=actFunctions, debug=debug);

		success = SAE.optimizeParameters(patches);

		# Visualize the learned bases
		weightMatrix = SAE.getWeights();
		Visualization.displayNetwork(np.transpose(weightMatrix));

	elif(example==2):

		print "----------------------"
		print "Linear Decoder Example"
		print "----------------------"

		if 1:
		  filename_data = '/home/cem/develop/UFL/data/stlSampledPatches.mat';
		else:
		  filename_data = 'C://develop//python//UFL//data//stlSampledPatches.mat';

		debug 			= 1;
		numPatches 		= 10000;
		patchWidth 		= 8;
		patchHeight 	= 8;
		imChannels		= 3;
		inputDim 		= patchWidth * patchHeight * imChannels;
		numFeatures 	= 400;
		lambda_w 		= 3e-3;		# weight decay parameter
		beta 			= 5;		# weight of sparsity penalty term
		sparsityParam 	= 0.035;
		actFunctions	= [ACTIVATION_FUNCTION_SIGMOID, ACTIVATION_FUNCTION_IDENTITY]
		epsilon 		= 0.1		# Epsilon for ZCA whitening

		# Read data from file
		data = scipy.io.loadmat(filename_data);
		patches = data['patches'][:, 0:numPatches];

		instance_pca = PCA.PCA(inputDim, epsilon=epsilon, debug=debug);
		ZCAWhite = instance_pca.computeZCAWhiteningMatrix(patches);
		patches_ZCAwhite = instance_pca.doZCAWhitening(patches);

		if debug>1:
			Visualization.displayColorNetwork(patches[:, 0:100]);

		if debug:
			print 'Number of samples: ', patches.shape[1];

		dimLayers = [inputDim, numFeatures, inputDim];

		SAE = SparseAutoencoder(dimLayers=dimLayers, lambda_w=lambda_w, beta=beta, sparsityParam=sparsityParam, actFunctions=actFunctions, debug=debug);

		success = SAE.optimizeParameters(patches_ZCAwhite);

		# Visualize the learned bases
		weightMatrix = SAE.getWeights();
		Visualization.displayColorNetwork(np.transpose(np.dot(weightMatrix, ZCAWhite)));

max_stars_repo_name: cemkaraoguz/UFLDLTutorial

max_issues_repo_name: cemkaraoguz/UFLDLTutorial

max_forks_repo_name: cemkaraoguz/UFLDLTutorial


text:
import numpy as np
import os
import urllib
import gzip
import struct

def download_data(url, force_download=True):
    return url
    fname = url.split("/")[-1]
    if force_download or not os.path.exists(fname):
        urllib.urlretrieve(url, fname)
    return fname

def read_data(label_url, image_url):
    with gzip.open(download_data(label_url)) as flbl:
        magic, num = struct.unpack(">II", flbl.read(8))
        label = np.fromstring(flbl.read(), dtype=np.int8)
    with gzip.open(download_data(image_url), 'rb') as fimg:
        magic, num, rows, cols = struct.unpack(">IIII", fimg.read(16))
        image = np.fromstring(fimg.read(), dtype=np.uint8).reshape(len(label), rows, cols)
    return (label, image)

path='/home/zxteloiv/dataset/mnist/'
(train_lbl, train_img) = read_data(
    path+'train-labels-idx1-ubyte.gz', path+'train-images-idx3-ubyte.gz')
(val_lbl, val_img) = read_data(
    path+'t10k-labels-idx1-ubyte.gz', path+'t10k-images-idx3-ubyte.gz')

ext: py

max_stars_repo_name: zxteloiv/mxnet-ex

max_issues_repo_name: zxteloiv/mxnet-ex

max_forks_repo_name: zxteloiv/mxnet-ex

reason: import numpy


text:
import numpy as np
import pandas as pd
import torch
import mlflow
import mlflow.pyfunc
import warnings
warnings.filterwarnings('ignore')
import pickle
import os
import shutil
import click
import tempfile
from pathlib import Path
# Load model module ... containing model creation and inference methods!
import model

# This will serve as an MLflow wrapper for the model
class ModelWrapper(mlflow.pyfunc.PythonModel):
    # Load in the model and all required artifacts
    # The context object is provided by the MLflow framework
    # It will contain all of the artifacts specified above
    def load_context(self, context):
        import torch
        import pickle

        # Load package meta
        with open(context.artifacts["meta"], 'rb') as handle:
            self.meta = pickle.load(handle)

        self.device = torch.device('cpu')

        # Load the model
        self.model = model.ThousandVoices(
            context.artifacts["model_state"],
            context.artifacts["model_dict"],
            context.artifacts["vocoder_state"],
            context.artifacts["vocoder_config"],
            context.artifacts["xvectors"],
            self.device)

    # Create a predict function for our models
    def predict(self, context, model_input):
        inference_type, args = model_input
        if inference_type == 'TTS':
            sentences, speaker = args
            with torch.no_grad():
                return self.model(sentences, speaker=speaker)

def pack_model(model_file: Path, addition_artifacts: dict, user_meta_data: Path, output_path: Path):
    # MLFlow is very picky about the output directory
    # Generate one temporal folder containing artifacts while generating the
    # mlflow package
    print(f'Writing packaged to {output_path} ...')
    shutil.rmtree(output_path, ignore_errors=True)
    #os.makedirs(f'{output_path}/tmp/', exist_ok=True)

    with tempfile.TemporaryDirectory() as tmp:
        # Serialize the label encoder
        # This will be required at inference time
        #le_path = f'{tmp}/label_encoder.pkl'
        #with open(le_path, 'wb') as handle:
        #    pickle.dump(label_encoder, handle, protocol=pickle.HIGHEST_PROTOCOL)

        # Build package meta dict
        package_meta = {
            'model_file': model_file,
            'user_meta': user_meta_data,
            }

        # Meta data about this package
        meta_path = f'{tmp}/meta.pkl'
        with open(meta_path, 'wb') as handle:
            pickle.dump(package_meta, handle, protocol=pickle.HIGHEST_PROTOCOL)

        # Serialize the models state_dict
        #state_dict_path = f'{tmp}/state_dict.pt'
        #torch.save(model.state_dict(), state_dict_path)

        # Here we will create an artifacts object
        # It will contain all of the data artifacts that we want to package with the model
        artifacts = {
            "meta": meta_path,
        }
        artifacts.update(addition_artifacts)
        env = mlflow.pyfunc.get_default_conda_env()

        # Package the model!
        mlflow.pyfunc.save_model(path=f'{output_path}',
                            python_model=ModelWrapper(),
                            artifacts=artifacts,
                            conda_env=env,
                            code_path=[model_file, user_meta_data])

def _pack(artifacts: Path, package: Path, user_meta: Path, model_file: Path):
    addition_artifacts = {
            "model_state": str(artifacts / "tacotron2/model.best"),
            "model_config": str(artifacts / "tacotron2/model.json"),
            "model_dict":  str(artifacts / "tacotron2/train_clean_460_units.txt"),
            "vocoder_state": str(artifacts / "parallelwavegan/model.pkl"),
            "vocoder_config":  str(artifacts / "parallelwavegan/config.yml"),
            "xvectors": str(artifacts / "xvectors.pkl"),
    }

    print(f'Packing ...')
    pack_model(model_file, addition_artifacts, user_meta, package)

@click.command()
@click.argument('artifacts')
@click.argument('packagepath')
@click.option('--user_meta', default='TTS/meta.json', help='Path to the meta file')
@click.option('--model_file', default='TTS/model.py', help='Path to the python module containing the model')
def pack(artifacts, packagepath, user_meta, model_file):
    _pack(Path(artifacts), Path(packagepath), Path(user_meta), Path(model_file))

if __name__ == '__main__':
    pack()


text:
"""
===============================================================================

Module featuring classes and functionality for unifying events collections.

Classes
-------
EventsUnion

Dependencies
------------
pandas, numpy, copy, warnings, functools

Development
-----------
Developed by:
Tariq Shihadah, tariq.shihadah@gmail.com

Created:
4/13/2022

Modified:
4/13/2022

===============================================================================
"""

################
# DEPENDENCIES #
################

import pandas as pd
import numpy as np
import copy, warnings
from functools import wraps
from rangel import RangeCollection

class EventsUnion(object):
    """
    Parameters
    ----------
    objs : list-like of EventsCollection instances
        A selection of EventsCollection object instances to be combined
        into a single instance based on the input parameters.
    **kwargs
        Keyword arguments to be passed to the initialization function for
        the new EventsCollection instance.
    """

    def __init__(self, objs, **kwargs):
        self.objs = objs

    @property
    def objs(self):
        return self._objs

    @objs.setter
    def objs(self, objs):
        # Validate input object types
        if check_compatibility(objs):
            self._objs = objs

    @property
    def group_keys_unique(self):
        return list(set(
            key for obj in self.objs for key in obj.group_keys_unique))

    def get_groups(self, keys, empty=True):
        """
        Retrieve unique groups of events from each related collection based on
        provided key values.

        Parameters
        ----------
        keys : key value, tuple of key values, or list of the same
            If only one key column is defined within the collections, a single
            column value may be provided. Otherwise, a tuple of column values
            must be provided in the same order as they appear in self.keys. To
            get multiple groups, a list of key values or tuples may be
            provided.
        empty : bool, default True
            Whether to allow for empty events groups to be returned when the
            provided keys are valid but are not associated with any actual
            events. If False, these cases will return a KeyError.
        """
        # Retrieve groups from all collections
        groups = []
        for obj in self._objs:
            group = obj.get_group(keys, empty=empty)
            groups.append(group)
        return groups

    def union(self, fill_gaps=False, suffix=True):
        """
        Combine multiple EventsCollection instances into a single instance,
        creating least common intervals among all collections and maintaining
        all event attributes. The resulting combined events will be used to
        create and return an EventsCollection modeled after the first indexed
        collection in self.objs.

        Parameters
        ----------
        fill_gaps : bool, default False
            Whether to fill gaps in the merged collection with empty events.
            These events would not be associated with any parent collection and
            would not be populated with any events attributes.
        suffix : bool, default True
            Whether to address repeating column labels by appending a suffix
            to all repeating labels indicating the order in which they appear.
            E.g., repeated column labels 'info' and 'info' would be converted
            into 'info_1' and 'info_2'. If False, no modifications will be made
            to column labels which may produce errors when instantiating or
            utilizing the resulting EventsCollection.
        """
        # Iterate over all unique group keys
        records = []
        for keys in self.group_keys_unique:
            # Iterate over groups
            groups = self.get_groups(keys, empty=True)
            ranges = []
            sources = []
            for i, group in enumerate(groups):
                # Create range data for group
                begs = group.df[group.beg].values
                ends = group.df[group.end].values
                rc = RangeCollection(
                    begs=begs, ends=ends, sort=False, copy=False)
                ranges.append(rc)
                # Prepare source data
                source = group.df[self._objs[i].others]
                source = \
                    source.append(pd.Series(dtype=float), ignore_index=True)
                sources.append(source)
            # Union ranges
            rc, indices = RangeCollection.union(
                ranges, fill_gaps=fill_gaps, return_index=True, null_index=-1)
            # Prepare dataframes
            basedata = {col: [val] * rc.num_ranges for col, val \
                in zip(self.objs[0].keys, keys)}
            dfs = [pd.DataFrame({**basedata,
                self.objs[0].beg: rc.begs,
                self.objs[0].end: rc.ends,
            })]
            dfs += [source.iloc[index].reset_index(drop=True) \
                for index, source in zip(indices, sources)]
            records.append(pd.concat(dfs, axis=1))
        # Concatenate records
        df = pd.concat(records)
        # Address repeating column labels
        if suffix:
            cols = list(df)
            counts = [cols.count(col) for col in cols]
            suffixes = ['_' + str(cols[:i+1].count(col)) \
                for i, col in enumerate(cols)]
            labels = [col + (suffix if count > 1 else '') \
                for col, suffix, count in zip(cols, suffixes, counts)]
            df.columns = labels
        # Create events collection based on objs[0]
        ec = self.objs[0].from_similar(df, geom=None)
        return ec

#####################
# LATE DEPENDENCIES #
#####################

from linref.events.collection import EventsCollection, check_compatibility

max_stars_repo_name: tariqshihadah/linref

max_issues_repo_name: tariqshihadah/linref

max_forks_repo_name: tariqshihadah/linref


text:
# %%
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#v = np.array([int(x) for x in '16,1,2,0,4,2,7,1,2,14'.split(',')])
v = np.array([int(x) for x in Path('data.txt').read_text().split(',')], dtype='int32')

v_set = set(v)

min_pos = -1
min_val = 100000000

for x in v_set:
	d = np.sum(np.abs(v - x))
	if d < min_val:
		min_val = d
		min_pos = x

print('# of elements:', len(v))
print('Min sum:', min_val)
print('pos:', min_pos)
# %%
# Problem B
def distance(d):
	return d*(d+1)/2

print('min:', np.min(v), 'max:', np.max(v))
q = np.zeros(np.max(v)+1, dtype='int32')
for x in range(np.max(v)+1):
	q[x] = np.sum(distance(np.abs(v - x)))
plt.plot(q)
print(min(q))


text:
import numpy as np
import h5py
from scipy import stats, mgrid, c_, reshape
from dateutil    import parser

from prime_utils import runningAvg, prediction_filename

def getKDE(spl,nskip=0,nthin=1,npts=100,bwfac=1.0):
    r"""
    Compute 1D and 2D marginal PDFs via Kernel Density Estimate

    Parameters
    ----------
    spl: numpy array
         MCMC chain [number of samples x number of parameters]
    nskip: int
           number of initial samples to skip when sampling the MCMC chain
    nthin: int
           use every 'nthin' samples
    npts: int
          number of grid points
    bwfac: double
          bandwidth factor

    Returns
    -------
    dict: dictionary with results
         'x1D': list of numpy arrays with grids for the 1D PDFs;
         'p1D': list of numpy arrays with 1D PDFs;
         'x2D': list of numpy arrays of x-axis grids for the 2D PDFs;
         'y2D': list of numpy arrays of y-axis grids for the 2D PDFs;
         'p2D': list of numpy arrays containing 2D PDFs
    """
    spl=spl[nskip::nthin,:]
    # Compute 1D pdf's
    x1D=[];
    p1D=[];
    for i in range(spl.shape[1]):
        print(" - Evaluate 1D pdf for variable",i+1)
        spls = spl[:,i]
        x1D.append(np.linspace(spls.min(),spls.max(),npts));
        kern=stats.kde.gaussian_kde(spls);
        p1D.append(kern(x1D[i]));
    # Compute 2D joint pdf's
    npts = 1j*npts
    x2D = []
    y2D = []
    p2D = []
    for j in range(1,spl.shape[1]):
        for i in range(j):
            print(" - Evaluate 2D pdf for variables",i+1,j+1)
            kern=stats.kde.gaussian_kde(c_[spl[:,i],spl[:,j]].T)
            kern.set_bandwidth(bw_method=kern.factor*bwfac)
            xmin = spl[:,i].min(); xmax = spl[:,i].max()
            ymin = spl[:,j].min(); ymax = spl[:,j].max()
            x,y = mgrid[xmin:xmax:npts,ymin:ymax:npts]
            z   = reshape(kern(c_[x.ravel(), y.ravel()].T).T, x.T.shape)
            x2D.append(x);
            y2D.append(y);
            p2D.append(z);
    return {'x1D':x1D,'p1D':p1D,'x2D':x2D,'y2D':y2D,'p2D':p2D}

def distcorr(spl):
    r"""
    Compute distance correlation between random vectors

    Parameters
    ----------

    spl: numpy array [number of samples x number of variables]
         first dimension is the number of samples,
         second dimension is the number of random vectors

    Returns
    -------
       Returns a 2D array of distance correlations between pairs of random vectors;
       only entries 0<=j<i<no. of random vectors are populated

    References:
       http://en.wikipedia.org/wiki/Distance_correlation
     """
    nspl,nvars = spl.shape
    if (nspl>5000):
        print('Warning ! This might be a lengthy calculation: nspl=',nspl)
    As=[]
    for i in range(nvars):
        Amat = np.matrix(np.zeros((nspl,nspl)))
        for i1 in range(nspl):
            for j1 in range(nspl):
                Amat[i1,j1] = abs(spl[i1,i]-spl[j1,i])
        # compute means
        Alin = np.array([np.mean(Amat[i1,:]) for i1 in range(nspl)])
        Acol = np.array([np.mean(Amat[:,j1]) for j1 in range(nspl)])
        Amn  = np.mean(Amat)
        #print Amat
        #for i in range(nspl):
        #    print i,Alin[i],Acol[i]
        # subtract/add means (linewise, columnwise, overall)
        Amat = Amat - Alin.reshape(nspl,1)
        Amat = (Amat.T - Acol.reshape(nspl,1)).T
        Amat = Amat+Amn
        #print Amat
        As.append(Amat.copy())
    dCor = np.zeros((nvars,nvars))
    dVarX = [np.sqrt(np.sum(np.multiply(As[i],As[i]))/(nspl*nspl))
       for i in range(nvars)]
    print("Variances:")
    print(dVarX)
    for i in range(1,nvars):
        for j in range(i):
            dCov   = np.sqrt(np.sum(np.multiply(As[i],As[j]))/(nspl*nspl))
            dCor[i,j] = dCov/np.sqrt(dVarX[i]*dVarX[j])
    As=[]
    return dCor

def computeAICandBIC(run_setup,verbose=0):
    """
    Compute Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC)

    Parameters
    ----------
    run_setup: dictionary with run settings; see the Examples section in the manual

    Returns
    -------
    AIC: float
    BIC: float
    """
    #-------------------------------------------------------
    # definitions
    fdata = run_setup["regioninfo"]["regionname"]+".dat"
    fchno = run_setup["regioninfo"]["fchain"]

    #-------------------------------------------------------
    # retrieve log likelihoods
    file = h5py.File(fchno, 'r')
    loglik = np.array(file["minfo"][:,1])
    file.close()

    if verbose>0:
        print(loglik.shape)

    #-------------------------------------------------------
    # compute AIC and BIC

    # Read in initial parameter guess to obtain number of parameters
    Nparam = len(run_setup["mcmcopts"]["cini"])

    # compute number of data points used (number of days)
    rawdata = np.loadtxt(fdata,dtype=str)
    ndays = rawdata.shape[0]

    AIC = 2 * Nparam - 2 * np.max(loglik)
    BIC = 2 * Nparam * np.log(ndays) - 2 * np.max(loglik)

    return AIC, BIC

def computeCRPS(run_setup):
    """
    Compute Continuous Rank Predictive Score (CRPS)

    Parameters
    ----------
    run_setup: dictionary with run settings; see the Examples section in the manual

    Returns
    -------
    CRPS: float
    """
    #-------------------------------------------------------
    # definitions
    fdata = run_setup["regioninfo"]["regionname"]+".dat"
    fchno = run_setup["regioninfo"]["fchain"]
    day0  = run_setup["regioninfo"]["day0"]

    #-------------------------------------------------------
    # extract data from raw data
    rawdata = np.loadtxt(fdata,dtype=str)
    ndays_data = rawdata.shape[0]
    days_since_day0 = np.array([(parser.parse(rawdata[i,0])-parser.parse(day0)).days for i in range(ndays_data)])
    new_cases       = runningAvg(np.array([float(rawdata[i,1]) for i in range(ndays_data)]),7)

    #------------------------------------------------------
    # read model predictions
    filename = prediction_filename(run_setup)
    file   = h5py.File(filename, 'r')
    pred   = np.array(file["predictions"])
    file.close()

    # only need predictions for days on which data is available
    pred = pred[:,:ndays_data]

    #------------------------------------------------------
    # compute integrals on each day
    CRPS_daily = np.zeros(ndays_data)

    for i in range(ndays_data):
        pred_i = pred[:,i]

        # compute empirical pdf and cdf
        y = np.sort(pred_i)
        ecdf = np.array(range(1,len(pred_i)+1)) / float(len(pred_i))

        # compute heaviside function
        h_i = np.ones(ecdf.shape) * (y > new_cases[i])

        # compute integrand
        integrand = (ecdf - h_i) ** 2

        # average for trapezoid rule integration
        integrand = (integrand[1:] + integrand[:-1]) / 2

        # compute inegral for day i
        CRPS_daily[i] = np.sum( (y[1:] - y[:-1]) * integrand)

    CRPS = np.sum(CRPS_daily) / ndays_data
    return CRPS

max_stars_repo_name: sandialabs/PRIME

max_issues_repo_name: sandialabs/PRIME

max_forks_repo_name: sandialabs/PRIME

reason: import numpy,from scipy


text:
# coding: utf-8
import numpy as np
from scipy.sparse import coo_matrix

from .generic_utils import handleKeyError

def take_centers(a):
    """Take a center values.

    Args:
        a (array) : Array

    Example:
        >>> from pycharmers.utils import take_centers
        >>> take_centers([0,1,2,3,4,5])
        array([0.5, 1.5, 2.5, 3.5, 4.5])
        >>> take_centers([[0,1,2,3,4,5],[1,2,3,4,5,6]])
        array([[0.5, 1.5, 2.5, 3.5, 4.5],
            [1.5, 2.5, 3.5, 4.5, 5.5]])
    """
    a = np.asarray(a, dtype=np.float64)
    return (a[Ellipsis,1:] + a[Ellipsis,:-1])/2

def confusion_matrix(y_true, y_pred):
    """Compute confusion matrix to evaluate the accuracy of a classification.

    By definition a confusion matrix :math:`C` is such that :math:`C_{i, j}` is
    equal to the number of observations known to be in group :math:`i` and predicted
    to be in group :math:`j`.

    Thus in binary classification, the count of true negatives is
    :math:`C_{0,0}`, false negatives is :math:`C_{1,0}`, true positives is
    :math:`C_{1,1}` and false positives is :math:`C_{0,1}`.

    Args:
        y_true (array): Ground truth (correct) target values.
        y_pred (array): Estimated targets as returned by a classifier.

    Returns:
        cm (array): Confusion matrix whose i-th row and j-th column entry indicates
                    the number of samples with true label being i-th class and prediced label being j-th class.

    References:
        `Wikipedia entry for the Confusion matrix <https://en.wikipedia.org/wiki/Confusion_matrix>`_ (Wikipedia and other references may use a different convention for axes)

    Examples:
        >>> from pycharmers import confusion_matrix
        >>> y_true = [2, 0, 2, 2, 0, 1]
        >>> y_pred = [0, 0, 2, 2, 0, 2]
        >>> confusion_matrix(y_true, y_pred)
        array([[2, 0, 0],
            [0, 0, 1],
            [1, 0, 2]])
        >>> # In the binary case, we can extract true positives, etc as follows:
        >>> tn, fp, fn, tp = confusion_matrix([0, 1, 0, 1], [1, 1, 1, 0]).ravel()
        >>> (tn, fp, fn, tp)
        (0, 2, 1, 1)
    """
    y_true = np.asarray(y_true).flatten(order='C')
    y_test = np.asarray(y_pred).flatten(order='C')
    num_labels = len(np.unique(np.concatenate([y_true, y_pred])))
    num_data   = len(y_true)
    cm = coo_matrix(
        (np.ones(num_data, dtype=np.int64), (y_true, y_pred)),
        shape=(num_labels, num_labels), dtype=np.int64,
    ).toarray()
    return cm

def rotate2d(a, theta):
    """Rotate 2d vectors using Rotation matrix :math:`R(\\theta)`

    .. math::

        R(\\theta) = \\left( \\begin{array}{c} \\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{array} \\right)

    Using this matrix, the rotation of the vector can be expressed as follows.

    .. math::

        \\left( \\begin{array}{c} x^{\\prime} \\\\ y^{\\prime} \\end{array} \\right) = \\left( \\begin{array}{c} \\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{array} \\right)\\left( \\begin{array}{c} x \\\\ y \\end{array} \\right)

    Args:
        a (ndarray)   : Array
        theta (float) : float value like `1/2*np.pi`
    """
    c, s = np.cos(theta), np.sin(theta)
    R = np.array([[c, -s],[s, c]])
    return (R @ a.reshape(2,1)).squeeze()

def replaceArray(a, old=(255,255,255), new=(0,0,0)):
    """Replace an Array from ``old`` to ``new``

    Args:
        old (tuple) : Old value.
        new (tuple) : New value.

    Returns:
        np.ndarray: New Array

    Examples:
        >>> import cv2
        >>> import numpy as np
        >>> from pycharmers.opencv import SAMPLE_LENA_IMG, cv2plot
        >>> from pycharmers.utils import replaceArray
        >>> img = cv2.imread(SAMPLE_LENA_IMG)
        >>> img = replaceArray(img, old=[77, 66, 176], new=[0,0,0]).astype(np.uint8)
        >>> cv2plot(img, is_cv2=True)
    """
    ch = a.shape[2]
    if not hasattr(old, "__len__"):
        old = [old]*ch
    if not hasattr(new, "__len__"):
        new = [new]*ch
    return np.where(np.expand_dims(np.all(a==old, axis=-1), axis=-1), new, a)

def fill_between_angle(arr, s, e, center=None, is_radian=True):
    """Fill Between ``s`` and ``e``.

    Args:
        arr (np.ndarray)          : Input array.
        s (Number)                : Start angle of fill.
        e (Number)                : End angle of fill.
        center (tuple)            : Center coordinates.
        is_radian (bool, optional): whether ``s`` and ``e`` are defined in radians.

    Returns:
        np.ndarray: Whether it is a place to be filled.

    Examples:
        >>> import numpy as np
        >>> from PIL import Image
        >>> from pycharmers.utils import fill_between_angle
        >>> arr = np.zeros(shape=(100,100,3)).astype(np.uint8)
        >>> flag = fill_between_angle(arr, s=30, e=120, is_radian=False)
        >>> Image.fromarray(np.where(flag, arr, 255))
    """
    H,W = arr.shape[:2]
    flag = np.zeros_like(arr, dtype=np.bool)
    if not is_radian:
        s = np.pi*(s/180)
        e = np.pi*(e/180)
    if center is None:
        center = (H//2,W//2)
    cx,cy = center
    for i in range(H):
        ry = H-i-cy
        for j in range(W):
            rx = j-cx
            radi = np.arctan2(ry, rx)
            if ry<0:
                radi+=2*np.pi
            flag[i,j] = s<=radi<=e
    return flag

max_stars_repo_name: iwasakishuto/Py-utils

max_issues_repo_name: iwasakishuto/Py-utils

max_forks_repo_name: iwasakishuto/Py-utils


text:
import fatd.holders
import numpy as np

def csv_loader(csv_path, header=None):
    data_matrix = np.loadtxt(csv_path, delimiter=',')

    if header is None:
        data = fatd.holders.Data(data_matrix[:,:-1],
                                 data_matrix[:,-1]
                                )
    else:
        data = fatd.holders.Data(data_matrix[:,:-1],
                                 data_matrix[:,-1],
                                 header[:-1],
                                 header[-1]
                                )

    return data

max_stars_repo_name: AnthropocentricAI/fat-dummy

max_issues_repo_name: AnthropocentricAI/fat-dummy

max_forks_repo_name: AnthropocentricAI/fat-dummy


text:
import numpy as np

def count_true_positive(truth, recommend):
    """Count number of true positives from given sets of samples.

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.

    Returns:
        int: Number of true positives.

    """
    tp = 0
    for r in recommend:
        if r in truth:
            tp += 1
    return tp

def recall(truth, recommend, k=None):
    """Recall@k.

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.
        k (int): Top-k items in `recommend` will be recommended.

    Returns:
        float: Recall@k.

    """
    if len(truth) == 0:
        if len(recommend) == 0:
            return 1.
        return 0.

    if k is None:
        k = len(recommend)
    return count_true_positive(truth, recommend[:k]) / float(truth.size)

def precision(truth, recommend, k=None):
    """Precision@k.

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.
        k (int): Top-k items in `recommend` will be recommended.

    Returns:
        float: Precision@k.

    """
    if len(recommend) == 0:
        if len(truth) == 0:
            return 1.
        return 0.

    if k is None:
        k = len(recommend)
    return count_true_positive(truth, recommend[:k]) / float(k)

def average_precision(truth, recommend):
    """Average Precision (AP).

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.

    Returns:
        float: AP.

    """
    if len(truth) == 0:
        if len(recommend) == 0:
            return 1.
        return 0.

    tp = accum = 0.
    for n in range(recommend.size):
        if recommend[n] in truth:
            tp += 1.
            accum += (tp / (n + 1.))
    return accum / truth.size

def auc(truth, recommend):
    """Area under the ROC curve (AUC).

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.

    Returns:
        float: AUC.

    """
    tp = correct = 0.
    for r in recommend:
        if r in truth:
            # keep track number of true positives placed before
            tp += 1.
        else:
            correct += tp
    # number of all possible tp-fp pairs
    pairs = tp * (recommend.size - tp)

    # if there is no TP (or no FP), it's meaningless for this metric (i.e., AUC=0.5)
    if pairs == 0:
        return 0.5

    return correct / pairs

def reciprocal_rank(truth, recommend):
    """Reciprocal Rank (RR).

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.

    Returns:
        float: RR.

    """
    for n in range(recommend.size):
        if recommend[n] in truth:
            return 1. / (n + 1)
    return 0.

def mpr(truth, recommend):
    """Mean Percentile Rank (MPR).

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.

    Returns:
        float: MPR.

    """
    if len(recommend) == 0 and len(truth) == 0:
        return 0.  # best
    elif len(truth) == 0 or len(truth) == 0:
        return 100.  # worst

    accum = 0.
    n_recommend = recommend.size
    for t in truth:
        r = np.where(recommend == t)[0][0] / float(n_recommend)
        accum += r
    return accum * 100. / truth.size

def ndcg(truth, recommend, k=None):
    """Normalized Discounted Cumulative Grain (NDCG).

    Args:
        truth (numpy 1d array): Set of truth samples.
        recommend (numpy 1d array): Ordered set of recommended samples.
        k (int): Top-k items in `recommend` will be recommended.

    Returns:
        float: NDCG.

    """
    if k is None:
        k = len(recommend)

    def idcg(n_possible_truth):
        res = 0.
        for n in range(n_possible_truth):
            res += 1. / np.log2(n + 2)
        return res

    dcg = 0.
    for n, r in enumerate(recommend[:k]):
        if r not in truth:
            continue
        dcg += 1. / np.log2(n + 2)

    res_idcg = idcg(np.min([truth.size, k]))
    if res_idcg == 0.:
        return 0.
    return dcg / res_idcg

max_stars_repo_name: KevinGetandGive/flurs

max_issues_repo_name: KevinGetandGive/flurs

max_forks_repo_name: KevinGetandGive/flurs


text:
import numpy as np
import astropy.io.fits as pf
from configobj import ConfigObj
import ipdb
import matplotlib.pyplot as plt
import medis.speckle_nulling.sn_hardware as hardware
from validate import Validator
import medis.speckle_nulling.sn_preprocessing as pre
import medis.speckle_nulling.sn_processing as pro
from scipy.interpolate import interp1d
import medis.speckle_nulling.sn_math as snm
import dm_registration as dm
import PID as pid

def fit_satellite_centers(image, spotcenters_init, window=20):
    """ Returns the center of the image given approximative spot centers """
    # fit satellite centers
    scs = np.zeros((len(spotcenters_init), 2))
    for idx,xy in enumerate(spotcenters_init):
        subim = pre.subimage(image, xy, window=window)
        popt = snm.image_centroid_gaussian1(subim)
        xcenp = popt[1]
        ycenp = popt[2]
        xcen = xy[0]-round(window/2)+xcenp
        ycen = xy[1]-round(window/2)+ycenp
        scs[idx,:] = xcen, ycen

    return scs

def get_delta_I(image, cx = None, cy = None,
                quad_width_pix = 7., inner_rad_pix = 0.,
                zone_type = None) :

    if zone_type == 'inner' :
        image = image * pro.circle(image, cx, cy, inner_rad_pix)
    elif zone_type == 'outer' :
        image = image * (1 - pro.circle(image, cx, cy, inner_rad_pix))

    window = quad_width_pix*2.
    xs, ys = np.meshgrid(np.arange(image.shape[0]), np.arange(image.shape[1]))
    subim = pre.subimage(image, (round(cx), round(cy)), window=window)
    subx = pre.subimage(xs, (round(cx), round(cy)), window=window)
    suby = pre.subimage(ys, (round(cx), round(cy)), window=window)

    cumx = np.cumsum(np.sum(subim, axis = 1))
    cumy = np.cumsum(np.sum(subim, axis = 0))

    f_interpx = interp1d(subx[0,:], cumx)
    f_interpy = interp1d(suby[:,0], cumy)

    deltaIx = max(cumx)-2*f_interpx(cx-1)
    deltaIy = max(cumy)-2*f_interpy(cy-1)

    #ipdb.set_trace()

    #plt.imshow(subim)
    #plt.show()

    #plt.plot(cumx)
    #plt.show()

    return deltaIy, deltaIx

def tiptiltestimator_circle(image, cx=None, cy= None, beta = None,
                     lambdaoverd=None, window = None):
    def cuberoot(x):
        if x<0:
            return -(-x)**(1.0/3)
        else:
            return x**(1.0/3)

    xs, ys = np.meshgrid(np.arange(image.shape[0]),
                         np.arange(image.shape[1]))
    subim = pre.subimage(image, (round(cx), round(cy)), window=window)
    subx = pre.subimage(xs, (round(cx), round(cy)), window=window)
    suby = pre.subimage(ys, (round(cx), round(cy)), window=window)

    cumx = np.cumsum(np.sum(subim, axis = 1))
    cumy = np.cumsum(np.sum(subim, axis = 0))

    f_interpx = interp1d(subx[0], cumx)
    f_interpy = interp1d(suby[:,0], cumy)

    deltaIx = max(cumx)-2*f_interpx(cx)
    deltaIy = max(cumy)-2*f_interpy(cy)

    Tx = cuberoot(deltaIx/beta)*cuberoot(deltaIx**2/(deltaIx**2+deltaIy**2))
    Ty = cuberoot(deltaIy/beta)*cuberoot(deltaIy**2/(deltaIx**2+deltaIy**2))
    return Tx, Ty

def tiptiltestimator(delat_i_x, delta_i_y, gamma = 0., rotangle = 0.):

    # tip-tilt estimation
    theta = np.arctan(delta_i_y/delta_i_x)
    delta_i_theta = np.sqrt(delta_i_x**2.+delta_i_y**2.)
    T_theta = delta_i_theta/gamma

    Tx = T_theta * np.cos(theta-rotangle)
    Ty = T_theta * np.sin(theta-rotangle)

    return Tx, Ty

if __name__ == "__main__":
    configfilename = 'qacits_config.ini'
    hardwareconfigfile = 'speckle_instruments.ini'
    configspecfile = 'qacits_config.spec'
    config = ConfigObj(configfilename, configspec=configspecfile)
    val = Validator()
    check = config.validate(val)

    pharo = hardware.fake_pharo()
    #p3k = hardware.P3K_COM('P3K', configfile = hardwareconfigfile)

    bgd = np.zeros((1024, 1024))
    flat = np.ones((1024, 1024))
    bpix = np.zeros((1024, 1024))

    # parameters defining the zone of interest
    spotcenters = np.resize(config['Image_params']['spotcenters_init'], (4,2))
    quad_width_pix = config['Image_params']['quad_width'] * config['Image_params']['lambdaoverd']
    inner_rad_pix = config['Image_params']['inner_rad'] * config['Image_params']['lambdaoverd']
    lambdaoverd_arc = config['Image_params']['lambdaoverd_arc']
    # reference values
    Itot_off = config['QACITS_params']['Itot_off']
    DIx_ref = config['QACITS_params']['DIx_ref']
    DIy_ref = config['QACITS_params']['DIy_ref']

    # PID loop gains
    Kp = config['PID']['Kp']
    Ki = config['PID']['Ki']
    Kd = config['PID']['Kd']

    #ipdb.set_trace()

    # PID loop
    p = pid.PID(P=np.array([Kp,Kp]),
                I=np.array([Ki,Ki]),
                D=np.array([Kd,Kd]),
                Deadband=.001)
    p.setPoint(np.array[0.,0.])

    c=1
    while True:
        img = pharo.get_image()
        img = pre.equalize_image(img,
                                 masterflat = flat,
                                 bkgd = bgd,
                                 badpix = bpix)

        # Derive center of the image from the satellite spots
#        if c == 1 :
#            spotcenters = dm.get_satellite_centroids(img)
#        else :
#            spotcenters = fit_satellite_centers(img, spotcenters, window=20)

        spotcenters = fit_satellite_centers(img, spotcenters, window=20)

        centerx, centery = np.mean(spotcenters, axis = 0)
        #centerx, centery = 512, 512

        print centerx, centery

        if c== 3 : ipdb.set_trace()

        delta_i_x, delta_i_y = get_delta_I(img, cx = centerx, cy=centery,
                                quad_width_pix = quad_width_pix,
                                inner_rad_pix = inner_rad_pix,
                                zone_type = "inner")

        delta_i_x = (delta_i_x - DIx_ref) / Itot_off
        delta_i_y = (delta_i_y - DIy_ref) / Itot_off

        # tip tilt estimator in lambda over D
        Tx, Ty = tiptiltestimator(delta_i_x, delta_i_y, gamma = 1.)
        # conversion in arcsec to feed the P3K tip tilt mirror
        rx = Tx*lambdaoverd_arc
        ry = Ty*lambdaoverd_arc

        offx, offy = snm.rotateXY(rx, ry, config['AO']['rotang'])

        # command value according to PID loop
        p.update([offx,offy])

        print offx, offy
        print '-----------'

        c=c+1
        #offx and offy that you feed to p3k needs to be in arc seconds
        #p3k.sci_offset_up(offy)
        #p3k.sci_offset_left(offx)

max_stars_repo_name: RupertDodkins/medis

max_issues_repo_name: RupertDodkins/medis

max_forks_repo_name: RupertDodkins/medis

reason: import numpy,from scipy,import astropy


text:
'''
#                           This is an 80 character line                       #
Compute the length of the cluster edge:
-Use Freud to find the complete system neighborlist
-Grab the largest cluster
-Mesh the system
-Compute which bins have largest cluster particles
-If adjacent bins are empty, the reference bin is an edge
-Multiply by bin size to get length
'''

import sys

# Run locally
sys.path.append('/Users/kolbt/Desktop/compiled/hoomd-blue/build')
sys.path.append('/Users/kolbt/Desktop/compiled/gsd/build')
# Run on the cpu
sys.path.append('/nas/longleaf/home/kolbt/programs/cpu-hoomd/hoomd-blue/build')
# Run on the gpu

sys.path.append('/nas/longleaf/home/kolbt/programs/hoomd_2.2.1/hoomd-blue/build')
sys.path.append('/nas/longleaf/home/kolbt/programs/gsd/build')

import gsd
from gsd import hoomd
from gsd import pygsd

import freud
from freud import parallel
from freud import box
from freud import density
from freud import cluster

import numpy as np
import math
import random
from scipy import stats

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.collections
from matplotlib import colors
import matplotlib.gridspec as gridspec
import matplotlib.patches as patches

def computeDist(x1, y1, x2, y2):
    '''Compute distance between two points'''
    return np.sqrt( ((x2-x1)**2) + ((y2 - y1)**2) )

def computeFLJ(r, x1, y1, x2, y2, eps):
    sig = 1.
    f = (24. * eps / r) * ( (2*((sig/r)**12)) - ((sig/r)**6) )
    fx = f * (x2 - x1) / r
    fy = f * (y2 - y1) / r
    return fx, fy

def computeTauPerTstep(epsilon, mindt=0.000001):
    '''Read in epsilon, output tauBrownian per timestep'''
#    if epsilon != 1.:
#        mindt=0.00001
    kBT = 1.0
    tstepPerTau = float(epsilon / (kBT * mindt))
    return 1. / tstepPerTau

def roundUp(n, decimals=0):
    '''Round up size of bins to account for floating point inaccuracy'''
    multiplier = 10 ** decimals
    return math.ceil(n * multiplier) / multiplier

def getNBins(length, minSz=(2**(1./6.))):
    "Given box size, return number of bins"
    initGuess = int(length) + 1
    nBins = initGuess
    # This loop only exits on function return
    while True:
        if length / nBins > minSz:
            return nBins
        else:
            nBins -= 1

def findBins(lookN, currentInd, maxInds):
    '''Get the surrounding bin indices'''
    maxInds -= 1
    left = currentInd - lookN
    right = currentInd + lookN
    binsList = []
    for i in range(left, right):
        ind = i
        if i > maxInds:
            ind -= maxInds
        binsList.append(ind)
    return binsList

# Get infile and open
inFile = str(sys.argv[1])
if inFile[0:7] == "cluster":
    add = 'cluster_'
else:
    add = ''

f = hoomd.open(name=inFile, mode='rb')
# Inside and outside activity from command line
peA = float(sys.argv[2])
peB = float(sys.argv[3])
parFrac = float(sys.argv[4])
eps = float(sys.argv[5])
try:
    phi = float(sys.argv[6])
    intPhi = int(phi)
    phi /= 100.
except:
    phi = 0.6
    intPhi = 60
try:
    dtau = float(sys.argv[7])
except:
    dtau = 0.000001

# Get number of particles for textfile name
with hoomd.open(name=inFile, mode='rb') as t:
    snap = t[0]
    typ = snap.particles.typeid
    partNum = len(typ)

# Outfile to write data to
base = add + 'pressure_pa' + str(peA) +\
       '_pb' + str(peB) +\
       '_xa' + str(parFrac) +\
       '_phi' + str(intPhi) +\
       '_ep' + '{0:.3f}'.format(eps) +\
       '_N' + str(partNum)
outFile = base + '.txt'
imgFile = base + '.png'

g = open(outFile, 'w') # write file headings
g.write('Timestep'.center(10) + ' ' +\
        'gasArea'.center(20) + ' ' +\
        'gasTrace'.center(20) + ' ' +\
        'gasPress'.center(20) + ' ' +\
        'bulkArea'.center(20) + ' ' +\
        'bulkTrace'.center(20) + ' ' +\
        'bulkPress'.center(20) + ' ' +\
        'SurfaceTense'.center(20) + ' ' +\
        'Length'.center(20) + ' ' +\
        'NDense'.center(20) + '\n')
g.close()

start = 0                   # first frame to process
dumps = int(f.__len__())    # get number of timesteps dumped
end = dumps                 # final frame to process
#start = end - 1

box_data = np.zeros((1), dtype=np.ndarray)  # box dimension holder
r_cut = 2**(1./6.)                          # potential cutoff
tauPerDT = computeTauPerTstep(epsilon=eps)  # brownian time per timestep

with hoomd.open(name=inFile, mode='rb') as t:
    snap = t[0]
    first_tstep = snap.configuration.step
    box_data = snap.configuration.box
    l_box = box_data[0]
    h_box = l_box / 2.
    typ = snap.particles.typeid
    partNum = len(typ)
    # Set up cluster computation using box
    f_box = box.Box(Lx=l_box, Ly=l_box, is2D=True)
    my_clust = cluster.Cluster()
    c_props = cluster.ClusterProperties()
    # Compute each mesh
    NBins = getNBins(l_box, r_cut)
    sizeBin = roundUp((l_box / NBins), 6)

    # Loop through each timestep
    for j in range(start, end):
        snap = t[j]
        # Easier accessors
        pos = snap.particles.position               # position
        pos[:,-1] = 0.0
        xy = np.delete(pos, 2, 1)
        typ = snap.particles.typeid                 # type
        tst = snap.configuration.step               # timestep
        tst -= first_tstep                          # normalize by first timestep
        tst *= dtau                                 # convert to Brownian time

        # Compute clusters for this timestep
        system = freud.AABBQuery(f_box, f_box.wrap(pos))

        # Compute neighbor list for only largest cluster
        my_clust.compute(system, neighbors={'r_max': 1.0})
        ids = my_clust.cluster_idx              # get id of each cluster
        c_props.compute(system, ids)            # find cluster properties
        clust_size = c_props.sizes              # find cluster sizes

        # We can also grab all clusters over a set size
        min_size = 5000

        # Get the positions of all particles in LC
        binParts = [[[] for b in range(NBins)] for a in range(NBins)]
        occParts = [[0 for b in range(NBins)] for a in range(NBins)]
        edgeBin = [[0 for b in range(NBins)] for a in range(NBins)]
        liqPos = []
        gasPos = []
        for k in range(0, len(ids)):
            # Convert position to be > 0 to place in list mesh
            tmp_posX = pos[k][0] + h_box
            tmp_posY = pos[k][1] + h_box
            x_ind = int(tmp_posX / sizeBin)
            y_ind = int(tmp_posY / sizeBin)
            # Append all particles to appropriate bin
            binParts[x_ind][y_ind].append(k)

            # Get sufficient cluster mesh as well
            if clust_size[ids[k]] >= min_size:
                liqPos.append(pos[k])
                occParts[x_ind][y_ind] = 1
            # Get a gas particle list as well
            elif clust_size[ids[k]] <= 100:
                gasPos.append(pos[k])

        # If sufficient neighbor bins are empty, we have an edge
        thresh = 1.5
        # Loop through x index of mesh
        for ix in range(0, len(occParts)):

            # If at right edge, wrap to left
            if (ix + 1) != NBins:
                lookx = [ix-1, ix, ix+1]
            else:
                lookx = [ix-1, ix, 0]

            # Loop through y index of mesh
            for iy in range(0, len(occParts[ix])):

                # Reset neighbor counter
                count = 0
                # If the bin is not occupied, skip it
                if occParts[ix][iy] == 0:
                    continue
                # If at top edge, wrap to bottom
                if (iy + 1) != NBins:
                    looky = [iy-1, iy, iy+1]
                else:
                    looky = [iy-1, iy, 0]

                # Loop through surrounding x-index
                for indx in lookx:
                    # Loop through surrounding y-index
                    for indy in looky:

                        # If neighbor bin is NOT occupied
                        if occParts[indx][indy] == 0:
                            # If neighbor bin shares a vertex
                            if indx != ix and indy != iy:
                                count += 0.5
                            # If neighbor bin shares a side
                            else:
                                count += 1

                # If sufficient neighbors are empty, we found an edge
                if count >= thresh:
                    edgeBin[indx][indy] = 1

        blurBin = [[0 for b in range(NBins)] for a in range(NBins)]
        Nedges = 0
        # Blur the edge bins a bit
        for ix in range(0, len(occParts)):
            for iy in range(0, len(occParts)):
                Nedges += edgeBin[ix][iy]
                if edgeBin[ix][iy] == 1:
                    # If at right edge, wrap to left
                    if (ix + 1) != NBins:
                        lookx = [ix-1, ix, ix+1]
                    else:
                        lookx = [ix-1, ix, 0]
                    # If at top edge, wrap to bottom
                    if (iy + 1) != NBins:
                        looky = [iy-1, iy, iy+1]
                    else:
                        looky = [iy-1, iy, 0]
                    # Loop through surrounding x-index
                    for indx in lookx:
                        # Loop through surrounding y-index
                        for indy in looky:
                            # Make all surrounding bins 'edge' bins
                            blurBin[indx][indy] = 1

        # Now let's compute the pressure
        bulkSigXX = 0
        bulkSigYY = 0
        gasSigXX = 0
        gasSigYY = 0
        print("Computing the pressure")
        for k in range(0, len(ids)):
            # Convert position to be > 0 to place in list mesh
            tmp_posX = pos[k][0] + h_box
            tmp_posY = pos[k][1] + h_box
            x_ind = int(tmp_posX / sizeBin)
            y_ind = int(tmp_posY / sizeBin)
            # Only compute pressure for non-edge bins
            if blurBin[x_ind][y_ind] != 1:
                # If at right edge, wrap to left
                if (x_ind + 1) != NBins:
                    lookx = [x_ind-1, x_ind, x_ind+1]
                else:
                    lookx = [x_ind-1, x_ind, 0]
                # If at top edge, wrap to bottom
                if (y_ind + 1) != NBins:
                    looky = [y_ind-1, y_ind, y_ind+1]
                else:
                    looky = [y_ind-1, y_ind, 0]
                # Reference particle position
                refx = pos[k][0]
                refy = pos[k][1]
                # Loop through surrounding x-index
                for indx in lookx:
                    # Loop through surrounding y-index
                    for indy in looky:
                        # Loop through all particles in that bin
                        for comp in binParts[indx][indy]:
                            # Compute the distance
                            dist = computeDist(refx, refy, pos[comp][0], pos[comp][1])
                            # If potential is on ...
                            if 0.1 < dist <= r_cut:
                                # Compute the x and y components of force
                                fx, fy = computeFLJ(dist, refx, refy, pos[comp][0], pos[comp][1], eps)
                                # This will go into the bulk pressure
                                if clust_size[ids[k]] >= min_size:
                                    bulkSigXX += (fx * (pos[comp][0] - refx))
                                    bulkSigYY += (fy * (pos[comp][1] - refy))
                                # This goes into the gas pressure
                                else:
                                    gasSigXX += (fx * (pos[comp][0] - refx))
                                    gasSigYY += (fy * (pos[comp][1] - refy))

        # Now let's get the area of each phase (by summing bin areas)
        gasBins = 0
        bulkBins = 0
        testIDs = [[0 for b in range(NBins)] for a in range(NBins)]
        for ix in range(0, len(occParts)):
            for iy in range(0, len(occParts)):
                # Is the bin an edge?
                if blurBin[ix][iy] == 1:
                    testIDs[ix][iy] = 0
                    continue
                # Does the bin belong to the dense phase?
                if len(binParts[ix][iy]) != 0:
                    if clust_size[ids[binParts[ix][iy][0]]] > min_size:
                        bulkBins += 1
                        testIDs[ix][iy] = 1
                        continue
                gasBins += 1
                testIDs[ix][iy] = 2

        # The edge length of sufficiently large clusters
        lEdge = Nedges * sizeBin
        # Divide by two because each pair is counted twice
        bulkTrace = (bulkSigXX + bulkSigYY)/2.
        gasTrace = (gasSigXX + gasSigYY)/2.
        # Area of a bin
        binArea = sizeBin * sizeBin
        # Area of each phase
        gasArea = binArea * gasBins
        bulkArea = binArea * bulkBins
        # Pressure of each phase
        gasPress = gasTrace / gasArea
        if bulkArea > 0:
            bulkPress = bulkTrace / bulkArea
            # Surface tension
            surfaceTense = (bulkPress - gasPress) * lEdge
        else:
            bulkPress = 0.
            surfaceTense = 0.

        print("Number of gas bins: ", gasBins)
        print("Gas phase area: ", gasArea)
        print("Number of bulk bins: ", bulkBins)
        print("Bulk phase area: ", bulkArea)
        print("Trace of gas stress tensor: ", gasTrace)
        print("Trace of bulk stress tensor: ", bulkTrace)
        print("Gas phase pressure: ", gasPress)
        print("Bulk phase pressure: ", bulkPress)
        print("Surface tension: ", surfaceTense)

#        # A sanity check on a perfect hcp circle
#        print(Nedges)
#        print(Nedges * sizeBin)
#        x = list(list(zip(*liqPos))[0])
#        y = list(list(zip(*liqPos))[1])
#        diam = max(x) - min(x)
#        circ = diam * np.pi
#        print(circ)
#        print(Nedges * sizeBin / circ)

        # Number of particles in dense phase
        ndense = max(clust_size)

        # Write this to a textfile with the timestep
        g = open(outFile, 'a')
        g.write('{0:.3f}'.format(tst).center(10) + ' ')
        g.write('{0:.3f}'.format(gasArea).center(20) + ' ')
        g.write('{0:.3f}'.format(gasTrace).center(20) + ' ')
        g.write('{0:.3f}'.format(gasPress).center(20) + ' ')
        g.write('{0:.3f}'.format(bulkArea).center(20) + ' ')
        g.write('{0:.3f}'.format(bulkTrace).center(20) + ' ')
        g.write('{0:.3f}'.format(bulkPress).center(20) + ' ')
        g.write('{0:.3f}'.format(surfaceTense).center(20) + ' ')
        g.write('{0:.1f}'.format(lEdge).center(20) + ' ')
        g.write('{0:.0f}'.format(ndense).center(20) + '\n')
        g.close()

# Output an image of the frame we're computing
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow([*zip(*testIDs)], extent=[0, l_box, 0, l_box], origin='lower', aspect='auto')
ax[0].set_aspect('equal')
ax[0].axes.set_xticks([])
ax[0].axes.set_yticks([])
# Plot gas particles
x = list(list(zip(*gasPos))[0])
y = list(list(zip(*gasPos))[1])
ax[1].scatter(x, y, edgecolor='none', s=0.05, c='b')
# Plot liquid particles
if len(liqPos) > 0:
    x = list(list(zip(*liqPos))[0])
    y = list(list(zip(*liqPos))[1])
    ax[1].scatter(x, y, edgecolor='none', s=0.05, c='g')
ax[1].set_xlim(-h_box, h_box)
ax[1].set_ylim(-h_box, h_box)
ax[1].axes.set_xticks([])
ax[1].axes.set_yticks([])
ax[1].set_aspect('equal')
plt.tight_layout(w_pad=0.1)
plt.savefig(imgFile, bbox_inches='tight', pad_inches=0., dpi=250)
plt.close()

max_stars_repo_name: kolbt/whingdingdilly

max_issues_repo_name: kolbt/whingdingdilly

max_forks_repo_name: kolbt/whingdingdilly


text:
'''
This Python file accompanies the manuscript titled
“Assessing the Relative Climate Impact of Carbon dioxide Utilization for Concrete, Chemical and Mineral Production”
authored by Dr. Dwarakanath Ravikumar, Dr. Gregory Keoleian, Dr. Shelie Miller and Dr. Volker Sick.

This file contains code to generate the scatter plot, the distribution of the Climate ROI and the heatmap
for the two pathways – P11: Formic Acid and P2: Concrete – in Figure 2(a), Figure 2(b) and Figure 2(c) of the manuscript, respectively

@author: Dwarakanath Ravikumar
'''
import CO2_ROI_Constants
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from matplotlib import rcParams
import pandas as pd
from matplotlib.ticker import PercentFormatter
from CO2_ROI_Functions import f_get_climate_ROI_for_chemical_list, \
f_get_climate_ROI_for_concrete_list,\
f_get_outranking_calculator_for_proof_of_concept

rcParams['font.family'] = 'Arial'

Climate_ROI_Concrete_Curing_OPC_only_array, percentage_change_comp_strength_Concrete_Curing_OPC_only_array\
                                       =f_get_climate_ROI_for_concrete_list(\
                                       CO2_ROI_Constants.excel_sheet_name_concrete,\
                                       CO2_ROI_Constants.dataset_category_Curing_OPC_only,\
                                       CO2_ROI_Constants.column_name_conventional_first_input_parameter, \
                                       CO2_ROI_Constants.column_name_conventional_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized,\
                                       CO2_ROI_Constants.column_name_CCU_first_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_CCU_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_CCU_concrete_CO2_utilized)

positive_values_Climate_ROI_Concrete_Curing_OPC_only_array = \
len(list(filter(lambda x: x>0,list(Climate_ROI_Concrete_Curing_OPC_only_array.flatten()))))

Climate_ROI_Concrete_Curing_OPC_plus_SCM_array,percentage_change_comp_strength_Concrete_Curing_OPC_plus_SCM_array\
                                       =f_get_climate_ROI_for_concrete_list(\
                                       CO2_ROI_Constants.excel_sheet_name_concrete,\
                                       CO2_ROI_Constants.dataset_category_Curing_OPC_plus_SCM,\
                                       CO2_ROI_Constants.column_name_conventional_first_input_parameter, \
                                       CO2_ROI_Constants.column_name_conventional_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized,\
                                       CO2_ROI_Constants.column_name_CCU_first_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_CCU_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_CCU_concrete_CO2_utilized)

positive_values_Climate_ROI_Concrete_Curing_OPC_plus_SCM_array = \
len(list(filter(lambda x: x>0,list(Climate_ROI_Concrete_Curing_OPC_plus_SCM_array.flatten()))))

Climate_ROI_Concrete_Mixing_OPC_only_array,percentage_change_comp_strength_Concrete_Mixing_OPC_only_array\
                                       =f_get_climate_ROI_for_concrete_list(\
                                       CO2_ROI_Constants.excel_sheet_name_concrete,\
                                       CO2_ROI_Constants.dataset_category_Mixing_OPC_only,\
                                       CO2_ROI_Constants.column_name_conventional_first_input_parameter, \
                                       CO2_ROI_Constants.column_name_conventional_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized,\
                                       CO2_ROI_Constants.column_name_CCU_first_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_CCU_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_CCU_concrete_CO2_utilized)
positive_values_Climate_ROI_Concrete_Mixing_OPC_only_array = \
len(list(filter(lambda x: x>0,list(Climate_ROI_Concrete_Mixing_OPC_only_array.flatten()))))

Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array, percentage_change_comp_strength_Concrete_Mixing_OPC_plus_SCM_array\
                                       =f_get_climate_ROI_for_concrete_list(\
                                       CO2_ROI_Constants.excel_sheet_name_concrete,\
                                       CO2_ROI_Constants.dataset_category_Mixing_OPC_plus_SCM,\
                                       CO2_ROI_Constants.column_name_conventional_first_input_parameter, \
                                       CO2_ROI_Constants.column_name_conventional_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_conventional_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_conventional_CO2_Utilized,\
                                       CO2_ROI_Constants.column_name_CCU_first_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_last_input_parameter,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_LB,\
                                       CO2_ROI_Constants.column_name_CCU_Compressive_strength_range_UB,\
                                       CO2_ROI_Constants.column_name_CCU_CO2_Utilized_Duplicate,\
                                       CO2_ROI_Constants.column_name_CCU_concrete_CO2_utilized)

positive_values_Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array = \
len(list(filter(lambda x: x>0,list(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten()))))

climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITHOUT_EOC,\
climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC\
                                        =f_get_climate_ROI_for_chemical_list(\
                                        CO2_ROI_Constants.formic_acid,\
                                        CO2_ROI_Constants.Catalytic_reduction_of_CO2_with_H2,\
                                        CO2_ROI_Constants.CO2_emissions_formic_acid_conventional_mean, \
                                        CO2_ROI_Constants.CO2_emissions_formic_acid_conventional_sd,\
                                        CO2_ROI_Constants.lognormal_distribution)

positive_values_climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array = \
len(list(filter(lambda x: x>0,list(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITHOUT_EOC.flatten()))))

positive_values=[]

#positive_values.append(positive_values_Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array)
positive_values.append(positive_values_Climate_ROI_Concrete_Mixing_OPC_only_array)
#positive_values.append(positive_values_Climate_ROI_Concrete_Curing_OPC_only_array)
#positive_values.append(positive_values_Climate_ROI_Concrete_Curing_OPC_plus_SCM_array)
positive_values.append(positive_values_climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array)

print(np.round(((np.array(positive_values)/CO2_ROI_Constants.number_of_samples)*100)))

#*********** SCATTER PLOT START ********************************

scatter_plot_concept_fig=plt.figure(figsize=(30,40))
x_axis_array=[]
#x_axis_array=list(percentage_change_comp_strength_Concrete_Curing_OPC_only_array.flatten())
#(x_axis_array).extend(list(percentage_change_comp_strength_Concrete_Curing_OPC_plus_SCM_array.flatten()))
#(x_axis_array).extend(list(percentage_change_comp_strength_Concrete_Mixing_OPC_only_array.flatten()))
(x_axis_array).extend(list(percentage_change_comp_strength_Concrete_Mixing_OPC_plus_SCM_array.flatten()))
(x_axis_array).extend(len(list(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten()))*[0])

y_axis_array=[]
#y_axis_array=list(Climate_ROI_Concrete_Curing_OPC_only_array.flatten())
#(y_axis_array).extend(list(Climate_ROI_Concrete_Curing_OPC_plus_SCM_array.flatten()))
#(y_axis_array).extend(list(Climate_ROI_Concrete_Mixing_OPC_only_array.flatten()))
(y_axis_array).extend(list(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten()))
(y_axis_array).extend(list(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten()))

color_array=[]
#color_array=len(Climate_ROI_Concrete_Curing_OPC_only_array.flatten())*[CO2_ROI_Constants.concrete_Curing_OPC_only_color]
#color_array.extend(len(Climate_ROI_Concrete_Curing_OPC_plus_SCM_array.flatten())*[CO2_ROI_Constants.concrete_Curing_OPC_plus_SCM_color])
#color_array.extend(len(Climate_ROI_Concrete_Mixing_OPC_only_array.flatten())*[CO2_ROI_Constants.concrete_Mixing_OPC_only_color])
color_array.extend(len(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten())*[CO2_ROI_Constants.concrete_Mixing_OPC_plus_SCM_color])
color_array.extend(len(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())*[CO2_ROI_Constants.formic_acid_Catalytic_reduction_of_CO2_with_H2_WITH_EOC_color])

size_array=[]
(size_array).extend(CO2_ROI_Constants.scatter_concrete_mixing_dot_size*np.ones(len(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten())))
(size_array).extend(CO2_ROI_Constants.scatter_formic_acid_dot_size*np.ones(len(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())))

fc_array=[]
fc_array.extend(len(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten())*[CO2_ROI_Constants.concrete_Mixing_OPC_plus_SCM_color])
fc_array.extend(len(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())*['none'])

ec_array=[]
ec_array.extend(len(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten())*[CO2_ROI_Constants.concrete_Mixing_OPC_plus_SCM_color])
ec_array.extend(len(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())*[CO2_ROI_Constants.formic_acid_Catalytic_reduction_of_CO2_with_H2_WITH_EOC_color])

Concrete_Mixing_OPC_only_array_scatter_plot = {\
         'x_axis_array': x_axis_array,\
        'y_axis_array': y_axis_array,\
        'color_array':color_array,\
        'size_array':size_array,\
        'fc_array':fc_array,\
        'ec_array':ec_array\
        }

df_scatter_plot_data = pd.DataFrame(Concrete_Mixing_OPC_only_array_scatter_plot, columns = ['x_axis_array', 'y_axis_array','color_array','size_array','fc_array','ec_array'])

plt.scatter(x='x_axis_array',\
               y='y_axis_array',\
               data=df_scatter_plot_data,\
               s='size_array',\
               facecolors='fc_array',\
               c='color_array',\
               marker="s")

ax = plt.axes()
ax.set_ylim(-100,100)
ax.set_xlim(-100,100)
ax.set_ylabel('')
ax.set_xlabel('')
plt.yticks(fontsize=60)
plt.xticks(fontsize=60)
ax.axhline(y=0, color='black',linewidth=0.5)
ax.axvline(x=0, color='black',linewidth=0.5)
ax.tick_params(pad=20)
scatter_plot_concept_fig.savefig('Proof_of_Concept_Scatter_Plot.png')

#*********** SCATTER PLOT END ********************************

#*************** KDE PLOT START ****************************

kde_plot_fig=plt.figure(figsize=(30,40))

#sns.kdeplot(Climate_ROI_Concrete_Curing_OPC_only_array.flatten(), vertical=True,linewidth=CO2_ROI_Constants.kde_curve_width, shade=True,color=CO2_ROI_Constants.concrete_Curing_OPC_only_color)
#sns.kdeplot(Climate_ROI_Concrete_Curing_OPC_plus_SCM_array.flatten(), vertical=True,linewidth=CO2_ROI_Constants.kde_curve_width, shade=True,color=CO2_ROI_Constants.concrete_Curing_OPC_plus_SCM_color)
#sns.kdeplot(Climate_ROI_Concrete_Mixing_OPC_only_array.flatten(), vertical=True,linewidth=CO2_ROI_Constants.kde_curve_width, shade=True,color=CO2_ROI_Constants.concrete_Mixing_OPC_only_color)
sns.kdeplot(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten(), vertical=True,linewidth=CO2_ROI_Constants.kde_curve_width, shade=True,color=CO2_ROI_Constants.concrete_Mixing_OPC_plus_SCM_color)
sns.kdeplot(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten(), vertical=True,linewidth=CO2_ROI_Constants.kde_curve_width, shade=True,color=CO2_ROI_Constants.formic_acid_Catalytic_reduction_of_CO2_with_H2_WITH_EOC_color)

ax = plt.axes()
ax.set_xticks([])
ax.set_ylim(-100,100)
plt.yticks(fontsize=60)
ax.axhline(y=0, color='black',linewidth=0.5)
#ax.axvline(y=0, color='black',linewidth=0.5)

kde_plot_fig.savefig('Proof_of_Concept_KDE_Plot.png')

#*************** KDE PLOT END ****************************

super_array=[]

#It is important to note that the ROI arrays below are appended in the descending order based on
#the climate ROI
#The  arrays of the products at the top (appended first) have the greatest ROI.
#Do NOT change this order as this order helps achieve the 'green' rows on top versus 'red'
#rows at the bottom for the heatmap

super_array.append(Climate_ROI_Concrete_Mixing_OPC_plus_SCM_array.flatten())
super_array.append(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())

#super_array.append(climate_ROI_formic_acid_Catalytic_reduction_of_CO2_with_H2_array_WITH_EOC.flatten())
#super_array.append(Climate_ROI_Concrete_Mixing_OPC_only_array.flatten())
#super_array.append(Climate_ROI_Concrete_Curing_OPC_only_array.flatten())
#super_array.append(Climate_ROI_Concrete_Curing_OPC_plus_SCM_array.flatten())

#get the outranking % scores which are included in the heatmap
outranking_array=f_get_outranking_calculator_for_proof_of_concept(super_array)

matrix = np.zeros(outranking_array.shape, dtype=bool)
np.fill_diagonal(matrix, True)

#Uncomment the following two lines for triangular heatmap
#matrix = np.zeros_like(outranking_array, dtype=np.bool)
#matrix[np.triu_indices_from(matrix)] = True

heatmap_fig=plt.figure(figsize=(30,10))

#Properties for the colorbar. Make the colorbar horizontal (orientation), make it the same
#dimensions as the heatmap (shrink), decrease space between x-axis of heatmap and colorbar (pad)
#and decrease height when compared to length of colorbar (aspect)
cbar_kws1={ 'orientation': 'horizontal', 'shrink': 1, 'pad':0.03, 'aspect':40}

ax3 = sns.heatmap(outranking_array, linewidth=2,linecolor='black', annot = True, \
                  cmap= 'RdYlGn',square=False, annot_kws={'size':52, 'weight': 'bold'},fmt=".0%",\
                  mask=matrix,cbar_kws=cbar_kws1)

#Set the increments for the colorbar in increments of 20% from 0 to 100%
cbar = ax3.collections[0].colorbar
cbar.set_ticks([0,0.1,0.2,0.3,0.4,0.5,0.6,0.71])
cbar.ax.xaxis.set_major_formatter(PercentFormatter(1, 0))
cbar.ax.tick_params(labelsize=52)
#cbar.ax.tick_params(weight='bold')

ax = plt.axes()
#ax3.xaxis.set_ticks_position('top')
ax.set_xticks([])
ax.set_yticks([])

heatmap_fig.savefig('Proof_of_Concept_heatmap.png')
#***************HEATMAP OUTRANKING  PLOT END ****************************