##########################################################################
# ClusteringControl.py:  Python source file for the Simcal1 module.      #
# Copyright (c) 2004 Daniel Dvorkin.                                     #
# Distributed under the Artistic License, WITHOUT "optional clause 8"    #
# (see http://www.opensource.org/licenses/artistic-license.php).         #
# All rights not granted in this license are reserved by the author.     #
##########################################################################

"""Defines functions do_find_gene_relationships, do_clustering, and
do_find_cluster_successions for control of clustering algorithm."""

# MASTER IMPORT STATEMENT

from Simcal1 import *
	
# CLUSTERING ALGORITHM FUNCTIONS

def do_find_gene_relationships():
	"""Find gene leaders and followers, determine successions, 
	and find neighboring genes up to MAX_RANK."""
	
	# prep work
	print "Begin gene relationship finding operations %s" % asctime()
	Gene.do_cleanup()
	Genes = Gene.get_group()
	
	# extract location information
	print "\tExtracting location information ... %s" % asctime()
	Locations = array([gene.location() for gene in Genes])
	genes = Locations.shape[0] # number of genes
	exps = Locations.shape[1] # time points per gene
	LateLocations = Locations[:, 1:]
	EarlyLocations = Locations[:, 0:exps-1]
	
	# find correlations, adjust for similarities (close genes should not 
	# lead or follow each other) and make sure no gene is marked as its
	# own neighbor, leader, or follower by zeroing diagonals
	print "\tFinding correlations ... %s" % asctime()
	print "\t\tFinding neighbor correlations ... %s" % asctime()
	Neighbors = max_zero(corr_coef(Locations, Locations))
	print "\t\tFinding raw succession correlations ... %s" % asctime()
	Leaders = corr_coef(LateLocations, EarlyLocations)
	print "\t\tAdjusting succession correlations ... %s" % asctime()
	Followers = transpose(Leaders)
	Leaders -= Neighbors
	Followers -= Neighbors
	print "\t\tRemoving self-association ... %s" % asctime()
	Mask = eye(genes) * 2
	Neighbors -= Mask
	Leaders -= Mask
	Followers -= Mask
	
	# find non-parametric relationships (we're flipping left-
	# to-right because we want to get highest arguments first)
	print "\tFinding relationships ... %s" % asctime()
	print "\t\tFinding neighbor relationships ... %s" % asctime()
	Neighbors = fliplr(argsort(Neighbors))
	print "\t\tFinding leader relationships ... %s" % asctime()
	Leaders = fliplr(argsort(Leaders))
	print "\t\tFinding follower relationships ... %s" % asctime()
	Followers = fliplr(argsort(Followers))
	
	# break out into lists of tuples for insertion
	print "\tExtracting relationships ... %s" % asctime()
	NeighborTuples = [(id, Neighbors[id, rank-1], rank) 
		for id in xrange(genes) for rank in prange(MAX_RANK)]
	SuccessionTuples = [(Leaders[id, rank-1], id, rank) 
		for id in xrange(genes) for rank in prange(MAX_RANK)]
	SuccessionTuples += [(id, Followers[id, rank-1], rank) 
		for id in xrange(genes) for rank in prange(MAX_RANK)]
		
	# save data, user notification
	print "\tSaving relationships ... %s" % asctime()
	Gene.set_neighbors(NeighborTuples)
	Gene.set_individual_successions(SuccessionTuples)
	
	# user notification
	print "\tCalculating final relationships ... %s" % asctime()
	Gene.set_successions()
		 
def do_clustering():
	
	# prep work
	print "Begin clustering operations %s" % asctime()
	Cluster.do_cleanup()
	GeneNeighbors = Gene.get_neighbors()
	NumOfGenes = GeneNeighbors.shape[1]
	GeneClusterIndexes = range(NumOfGenes)
	ClusterGenes = [[i] for i in GeneClusterIndexes]
	
	# cluster assignment loop
	for rank in prange(MAX_RANK):
		
		print "\tAssigning clusters, rank %s ... %s" % (rank, asctime())
	
		# cyle through genes
		for GeneId in xrange(NumOfGenes):
		
			# get ID's into easy forms
			i = GeneId
			j = GeneNeighbors[rank-1, i]
			
			# cluster assignment
			if GeneClusterIndexes[j] != GeneClusterIndexes[i]:
				# merge higher-ID cluster with lower-ID cluster
				localGCI = array([GeneClusterIndexes[i], 
								  GeneClusterIndexes[j]])
				live = min(localGCI)
				dead = max(localGCI)
			 	# update GeneClusterIndexes to reflect merge
				for id in ClusterGenes[dead]:
					GeneClusterIndexes[id] = GeneClusterIndexes[live]
				# merge old cluster into new
				ClusterGenes[live] += ClusterGenes[dead]
				# empty out old cluster
				ClusterGenes[dead] = []
			# else:  already in the same cluster
			
		# create clusters, terminate if there's no reason to go on
		if Cluster.do_create(ClusterGenes, rank) == 1:
			break
		
def do_find_succession(): 
	"Find cluster successions."""
	
	# prep work
	print "Begin cluster succession operations %s" % asctime()
	Cluster.do_succession_cleanup()
	ClusterGenes = Cluster.get_genes()
	NumOfClusterGenes = len(ClusterGenes)
	
	# cluster succession loop
	print "\tBegin succession finding loop %s" % asctime()
	pctbreaks = ceil(NumOfClusterGenes / 10.0)
	for i in xrange(NumOfClusterGenes):
	
		Cluster.set_individual_successions(ClusterGenes[i])
		
		# user notification of progress, as this can take a while
		if i > 0 and i % pctbreaks == 0:
			pct = 10 * i / pctbreaks
			print "\t\tSuccession finding loop %i%% done, %s" \
				  % (pct, asctime())  
				  
	print "\t\tSuccession finding loop 100%% done, %s" % asctime()
	
	# save succession information
	print "\tSaving cluster successions ... %s" % asctime()
	Cluster.set_successions()