#!/usr/bin/python

"""Module Description

  This module is used to measure various patterns of the dynamic nucleosome regions, including:
    1. occupancy change: the density changing of this region
    2. reposition variation: the maxium position changing nucleosome of this region ( use postioning degree )
    3. positioning degree change: use sliding window to calculate region's positioning degree and the nucleosome which has
    the largest positioning degree

Copyright (c) 2011 Kai Fu  < kelvinfu.tju@gmail.com >

This code is free software; you can redistribute it and/or modify it
under the terms of the Artistic License (see the file COPYING included
with the distribution).

@status:  alpha
@version: $Released$
@author:  Kai Fu
@contact: kelvinfu.tju@gmail.com
"""

# -----------------------------------
# own python modules
# -----------------------------------

from DINUP.Peak_Detect import *
from DINUP.ks_test import ks_2samp


class Peak_Pattern:
    
    """ Class to assign different physical properties to a dynamic nucleosome region.
    Includes: occupancy change, repostioned variation, positioning degree change (fuzziness to phasing)
    """
    
    def __init__(self, region_start = None, region_end = None, treat = None, control = None, coverage_treat = None, coverage_control = None, pvalue_location = None):
	
	""" Initilize the Peak_Pattern object
	"""
	
	self.region_start = region_start
	self.region_end = region_end
	self.treat = treat
	self.control = control
	self.coverage_treat = coverage_treat
	self.coverage_control = coverage_control
	self.pvalue_location = pvalue_location
    
    def get_region_reads(self):
	
	""" get reads of the identified differential nucleosome positioning regions.
	"""
	
	treat_region_reads = {}
	control_region_reads = {}
	
	for key in self.region_start.keys():
	    treat_region_reads[key] = [[] for k in range(len(self.region_start[key]))]
	    control_region_reads[key] = [[] for k in range(len(self.region_start[key]))]
	
	for key in self.region_start.keys():
	    for i in range(len(self.region_start[key])):
		for j in range(len(self.treat[key])):
		    if self.treat[key][j] >= self.region_start[key][i] and self.treat[key][j] <= self.region_end[key][i]:
			treat_region_reads[key][i].append(self.treat[key][j])
		
		for k in range(len(self.control[key])):
		    if self.control[key][k] >= self.region_start[key][i] and self.control[key][k] <= self.region_end[key][i]:
			control_region_reads[key][i].append(self.control[key][k])
	
	return treat_region_reads, control_region_reads

    
    def occupancy_change(self, treat_region_reads, control_region_reads):
	
	""" fold change of reads number between treatment and control
	"""
	
	fold_change = {}
	
	for key in treat_region_reads.keys():
	    for i in xrange(len(treat_region_reads[key])):
		### the reads within a region should be max(reads by change, reads occur in this region)  ###
		fold = ((float(max(len(treat_region_reads[key][i]), self.coverage_treat, 1)))/max(len(control_region_reads[key][i]),self.coverage_control,1))/(float(max(self.coverage_treat,1))/max(self.coverage_control,1))
		
		if fold_change.has_key(key):
		    fold_change[key].append(fold)
		else:
		    fold_change[key] = [fold]
		    
	return fold_change

    
    def positioning_degree(self, treat_region_reads, control_region_reads):
	
	""" measures nucleosome phasing to fuzziness or fuzziness to phasing
	"""
	
	window_160bp_treat_reads = []
	window_160bp_control_reads = []
	window_20bp_treat_reads = []
	window_20bp_control_reads = []
	degree_20bp_window_treat = []
	degree_20bp_window_control = []
	degree_160bp_window_treat = []
	degree_160bp_window_control = []

	degree_region = {}
	
	for key in self.region_start.keys():
	    for i in xrange(len(self.region_start[key])):
		for j in xrange(self.region_start[key][i], self.region_end[key][i], 10):
		    window_160bp_start = j
		    window_160bp_end = j + 160
		    for l in xrange(len(treat_region_reads[key][i])):
			if treat_region_reads[key][i][l] >= j and treat_region_reads[key][i][l] <= j+160:
			    window_160bp_treat_reads.append(treat_region_reads[key][i][l])
		    for l in xrange(len(control_region_reads[key][i])):
			if control_region_reads[key][i][l] >= j and control_region_reads[key][i][l] <= j+160:
			    window_160bp_control_reads.append(control_region_reads[key][i][l])

		    for k in xrange(window_160bp_start, window_160bp_end, 20):
			window_20bp_start = k
			window_20bp_end = k + 20
			for m in xrange(len(window_160bp_treat_reads)):
			    if window_160bp_treat_reads[m] >= k and window_160bp_treat_reads[m] <= k+20:
				window_20bp_treat_reads.append(window_160bp_treat_reads[m])
			
			degree_20bp_window_treat.append(float(len(window_20bp_treat_reads))/max(len(window_160bp_treat_reads), 2*self.coverage_treat,1))
			for m in xrange(len(window_160bp_control_reads)):
			    if window_160bp_control_reads[m] >= k and window_160bp_control_reads[m] <= k+20:
				window_20bp_control_reads.append(window_160bp_control_reads[m])
			degree_20bp_window_control.append(float(len(window_20bp_control_reads))/max(len(window_160bp_control_reads), 2*self.coverage_control,1))
			window_20bp_treat_reads = []
			window_20bp_control_reads = []
		    
		    degree_20bp_window_treat.sort()
		    degree_20bp_window_control.sort()
		    degree_160bp_window_treat.append(degree_20bp_window_treat[-1])
		    degree_160bp_window_control.append(degree_20bp_window_control[-1])
		    
		    window_160bp_treat_reads = []
		    window_160bp_control_reads = []
		    degree_20bp_window_treat = []
		    degree_20bp_window_control = []
		
		if degree_region.has_key(key):
		    degree_region[key].append(sum(degree_160bp_window_treat)/len(degree_160bp_window_treat)-sum(degree_160bp_window_control)/len(degree_160bp_window_control))
		    #degree_region_control[key].append(sum(degree_160bp_window_control)/len(degree_160bp_window_control))
		else:
		    degree_region[key] = [sum(degree_160bp_window_treat)/len(degree_160bp_window_treat)-sum(degree_160bp_window_control)/len(degree_160bp_window_control)]
		    #degree_region_control[key] = [sum(degree_160bp_window_control)/len(degree_160bp_window_control)]
		
		degree_160bp_window_treat = []
		degree_160bp_window_control = []
		
	return degree_region

    
    def reposition_variation(self, treat_region_reads, control_region_reads):
	
	""" measure the maxium location change of the nucleosomes in the indentified regions
	"""
	
	window_160bp_treat_reads = []
	window_160bp_control_reads = []
	
	window_20bp_treat_reads = []
	window_20bp_control_reads = []
	
	position_degree_treat = []
	position_degree_control = []
	position_degree_treat_local = []
	position_degree_control_local = []
	
	max_degree_treat = []
	max_degree_control = []
	variation = []
	
	position_variation = {}
	
	for key in self.region_start.keys():
	    for i in xrange(len(self.region_start[key])):
		
		### First choose the most well positioned nucleosome (the 160bp window has most reads) in treat and control region ###
		for j in xrange(self.region_start[key][i], self.region_end[key][i]-160, 10):
		    window_20bp_start = j + 70
		    window_20bp_end = j + 90          # this is 40bp small sliding window
		    window_160bp_start = j
		    window_160bp_end = j + 160
		    
		    for l in xrange(len(treat_region_reads[key][i])):
			if treat_region_reads[key][i][l] >= window_20bp_start and treat_region_reads[key][i][l] <= window_20bp_end:
			    window_20bp_treat_reads.append(treat_region_reads[key][i][l])
		    for l in xrange(len(control_region_reads[key][i])):
			if control_region_reads[key][i][l] >= window_20bp_start and control_region_reads[key][i][l] <= window_20bp_end:
			    window_20bp_control_reads.append(control_region_reads[key][i][l])
		    for l in xrange(len(treat_region_reads[key][i])):
			if treat_region_reads[key][i][l] >= window_160bp_start and treat_region_reads[key][i][l] <= window_160bp_end:
			    window_160bp_treat_reads.append(treat_region_reads[key][i][l])
		    for l in xrange(len(control_region_reads[key][i])):
			if control_region_reads[key][i][l] >= window_160bp_start and control_region_reads[key][i][l] <= window_160bp_end:
			    window_160bp_control_reads.append(control_region_reads[key][i][l])
		    
		    position_degree_treat.append(1.0*len(window_20bp_treat_reads)/max(len(window_160bp_treat_reads),2*self.coverage_treat))
		    position_degree_control.append(1.0*len(window_20bp_control_reads)/max(len(window_160bp_control_reads), 2*self.coverage_control))
		    
		    window_20bp_treat_reads = []
		    window_20bp_control_reads = []
		    window_160bp_treat_reads = []
		    window_160bp_control_reads = []
		
		max_number_treat = 0
		max_number_control = 0
		
		for k in xrange(len(position_degree_treat)):
		    
		    if position_degree_treat[k] >= 0.2: 
			position_degree_treat_local.append(position_degree_treat[k])
		    elif position_degree_treat[k] < 0.2 and len(position_degree_treat_local) >= 3:
			for n in xrange(len(position_degree_treat_local)):
			    if position_degree_treat_local[n] > max_number_treat:
				max_degree_treat_local = n
				max_number_treat = position_degree_treat_local[n]
			    else:
				pass
			max_degree_treat.append(k + n - len(position_degree_treat_local))  # the location of the local most positioned nucleosome
			position_degree_treat_local = []
		    else:
			position_degree_treat_local = []
		
		
		for m in xrange(len(position_degree_control)):
		    if position_degree_control[m] >= 0.2:
			position_degree_control_local.append(position_degree_control[m])
		    elif position_degree_control[m] < 0.2 and len(position_degree_control_local) >= 3:
			for n in xrange(len(position_degree_control_local)):
			    if position_degree_control_local[n] > max_number_control:
				max_degree_control_local = n
				max_number_control = position_degree_control_local[n]
			    else:
				pass
			max_degree_control.append(m + n - len(position_degree_control_local))  # the location of the local most positioned nucleosome
			position_degree_control_local = []
		    else:
			position_degree_control_local = []
		
		#print max_degree_treat
		#print max_degree_control
		
		for p in xrange(len(max_degree_treat)):
		    for q in xrange(len(max_degree_control)):
			if abs(max_degree_treat[p] - max_degree_control[q]) <= 10:
			    variation.append(abs(max_degree_treat[p] - max_degree_control[q]))
			else:
			    pass
		
		variation.sort()
		if position_variation.has_key(key):
		    if len(variation) == 0:
			position_variation[key].append(0)
		    else:
			position_variation[key].append(variation[-1]*10)
		else:
		    if len(variation) == 0:
			position_variation[key] = [0]
		    else:
			position_variation[key] = [variation[-1]*10]
		
		position_degree_treat = []
		position_degree_control = []
		
		max_degree_treat = []
		max_degree_control = []
		
		variation = []
	
	return position_variation