'''
This program:
- Reads PPG data from a csv file and produces a scatterplot
- Generates a linear combination f(t) = A1*sin(w1*t+c1) + A2*cos(w2*t+c2) + h
  to match the PPG data points
- Graph data points and function
'''

import csv
import matplotlib.pyplot as plt
import numpy as np
import math
from matplotlib.widgets import Slider,RadioButtons 
from math import pi

#Read PPG data in csv file 
def read_csv(filename):
    global h_init, h_max, h_min
    global a1, a2, freq
    times=[]
    amps=[]
    xmax=[]
    xmin=[]
    n = 0
    #Read PPG data from csv file
    with open(filename) as f:
        reader=csv.reader(f)
        next(reader)
        for row in reader:
            amps.append(float(row[0]))
            times.append(float(row[2]))
            n = n + 1
    #Calculates harmonic function vertical shift
    h_init = sum(amps)/len(amps)
    #Calculates Initial amplitudes for Sine and Cosine  
    a1 = (max(amps) - min(amps))/2
    a2 = a1/2
    #Calculates upper and lower limits for function vertical shift
    h_max = int(h_init + a1)
    h_min = int(h_init - a1)

    #Begin process to determine the PPG data average frequency 
    a = amps[1] - amps[0]
    #Sweep points looking for max and mins
    for i in range(1,n-1):
        #Find the signs of the slopes between 3-contigous points 
        b = amps[i+1] - amps[i]
        sa = np.sign(a)
        sb = np.sign(b)
        sn = sa*sb
        #Sign change to determine max or min & to include them in a list
        if sn < 0 and sa > 0:
            xmax.append(times[i])
        if sn < 0 and sa < 0:
            xmin.append(times[i])
        #Test when contiguous points have same y-value
        if sn == 0:
            while sb == 0:
                i = i+1
                sb = np.sign(amps[i+1] - amps[i])
            sn = sa*sb
            if sn < 0 and sa > 0:
                xmax.append(times[i])
            if sn < 0 and sa < 0:
                xmin.append(times[i])
        a = b       

    dif_max = []
    dif_min = []
    #Finding distances between maximums and minimums
    nmax = len(xmax) - 1
    nmin = len(xmin) - 1
    for j in range(0,nmax):
        dif_max.append(xmax[j+1]-xmax[j])
    for k in range(0,nmin):
        dif_min.append(xmin[k+1]-xmin[k])

    #computing average period for PPG data
    period_max = sum(dif_max)/len(dif_max)
    period_min = sum(dif_min)/len(dif_min)
    period = (period_max + period_min)/2
    #computing average frequency for PPG data
    freq = 1/period

    return times,amps

#Generate points to graph f(t) = A1*sin(w1*t+c1) + A2*cos(w2*t+c2) + h
def generate_f():
    global xr
    #Generate function x-values
    xr=np.linspace(25, 31, 300)
    #Initialize Lists and function parameters
    xi=[]
    fxi=[]

    for n in xr:
        xn=n
        xi.append(xn)
        fn=a1*np.sin(2*pi*freq*xn+c1)+a2*np.cos(4*pi*freq*xn+c2)+h_init
        fxi.append(fn)


    return xi,fxi

#Display Graphs and labels graphing window
def display_graphs():
    global fun_plot, fig
    global vshift_slider
    global amp_slider_sin, freq_slider_sin, hshift_slider_sin
    global amp_slider_cos, freq_slider_cos, hshift_slider_cos
    
    #Define size in inches of the whole window
    fig = plt.figure(figsize=(16,8),facecolor='y')

    #Define size & position of graphing window
    graphs = plt.axes([0.05, 0.25, 0.925, 0.7])

    #Define positions for sliders
    sl_vshift = plt.axes([0.385, 0.145, 0.25, 0.025])
    sl_amp_sin = plt.axes([0.0725, 0.095, 0.375, 0.025])
    sl_freq_sin = plt.axes([0.0725, 0.055, 0.375, 0.025])
    sl_hshift_sin = plt.axes([0.0725, 0.015, 0.375, 0.025])
    sl_amp_cos = plt.axes([0.575, 0.095, 0.375, 0.025])
    sl_freq_cos = plt.axes([0.575, 0.055, 0.375, 0.025])
    sl_hshift_cos = plt.axes([0.575, 0.015, 0.375, 0.025])

    #Define the sliders
    vshift_slider = Slider(sl_vshift,'v_shift',h_min,h_max,valinit=h_init)
    #Correct this 40
    amp_slider_sin = Slider(sl_amp_sin,'Amp_Sin',0,1.5*a1,valinit=a1,facecolor='r')
    freq_slider_sin = Slider(sl_freq_sin,'freq_Sin',0,2*freq,valinit=freq,facecolor='r')
    hshift_slider_sin = Slider(sl_hshift_sin,'Shift_Sin',-2*pi,2*pi,valinit=c1,facecolor='r')
    #correct this 40
    amp_slider_cos = Slider(sl_amp_cos,'Amp_Cos',0,1.5*a2,valinit=a2,facecolor='g')
    freq_slider_cos = Slider(sl_freq_cos,'freq_Cos',0,4*freq,valinit=2*freq,facecolor='g')
    hshift_slider_cos = Slider(sl_hshift_cos,'Shift_Cos',-2*pi,2*pi,valinit=c2,facecolor='g')

    #Define size & position of graphing window
    graphs = plt.axes([0.05, 0.25, 0.925, 0.7])

    #Set labels for graphing window
    plt.xlabel('Time(seconds)')
    plt.ylabel('Amplitude(mV)')
    plt.title('Photoplethysmography Data & Harmonic Function')    
    #fun_plot.set_ydata(a1*np.sin(w1*pi*xr+c1)+a2*np.cos(w2*pi*xr+c2)+h)
    #Plot PPG points and Intial Harmonic Function
    plt.plot(times,amps,'o')
    fun_plot, = plt.plot(xi,fxi)

    #Modify Harmonic Functions parameters with Sliders
    vshift_slider.on_changed(update)
    amp_slider_sin.on_changed(update)
    freq_slider_sin.on_changed(update)
    hshift_slider_sin.on_changed(update)
    amp_slider_cos.on_changed(update)
    freq_slider_cos.on_changed(update)
    hshift_slider_cos.on_changed(update)    

    plt.show()

#Update the harmonic function with the slliders values 
def update(val):
    global w1,c1,w2,c2
    h=vshift_slider.val
    a1=amp_slider_sin.val
    w1=freq_slider_sin.val
    c1=hshift_slider_sin.val
    a2=amp_slider_cos.val
    w2=freq_slider_cos.val
    c2=hshift_slider_cos.val    

    #Update harmonic function
    fun_plot.set_ydata(a1*np.sin(2*pi*w1*xr+c1)+a2*np.cos(2*pi*w2*xr+c2)+h)
    fig.canvas.draw_idle()


###########################################################################
if __name__=='__main__':

    c1=0
    c2=0.5

    times, amps = read_csv('ppg_short.csv')
    xi,fxi = generate_f()
    display_graphs()