__author__ = 'Joshua Zosky'

"""
    Copyright 2015 Joshua Zosky
    joshua.e.zosky@gmail.com

    This file is part of "RetroTS".
    "RetroTS" is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    "RetroTS" is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with "RetroTS".  If not, see <http://www.gnu.org/licenses/>.
"""

from numpy import zeros, ones, nonzero, pi, argmin, sin, cos
from numpy import size, arange, clip, histogram, r_, Inf, divide, append, delete, array
from zscale import z_scale
import matplotlib.pyplot as plt


def my_hist(x, bin_centers):
    """
    This frivolous yet convenient conversion from bin-edges to bin-centers is from Stack Overflow user Bas Swinckels
    http://stackoverflow.com/questions/18065951/why-does-numpy-histogram-python-leave-off-one-element-as-compared-to-hist-in-m
    :param x:dataset
    :param bin_centers:bin values in a list to be moved from edges to centers
    :return: counts = the data in bin centers ready for pyplot.bar
    """
    bin_edges = r_[-Inf, 0.5 * (bin_centers[:-1] + bin_centers[1:]), Inf]
    counts, edges = histogram(x, bin_edges)
    return counts

def phase_base(amp_type, phasee):
    """
    

    :type phasee: object
    :param amp_type:    if 0, it is a time-based phase estimation
                        if 1, it is an amplitude-based phase estimation
    :param phasee: phasee information
    :return:
    """
    if amp_type == 0:
        # Calculate the phase of the trace, with the peak to be the start of the phase
        nptrc = len(phasee['tp_trace'])
        phasee['phase'] = -2 * ones(size(phasee['t']))
        i = 0
        j = 0
        while i <= (nptrc - 2):
            while phasee['t'][j] < phasee['tp_trace'][i + 1]:
                if phasee['t'][j] >= phasee['tp_trace'][i]:
                    # Note: Using a constant 244 period for each interval
                    # causes slope discontinuity within a period.
                    # One should resample period[i] so that it is
                    # estimated at each time in phasee['t'][j],
                    # dunno if that makes much of a difference in the end however.
                    if j == 10975:
                        pass
                    phasee['phase'][j] = (phasee['t'][j] - phasee['tp_trace'][i]) / phasee['prd'][i]\
                                         + phasee['zero_phase_offset']
                    if phasee['phase'][j] < 0:
                        phasee['phase'][j] = -phasee['phase'][j]
                    if phasee['phase'][j] > 1:
                        phasee['phase'][j] -= 1
                j += 1
            if i == 124:
                pass
            i += 1

        # Remove the points flagged as unset
        temp = nonzero(phasee['phase'] < -1)
        phasee['phase'][temp] = 0.0
        # Change phase to radians
        phasee['phase'] = phasee['phase'] * 2 * pi
    else:  # phase based on amplitude
        # at first scale to the max
        mxamp = max(phasee['p_trace'])
        phasee['phase_pol'] = []
        gR = z_scale(phasee['v'], 0, mxamp) # Scale, per Glover 2000's paper
        bins = arange(.01, 1.01, .01) * mxamp
        hb_value = my_hist(gR,bins)
        #hb_value = histogram(gR, bins)
        if phasee['show_graphs'] == 1:
            center = (bins[:-1] + bins[1:]) / 2
            plt.bar(center, hb_value[:len(hb_value)-1])#, align='center')
            plt.show()
        #find the polarity of each time point in v
        i = 0
        itp = 0
        inp = 0
        tp = phasee['tp_trace'][0]
        tn = phasee['tn_trace'][0]
        while (i <= len(phasee['v'])) and (phasee['t'][i] < tp) and (phasee['t'][i] < tn):
            phasee['phase_pol'].append(0)
            i += 1
        if tp < tn:
            # Expiring phase (peak behind us)
            cpol = -1
            itp = 1
        else:
            # Inspiring phase (bottom behind us)
            cpol = 1
            inp = 1
        phasee['phase_pol'] = zeros(size(phasee['v']))  # Not sure why you would replace the
                                                        # list that you created 10 lines prior to this
        # Add a fake point to tptrace and tntrace to avoid ugly if statements
        phasee['tp_trace'] = append(phasee['tp_trace'], phasee['t'][-1])
        phasee['tn_trace'] = append(phasee['tn_trace'], phasee['t'][-1])
        while i < len(phasee['v']):
            phasee['phase_pol'][i] = cpol
            if phasee['t'][i] == phasee['tp_trace'][itp]:
                cpol = -1
                itp = min((itp + 1), (len(phasee['tp_trace']) - 1))
            elif phasee['t'][i] == phasee['tn_trace'][inp]:
                cpol = 1
                inp = min((inp + 1), (len(phasee['tn_trace']) - 1))
            # cpol, inp, itp, i, R
            i += 1
        phasee['tp_trace'] = delete(phasee['tp_trace'], -1)
        phasee['tn_trace'] = delete(phasee['tn_trace'], -1)
        if phasee['show_graphs'] == 1:
            #clf
            plt.plot(phasee['t'], gR,'b')
            ipositive = nonzero(phasee['phase_pol'] > 0)
            ipositive = ipositive[0]
            ipositive_x = []
            for i in ipositive:
                ipositive_x.append(phasee['t'][i])
            ipositive_y = zeros(size(ipositive_x))
            ipositive_y.fill(0.55 * mxamp)
            plt.plot(ipositive_x, ipositive_y, 'r.')
            inegative = nonzero(phasee['phase_pol'] < 0)
            inegative = inegative[0]
            inegative_x = []
            for i in inegative:
                inegative_x.append(phasee['t'][i])
            inegative_y = zeros(size(inegative_x))
            inegative_y.fill(0.45 * mxamp)
            plt.plot(inegative_x, inegative_y, 'g.')
            plt.show()
        # Now that we have the polarity, without computing sign(dR/dt)
        #   as in Glover et al 2000, calculate the phase per eq. 3 of that paper
        # First the sum in the numerator
        for i, val in enumerate(gR):
            gR[i] = (round(val / mxamp * 100) + 1)
        gR = clip(gR, 0, 99)
        shb = sum(hb_value)
        hbsum = []
        hbsum.append(float(hb_value[0]) / shb)
        for i in range(1, 100):
            hbsum.append(hbsum[i - 1] + (float(hb_value[i]) / shb))
        for i in range(len(phasee['t'])):
            phasee['phase'].append(pi * hbsum[int(gR[i]) - 1] * phasee['phase_pol'][i])
        phasee['phase'] = array(phasee['phase'])

    # Time series time vector
    phasee['time_series_time'] = arange(0, (max(phasee['t']) - 0.5 * phasee['volume_tr']), phasee['volume_tr'])
    # Python uses half open ranges, so we need to catch the case when the stop
    # is evenly divisible by the step and add one more to the time series in
    # order to match Matlab, which uses closed ranges  1 Jun 2017 [D Nielson]
    if (max(phasee['t']) - 0.5 * phasee['volume_tr'])%phasee['volume_tr'] == 0:
        phasee['time_series_time'] = append(phasee['time_series_time'],[phasee['time_series_time'][-1]+phasee['volume_tr']])
    phasee['phase_slice'] = zeros((len(phasee['time_series_time']), phasee['number_of_slices']))
    phasee['phase_slice_reg'] = zeros((len(phasee['time_series_time']), 4, phasee['number_of_slices']))
    for i_slice in range(phasee['number_of_slices']):
        tslc = phasee['time_series_time'] + phasee['slice_offset'][i_slice]
        for i in range(len(phasee['time_series_time'])):
            imin = argmin(abs(tslc[i] - phasee['t']))
            #mi = abs(tslc[i] - phasee['t']) # probably not needed
            phasee['phase_slice'][i, i_slice] = phasee['phase'][imin]
        # Make regressors for each slice
        phasee['phase_slice_reg'][:, 0, i_slice] = sin(phasee['phase_slice'][:,i_slice])
        phasee['phase_slice_reg'][:, 1, i_slice] = cos(phasee['phase_slice'][:,i_slice])
        phasee['phase_slice_reg'][:, 2, i_slice] = sin(2 * phasee['phase_slice'][:,i_slice])
        phasee['phase_slice_reg'][:, 3, i_slice] = cos(2 * phasee['phase_slice'][:,i_slice])

    if phasee['quiet'] == 0 and phasee['show_graphs'] == 1:
        print '--> Calculated phase'
        plt.subplot(413)
        a = divide(divide(phasee['phase'], 2), pi)
        plt.plot(phasee['t'], divide(divide(phasee['phase'], 2), pi), 'm')
        if 'phase_r' in phasee:
            plt.plot(phasee['tR'], divide(divide(phasee['phase_r'], 2), pi), 'm-.')
        plt.subplot(414)
        plt.plot(phasee['time_series_time'], phasee['phase_slice'][:,1], 'ro',
                phasee['time_series_time'], phasee['phase_slice'][:,2], 'bo',
                phasee['time_series_time'], phasee['phase_slice'][:,2], 'b-')
        plt.plot(phasee['t'], phasee['phase'], 'k')
        #grid on
        # title it
        plt.title(phasee['v_name'])
        plt.show()
        # Need to implement this yet
        #if phasee['Demo']:
            #uiwait(msgbox('Press button to resume', 'Pausing', 'modal'))
    return phasee


