# -*- coding: utf-8 -*-
"""
Created on Sat Apr  5 15:28:13 2014

@author: deokwooj
"""
from __future__ import division # To forace float point division
import numpy as np
from scipy import stats
from scipy.interpolate import interp1d
from sklearn import mixture
#from sklearn.cluster import Ward
from sklearn.cluster import KMeans
import time
##################################################################
# Custom library
##################################################################
from data_tools import *
from shared_constants import *
import pprint
import lib_bnlearn as rbn    

def X_INPUT_to_states(xinput,CORR_VAL_OUT=0,PARALLEL = False):
    #import pdb;pdb.set_trace()
    sinput=np.zeros(xinput.shape)
    num_samples=xinput.shape[0]
    num_sensors=xinput.shape[1]
    if num_samples <num_sensors:
        print 'Warning number of samplesa are smaller than number of sensors'

    print 'Mapping',xinput.shape, ' marix to discrete states '
    for k,samples in enumerate(xinput.T):
        obs=samples[:,np.newaxis]
        label,opt_num_cluster,model,score,score_err_sum= state_retrieval(obs,max_num_cluster=6,est_method='kmean',PARALLEL=PARALLEL,PRINTSHOW=False)
        high_peak_label_idx=np.argmax(model.cluster_centers_)
        low_peak_label_idx=np.argmin(model.cluster_centers_)
        high_peak_idx=np.nonzero(label==high_peak_label_idx)[0]
        sinput[high_peak_idx,k]=1
        low_peak_idx=np.nonzero(label==low_peak_label_idx)[0]
        sinput[low_peak_idx,k]=-1
    
    corr_state_val=[]    
    if CORR_VAL_OUT==1:
        print 'Compute Correlation Score....'
        for k,(row1, row2) in enumerate(zip(sinput.T, xinput.T)):
            corr_state_val.append(round(stats.pearsonr(row1,row2)[0],3))
    corr_state_val=np.array(corr_state_val)
    return sinput,corr_state_val
    
def interpolation_measurement(data_dict,input_names,err_rate=1,sgm_bnd=20):
    print 'interploattion starts....'
    measurement_point_set=[]
    num_of_discrete_val=[]
    sampling_interval_set=[]
    num_type_set=[]
    err_rate=1;sgm_bnd=20
    """
    try:
        import pdb;pdb.set_trace()
    except ValueError:
        import pdb;pdb.set_trace()
    """
    for i,key_name in enumerate(input_names):
        print key_name,'.....'
        t_=np.array(data_dict[key_name][2][0])
        
        if len(t_) == 0:
            continue
        
        intpl_intv=np.ceil((t_[-1]-t_[0]) /len(t_))
        sampling_interval_set.append(intpl_intv)
        val_=np.array(data_dict[key_name][2][1])
        num_of_discrete_val_temp=len(set(val_))
        num_of_discrete_val.append(num_of_discrete_val_temp)
        # filtering outlier
        # assuming 1% of errors and 30 x standard deviation rules
        outlier_idx=outlier_detect(val_,err_rate,sgm_bnd)
        if len(outlier_idx)>0:
            print 'outlier samples are detected: ', 'outlier_idx:', outlier_idx
            t_=np.delete(t_,outlier_idx)
            val_=np.delete(val_,outlier_idx)
        t_new=np.r_[t_[0]:t_[-1]:intpl_intv]
        """        
        if num_of_discrete_val_temp<MIN_NUM_VAL_FOR_FLOAT:
            num_type=INT_TYPE
            val_new=fast_nearest_interp(t_new, t_,val_)
        else:
            num_type=FLOAT_TYPE
            val_new = np.interp(t_new, t_,val_)
        """
        num_type=check_data_type(data_dict[key_name][2][1])
        if num_type==INT_TYPE:
            val_new=fast_nearest_interp(t_new, t_,val_)
        else:
            #num_type=FLOAT_TYPE
            val_new = np.interp(t_new, t_,val_)
        
        c=np.vstack([t_new,val_new])
        measurement_point_set.append(c)
        num_type_set.append(num_type)
        print '-----------------------------------------------------------------'
    #return measurement_point_set,num_type_set,num_of_discrete_val,sampling_interval_set
    return measurement_point_set,np.array(num_type_set)


