mi@0: #!/usr/bin/env python
mi@0: # encoding: utf-8
mi@0: """
mi@0: novelty.py
mi@0: 
mi@0: Created by mi tian on 2015-04-02.
mi@0: Copyright (c) 2015 __MyCompanyName__. All rights reserved.
mi@0: """
mi@0: 
mi@0: import sys, os
mi@0: import numpy as np
mi@0: from scipy.signal import correlate2d, convolve2d
mitian@1: import matplotlib.pyplot as plt
mitian@7: from utils.PeakPickerUtil import PeakPicker
mi@0: 
mitian@7: 	
mitian@7: peak_picker = PeakPicker()
mitian@7: peak_picker.params.alpha = 9.0 # Alpha norm
mitian@7: peak_picker.params.delta = 0.0 # Adaptive thresholding delta
mitian@7: peak_picker.params.QuadThresh_a = (100 - 20.0) / 1000.0
mitian@7: peak_picker.params.QuadThresh_b = 0.0
mitian@7: peak_picker.params.QuadThresh_c = (100 - 20.0) / 1500.0
mitian@7: peak_picker.params.rawSensitivity = 20
mitian@7: peak_picker.params.aCoeffs = [1.0000, -0.5949, 0.2348]
mitian@7: peak_picker.params.bCoeffs = [0.1600,	 0.3200, 0.1600]
mitian@7: peak_picker.params.preWin = 5
mitian@7: peak_picker.params.postWin = 5 + 1
mitian@7: peak_picker.params.LP_on = True
mitian@7: peak_picker.params.Medfilt_on = True
mitian@7: peak_picker.params.Polyfit_on = True
mitian@7: peak_picker.params.isMedianPositive = False
mitian@18: peak_picker.params.gt = None
mitian@18: 
mi@0: def getNoveltyCurve(ssm, kernel_size, normalise=False):
mi@0: 	'''Return novelty score from ssm.'''
mi@0: 
mi@0: 	kernel_size = int(np.floor(kernel_size/2.0) + 1)
mi@0: 	stripe = getDiagonalSlice(ssm, kernel_size)
mi@0: 	kernel = gaussian_kernel(kernel_size)
mi@0: 	xc = convolve2d(stripe,kernel,mode='same')
mi@0: 	xc[abs(xc)>1e+10]=0.00001			
mi@0: 	
mi@0: 	novelty = xc[int(np.floor(xc.shape[0]/2.0)),:]	
mi@0: 	novelty = [0.0 if (np.isnan(x) or np.isinf(x) or x > 1e+100) else x for x in novelty]
mi@0: 	
mi@0: 	if normalise:
mi@0: 		novelty = (novelty - np.min(novelty)) / (np.max(novelty) - np.min(novelty))
mi@0: 	return novelty	
mi@0: 
mitian@8: def getNoveltyDiff(novelty, N=1, relative=False):
mitian@8: 	'''Return the second order differece in the novelty curve.'''
mitian@8: 	
mitian@8: 	diff = np.zeros_like(novelty)
mitian@8: 	diff[:-N] = np.diff(novelty, n=N)
mitian@8: 	
mitian@8: 	if relative:
mitian@8: 		diff /= novelty 
mitian@8: 		
mitian@8: 	return diff
mitian@8: 	
mi@0: def getDiagonalSlice(ssm, width):
mi@0: 	''' Return a diagonal stripe of the ssm given its width, with 45 degrees rotation. 
mi@0: 	Note: requres 45 degrees rotated kernel also.'''
mi@0: 	w = int(np.floor(width/2.0))
mi@0: 	length = len(np.diagonal(ssm))
mi@0: 	stripe = np.zeros((2*w+1,length))
mi@0: 	# print 'diagonal', length, w, stripe.shape
mi@0: 	for i in xrange(-w, w+1) :		
mi@0: 		stripe[w+i,:] = np.hstack(( np.zeros(int(np.floor(abs(i)/2.0))), np.diagonal(ssm,i), np.zeros(int(np.ceil(abs(i)/2.0))) ))
mi@0: 	return stripe
mi@0: 
mi@0: def gaussian_kernel(size):
mi@0: 	'''Create a gaussian tapered 45 degrees rotated checkerboard kernel. 
mi@0: 	TODO: Unit testing: Should produce this with kernel size 3:
mi@0: 	0.1353	 -0.3679	0.1353
mi@0: 	0.3679	  1.0000	0.3679
mi@0: 	0.1353	 -0.3679	0.1353
mi@0: 	'''
mi@0: 	n = float(np.ceil(size / 2.0))
mi@0: 	kernel = np.zeros((size,size))
mi@0: 	for i in xrange(1,size+1) :
mi@0: 		for j in xrange(1,size+1) :
mi@0: 			gauss = np.exp( -4.0 * (np.square( (i-n)/n ) + np.square( (j-n)/n )) )
mi@0: 			# gauss = 1			
mi@0: 			if np.logical_xor( j - n > np.floor((i-n) / 2.0), j - n > np.floor((n-i) / 2.0) ) :
mi@0: 				kernel[i-1,j-1] = -gauss
mi@0: 			else:
mi@0: 				kernel[i-1,j-1] = gauss
mi@0: 				
mi@0: 	return kernel	
mi@0: 	
mitian@17: def process(ssm, peak_picker=peak_picker, kernel_size=48, normalise=False, plot=False):
mi@0: 	'''Detect segment boundaries in the ssm.'''
mitian@7: 	# peak_picker for the 1st round boudary detection
mitian@7: 
mitian@7: 	
mi@0: 	novelty = getNoveltyCurve(ssm, kernel_size, normalise=False)
mi@0: 	smoothed_novelty, novelty_peaks = peak_picker.process(novelty)
mi@0: 	
mitian@1: 	if plot:
mitian@1: 		plot_detection(smoothed_novelty, novelty_peaks)
mitian@1: 	return novelty, smoothed_novelty, novelty_peaks
mitian@1: 	
mitian@1: def plot_detection(smoothed_novelty, novelty_peaks):
mitian@17: 	pass