def phase_estimator(amp_phase, phase_info):
    '''
    v_name='',
    amp_phase=0,
    t=[],
    x=[],
    iz=[],   # zero crossing (peak) locations
    p_trace=[],
    tp_trace=[],
    n_trace=[],
    tn_trace=[],
    prd=[],
    t_mid_prd=[],
    p_trace_mid_prd=[],
    phase=[],
    rv=[],
    rvt=[],
    var_vector=[],
    phys_fs=(1 / 0.025),
    zero_phase_offset=0.5,
    quiet=0,
    resample_fs=(1 / 0.025),
    f_cutoff=10,
    fir_order=80,
    resample_kernel='linear',
    demo=0,
    as_window_width=0,
    as_percover=0,
    as_fftwin=0,
    sep_dups=0,
    phasee_list=0,
    show_graphs=0
    '''
    """
    Example: PhaseEstimator.phase_estimator(amp_phase, info_dictionary)
    or PhaseEstimator.phase_estimator(v) where v is a list
    if v is a matrix, each column is processed separately.
    :param var_vector: column vector--list of list(s)
    :param phys_fs: Sampling frequency
    :param zero_phase_offset: Fraction of the period that corresponds to a phase of 0
                                0.5 means the middle of the period, 0 means the 1st peak
    :param quiet:
    :param resample_fs:
    :param frequency_cutoff:
    :param fir_order: BC ???
    :param resample_kernel:
    :param demo:
    :param as_window_width:
    :param as_percover:
    :param fftwin:
    :param sep_dups:
    :return: *_phased: phase estimation of input signal
    """
    phasee = dict(v_name='',
                  t=[],
                  x=[],
                  iz=[],   # zero crossing (peak) locations
                  volume_tr=2,
                  p_trace=[],
                  tp_trace=[],
                  n_trace=[],
                  tn_trace=[],
                  prd=[],
                  t_mid_prd=[],
                  p_trace_mid_prd=[],
                  phase=[],
                  rv=[],
                  rvt=[],
                  var_vector=[],
                  phys_fs=(1 / 0.025),
                  zero_phase_offset=0.5,
                  quiet=0,
                  resample_fs=(1 / 0.025),
                  frequency_cutoff=10,
                  fir_order=80,
                  resample_kernel='linear',
                  demo=0,
                  as_window_width=0,
                  as_percover=0,
                  as_fftwin=0,
                  sep_dups=0,
                  phasee_list=0,
                  show_graphs=0,
                  number_of_slices=0)
    phasee.update(phase_info)
    if isinstance(phasee['phasee_list'], type([])):
        return_phase_list = []
        for phasee_column in phasee['phasee_list']:
            return_phase.append(phase_base(amp_phase, phasee_column))
        return return_phase_list
    else:
        return_phase = phase_base(amp_phase, phasee)
        return return_phase