##########################################################################
# DataImportUtils.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.     #
##########################################################################

"""Provides functions for importing data for use in clustering.

Major functions (populate_raw_replicants, populate_normal_replicants,
and populate_genes) are used to import data from the text files in the 
TextFiles directory and transform the data into a form usable by the
clustering algorithm.  A minor function, load_from_file, is a specialized
function called from populate_raw_replicants and should not be called from
anywhere else."""

# MASTER IMPORT STATEMENT

from Simcal1 import *
	
# MAJOR FUNCTIONS

def populate_raw_replicants(pathname):
	"""Get data from text files in the directory specified by pathname.
	
	The data from the files will be run through the RawDataIn table (note
	that this table is a utility table only, and will not permanently hold
	any useful data) and inserted in the RawReplicants table."""
	
	### function variable definitions
	
	replicants = ("1", "2", "3")
	types = ("control", "iso", "psr")
	times = ("05", "20", "85")
						   
	### populate table
	
	# user notification
	print "\nBeginning RawReplicants population ... " + asctime()
	
	# database pre-cleanup
	setq("delete from RawReplicants");
	
	# main loop -- form names like "rep2_iso_05.txt"
	for replicant in replicants:
		basefilename = "rep" + replicant + "_"
		for type in types:
			if type == "control":
				filename = basefilename + type + ".txt"
				load_from_file(pathname, filename, replicant, type)
			else:
				for time in times:
					filename = basefilename + type + "_" + time + ".txt"
					load_from_file \
						(pathname, filename, replicant, type, time)
					
	### notify user and get us out of here
	
	print "Raw replicants populated: " + asctime() + "\n"
	return
	
	
def populate_normal_replicants():
	"""Normalize data in RawReplicants and put it in NormalReplicants."""

	# table cleanup
	setq("delete from NormalReplicants")
						
	# get rid of unreliable data
	setq("delete from RawReplicants where Pvalue > 0.05")
	
	# user notification
	print "\nBeginning NormalReplicants population ... " + asctime()
	
	# normalization and insertion
	# "1-Pvalue" is raised to the power log_0.95(0.05) because we want
	# genes with p-values of 0.05 to have reliabilities of 0.05 (5%).
	setq("""insert into NormalReplicants
			select Replicant, 
				   Id, 
				   AffyName,
				   pow((1-Pvalue), log(0.95, 0.05)) as Reliability,
				   log2(pow(Psr1,2)/(Control*Iso1)) as Exp1,
				   log2(pow(Psr2,2)/(Control*Iso2)) as Exp2,
				   log2(pow(Psr3,2)/(Control*Iso3)) as Exp3,
				   AffyAnnotation 
		  	from   RawReplicants
		  	order  by Replicant, Id""")
					
	# notify user and get us out of here
	print "Normalized replicants populated: " + asctime() + "\n"
	return
	
	
def populate_genes():
	"""Average data from RawReplicants and insert results in Genes."""
	
	### user notification
	
	print "\nBeginning Genes population ... " + asctime()
	
	### table cleanup
	
	setq("delete from Genes")
	
	### averaging and insertion
		
	# modify table to do auto-incrementing of ID's
	setq("alter table Genes modify Id int auto_increment not null");
	setq("alter table Genes auto_increment = 1");
	
	# table population, ordered by variability
	setq("""insert into Genes
		  	select NULL,
				   sum(Reliability*Exp1)/sum(Reliability) as Exp1,
				   sum(Reliability*Exp2)/sum(Reliability) as Exp2,
				   sum(Reliability*Exp3)/sum(Reliability) as Exp3,
				   (abs(sum(Reliability*Exp1)) +
				    abs(sum(Reliability*Exp2)) +
				    abs(sum(Reliability*Exp3))) /
				   (3.0 * abs(sum(Reliability))) as Variability,
				   "" as UserAnnotation,
				   avg(Reliability) as Reliability,
				   AffyName,
				   AffyAnnotation
		    from   NormalReplicants
		    group  by Id
		    order  by Variability""")
		    
	# modify table back to regular int ID's, start ID's at 0	    
	setq("alter table Genes modify Id int not null")
	setq("update Genes set Id = Id - 1");
					
	### notify user and get us out of here
	
	print ("Genes populated -- ready for clustering: " + 
		   asctime() + "\n")
	return


# MINOR FUNCTION

def load_from_file(pathname, filename, replicant, type, time = None):

	# function variable definitions
	id = 0 # record id baseline (incremented before each use)
	nameDict = {"iso":"Iso", "psr":"Psr",
				"05":"1", "20":"2", "85":"3"}
		# what nameDict does is establish mappings from the names used
		# for the text files to those used for fields in RawReplicants

	# do pre-cleanup and put data in
	fullname = pathname + filename
	setq("delete from RawDataIn");
	setq("load data infile '" + fullname + "' into table RawDataIn");
	
	# get records
	recs = mobq("""select *
				   from	  RawDataIn
				   order  by AffyName""") 
					  
	# FUDGE_FACTOR is used here because the normalization algorithm, 
	# used in populate_normal_replicants, can't handle expression values
	# of 0 (of which there are a few).
	
	if type == "control": # start new row for each record
				 	 
		# insert records
		for rec in recs:
			id += 1 # start with Id 1, not 0
			setq("""insert into RawReplicants
					set Replicant = %s,
						Id = %s,
						AffyName = %s,
						Pvalue = %s,
						Control = %s,
						Iso1 = 0,
						Iso2 = 0,
						Iso3 = 0,
						Psr1 = 0,
						Psr2 = 0,
						Psr3 = 0,
						AffyAnnotation = %s""",
				 (replicant, id, rec.AffyName, rec.Pvalue,
				  rec.Expression + FUDGE_FACTOR, rec.AffyAnnotation))
				
	else: # update records with appropriate values
	
		# set up utility variable
		typeTime = nameDict[type] + nameDict[time]
		
		# update records -- note that for any replicant, its Pvalue is
		# defined as the greatest Pvalue of any of its measurements
		for rec in recs:
			id += 1 # start with Id 1, not 0
			setq("""update RawReplicants
					set    """ + typeTime + """ = %s,
						   Pvalue = greatest(Pvalue, %s)
					where  Replicant = %s and 
						   Id = %s""", 
			   	 (rec.Expression + FUDGE_FACTOR, 
			   	  rec.Pvalue, replicant, id))
						    
	# notify user and get us out of here
	print "\t" + filename + " loaded: " + asctime()
	return