def rolling_window(a, window):
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def rolling_window_label_mode(label,r_window):
    if (r_window/2)==int(r_window/2):
        r_window=int(r_window+1)
        #raise NameError('length of window size must be odd')
    offset=int(r_window/2)
    rw_label_temp=stats.mode(rolling_window(label, r_window),1)[0]
    head= rw_label_temp[0]*np.ones([offset,1])
    body=rw_label_temp
    tail= rw_label_temp[-1]*np.ones([offset,1])
    rw_label=np.r_[head,body,tail]
    return rw_label

def rolling_window_label_binary(label,r_window):
    if (r_window/2)==int(r_window/2):
        r_window=int(r_window+1)
        #raise NameError('length of window size must be odd')
    offset=int(r_window/2)
    rw_label_temp=np.array([ np.sum(temp)/r_window for temp in rolling_window(label, r_window)])
    #import pdb;pdb.set_trace()
    # rw_label_temp=stats.mode(rolling
_window(label, r_window),1)[0]
    head= rw_label_temp[0]*np.ones([offset,1])
    body=rw_label_temp
    tail= rw_label_temp[-1]*np.ones([offset,1])
    rw_label=np.r_[head,body[:,np.newaxis],tail]
    return rw_label

"""


def state_retrieval(obs,max_num_cluster=6,est_method='kmean'):
    #print '========================================================================='    
    #print 'Retrieving discrete states from data using ',est_method, ' model...'
    #print '=====================================
===================================='
    score=np.zeros(max_num_cluster)
    model_set=[]
    #print 'try ',max_num_cluster, ' clusters..... '
    for num_cluster in range(max_num_cluster):
        #print 'Try ',num_cluster+1, ' clusters '
        #print '-----------------------------------'
        if est_method=='kmean':
            kmean=KMeans(n_clusters=num_cluster+1).fit(obs)
            model_set.append(kmean)
            #import pdb;pdb.set_trace()
            score[num_cluster]=np.sum(kmean.score(obs))
        elif est_method=='gmm':
            gmm = mixture.GMM(n_components=num_cluster+1).fit(obs) 
            model_set.append(gmm)
            score[num_cluster]=np.sum(gmm.score(obs))
        else:
            raise NameError('not supported est_method')
    score_err_sum=np.zeros(max_num_cluster)
    #print 'Finding knee points of log likelihood...'
    for i in range(max_num_cluster):
        a_0=score[:(i)]
        if len(a_0)>1:
            slope, intercept, r_value, p_value, std_err = stats.linregress(range(len(a_0)),a_0)
            sqr_sum_err0=sum(((slope*np.arange(len(a_0))+ intercept)-a_0)**2)
        else:
            sqr_sum_err0=0
        a_1=score[(i):]
        if len(a_1)>1:
            slope, intercept, r_value, p_value, std_err = stats.linregress(range(len(a_1)),a_1)
            sqr_sum_err1=sum(((slope*np.arange(len(a_1))+ intercept)-a_1)**2)
        else:
            sqr_sum_err1=0
        score_err_sum[i]=sqr_sum_err0+sqr_sum_err1
    # Optimum number of clusters. 
    min_idx=np.argmin(score_err_sum)
    opt_num_cluster=min_idx+1
    #print 'opt_num_cluster: ' , opt_num_cluster
    if est_method=='kmean':
        label=model_set[min_idx].labels_
    elif est_method=='gmm':
        label=model_set[min_idx].predict(obs)
    else:
        raise NameError('not supported est_method')
    return label,opt_num_cluster, model_set[min_idx],score,score_err_sum

"""


def cluster_state_retrieval(tup):
    obs = tup[0]
    num_clusters = tup[1]
    est_method = tup[2]
    #print 'num clusters = ' + str(num_clusters)
    if est_method=='kmean':
        kmean=KMeans(n_clusters=num_clusters).fit(obs)
        model = kmean
        score=compute_log_ll(kmean.labels_,obs)
        #score=-1*np.log(-1*np.sum(kmean.score(obs)))
    elif est_method=='gmm':
        gmm = mixture.GMM(n_components=num_clusters).fit(obs)
        model = gmm
        score=np.sum(gmm.score(obs))
    #print 'Done ' + str(num_clusters)        
    return (num_clusters-1, [model,score])
from multiprocessing import Pool    


def compute_log_ll(label_in,obs_in):
    log_ll_sum=0
    for i in range(label_in.max()+1):
        idx=np.nonzero(label_in==i)[0]
        val_set=obs_in[idx]
        log_val=stats.n