# -*- coding: utf-8 -*-
"""
Created on Wed Jun 17 14:19:08 2020

@author: kblackw1
"""

import numpy as np
single_linear=True

#This bit of code tests whether ppERK increases linearly with duration or concentration
from scipy import optimize
import statsmodels.api as sm
from scipy.stats import pearsonr
import pylab as plt
plt.ion()


def func_linear(x,m,b):                                     
    return m*x+b
def func_hill(x,pow,Kd):
    return x**pow/(x**pow+Kd**pow)
def show_curvefit(X,Y,params,fname):
    plt.figure()
    plt.plot(X,Y/np.max(Y),'*',label='data')
    estimate=func_hill(X,params[0],params[1])
    plt.plot(X,estimate,label='estimate')
    plt.legend()
    plt.title(fname)
    plt.text(0.1,0.1,' '.join([str(round(p,2)) for p in params]),horizontalalignment='center')
                  
def ols_fit(X,Y,constant = None,printR = None): #ordinary least squares
    if constant:                   #by default no constant or y-intercept fitted to equation
        X = sm.add_constant(X)      
    model = sm.OLS(Y,X).fit()
    if printR:
        print(model.summary())     #by default no output printed within function
    return model,model.predict(X)

def regress_val(model, split = None, r = None, constant=False): 
    values={}
    #values['fstat']=round(model.fvalue,2)
    values['adjR']=round(model.rsquared_adj,4)
    #values['f_pvalue']=model.f_pvalue
    '''
    if constant==True:
        values['slope']=model.params[1]
        values['tpvalue']=model.pvalues[1] #significant of the slope
    else:
        values['slope']=model.params
        values['tpvalue']=model.pvalues
    '''
    values['aic']=model.aic
    return values

def fit_hill(dur,ppERK, add_constant=False):
    ppERKnorm=[erk/np.max(ppERK) for erk in ppERK]
    #print(ppERKnorm)
    popt,pcov=optimize.curve_fit(func_hill,dur,ppERKnorm,bounds=([0,0],[10,10000]))
    pow=popt[0]
    Kd=popt[1]
    #print('***********************''Kd',Kd,'*********************','pow',pow)
    hill=[d**pow/(d**pow+Kd**pow) for d in dur]
    model,predictions = ols_fit(hill,ppERKnorm,constant=add_constant)
    return regress_val(model),popt

def fit_log(dur,ppERK,add_constant=False):
    logdur=[np.log10(d) for d in dur]
    model,predictions = ols_fit(logdur,ppERK,constant=add_constant)
    return regress_val(model)


if single_linear==True:
    import glob
    path='./'
    pattern=path+'*_nonLTP_lineartest.txt*'
    add_constant=True
    files=glob.glob(pattern)
    regress={}
    for filename in files:
        fname=filename.split('/')[-1].split('.')[0]
        types=fname.split('_')[1]
        with open(fname+'.txt') as f:
            data = []
            for line in f.readlines():
                data.append(line.split())
            if not data[1][1].replace('.','',1).isdigit():
                types='Dur'
                dur=[float(dat[0]) for dat in data[1:]]
                trials=[dat[1] for dat in data[1:]]
                ppERK=[float(dat[2]) for dat in data[1:]]
            else:
                conc=[float(dat[0]) for dat in data[1:]]
                dur=[float(dat[1]) for dat in data[1:]]
                trials=[dat[2] for dat in data[1:]]
                ppERK=[float(dat[3]) for dat in data[1:]]
        key=fname.split('-')[0]
        if types=="Dur":
            model,predictions = ols_fit(dur,ppERK,constant=add_constant)
            print(fname,'PR',pearsonr(dur,ppERK),'PR2',round(pearsonr(dur,ppERK)[0]**2,3))
        else:
            model,predictions = ols_fit(conc,ppERK,constant=add_constant)
            print(fname,'PR',pearsonr(conc,ppERK),'PR2',round(pearsonr(conc,ppERK)[0]**2,3))
        regress[key] = regress_val(model, constant=add_constant)

        if types=="Dur":
            regress[key+'_log']=fit_log(dur,ppERK)
            regress[key+'_hill'],hill_coef=fit_hill(dur,ppERK)
            show_curvefit(dur,ppERK,hill_coef,fname)
        else:
            regress[key+'_log']=fit_log(conc,ppERK)
            regress[key+'_hill'],hill_coef=fit_hill(conc,ppERK)
            show_curvefit(conc,ppERK,hill_coef,fname)
    for k,vals in regress.items():
        print(k,vals)

    with open('linear.txt', 'w') as f:
        for k,vals in regress.items():
            f.write(k+'  '+'  '.join([str(a)+'='+str(b) for a,b in vals.items()])+'\n')
  
    '''
        popt, pcov = optimize.curve_fit(func_linear,dur,ppERK)    #fits equation to choosen function if possible
        print('This is fit estimate for : ', fname)
        print('popt',popt)
        print('pcov',np.sqrt(np.diag(pcov))) 
    '''

##########This bit of code tests whether ppERK combination = summation 
else:
    import pandas as pd
    from statsmodels.formula.api import ols
    import os
    import matplotlib.pyplot as plt
    import glob
    from contextlib import redirect_stdout
    #fname='cAMPGbgCa_lineartest.txt'
    results={}
    PATH='./'
    Pattern=PATH+'*Combo'+"*-lineartest.txt"
    files=sorted(glob.glob(Pattern))
    #files=['allCombo_lineartest']

    for filename in files:
        fname=filename.split('/')[-1].split('.')[0]
        with open(fname+'.txt') as f:
            data = []
            for line in f.readlines():
                data.append(line.split())
                dur_or_Conc_or_ITI=[float(dat[0]) for dat in data[1:]]
                trials=[dat[1] for dat in data[1:]]
                input1_ppERK=np.array([float(dat[2]) for dat in data[1:]])
                input2_ppERK=np.array([float(dat[3]) for dat in data[1:]])
                if len(data[0])>5:
                    input3_ppERK=np.array([float(dat[4]) for dat in data[1:]])
                    combo_ppERK=np.array([float(dat[5]) for dat in data[1:]])
                    sum_ppERK=input1_ppERK+input2_ppERK+input3_ppERK
                else:
                  combo_ppERK=np.array([float(dat[4]) for dat in data[1:]])
                  sum_ppERK=input1_ppERK+input2_ppERK   
                #reformat the data
        sum_or_combo=['sum' for i in sum_ppERK]+['combo' for i in sum_ppERK]
        alltrials=trials+trials
        all_dur_conc=dur_or_Conc_or_ITI+dur_or_Conc_or_ITI
        ppERK=np.zeros(2*len(sum_ppERK))
        ppERK[0:len(sum_ppERK)]=sum_ppERK
        ppERK[len(sum_ppERK):]=combo_ppERK
        df_data={'sum_combo':sum_or_combo,'trial':alltrials,'dur_conc_ITI':all_dur_conc,'ppERK':ppERK}
        df=pd.DataFrame(df_data)
        results=ols('ppERK ~C(sum_or_combo)+ dur_conc_ITI',data=df).fit()#ols('ppERK ~C(sum_or_combo)* dur_conc_ITI',data=df).fit()
        #print(fname,results.summary())
        plt.ion()
        plt.figure()
        plt.plot(dur_or_Conc_or_ITI,sum_ppERK,'ro',label=fname+'sum')
        plt.plot(dur_or_Conc_or_ITI,combo_ppERK,'bs',label=fname+'combo')
        plt.legend()
        with open(fname.split('-')[0]+'-NOVA.txt','w') as f:
            with redirect_stdout(f):
                print(fname,results.summary())
        


