import numpy as np
import sys
sys.path.append('../')
from transfer_functions.load_config import load_transfer_functions
from scipy.integrate import odeint


def build_up_differential_operator_first_order(TF1, TF2, T=5e-3):
    """
    simple first order system
    """
    def A0(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)-V[0])
    
    def A1(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(TF2(V[0]+exc_aff, V[1]+inh_aff)-V[1])
    
    def Diff_OP(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return np.array([A0(V, exc_aff=exc_aff,inh_aff=inh_aff, pure_exc_aff=pure_exc_aff),\
                         A1(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff)])
    return Diff_OP
    

def build_up_differential_operator(TF1, TF2,\
                                   Ne=8000, Ni=2000, T=5e-3):
    """
    Implements Equation (3.16) in El BOustani & Destexhe 2009
    in the case of a network of two populations:
    one excitatory and one inhibitory
    
    Each neuronal population has the same transfer function
    this 2 order formalism computes the impact of finite size effects
    T : is the bin for the Markovian formalism to apply

    the time dependent vector vector is V=[fe,fi, sfe, sfi, sfefi]
    the function returns Diff_OP
    and d(V)/dt = Diff_OP(V)
    """
    
    # we have the transfer function, now we also get its derivatives
    # TF, diff_fe, diff_fi, diff2_fe_fe, diff2_fe_fi, diff2_fi_fi, values = \
    #                         get_derivatives_of_TF(params)
    
    def A0(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(\
                .5*V[2]*diff2_fe_fe(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                .5*V[3]*diff2_fe_fi(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                .5*V[3]*diff2_fi_fe(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                .5*V[4]*diff2_fi_fi(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)-V[0])
    
    def A1(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(\
                .5*V[2]*diff2_fe_fe(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                .5*V[3]*diff2_fe_fi(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                .5*V[3]*diff2_fi_fe(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                .5*V[4]*diff2_fi_fi(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                TF2(V[0]+exc_aff, V[1]+inh_aff)-V[1])
    
    def A2(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(\
                1./Ne*TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)*(1./T-TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff))+\
                (TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)-V[0])**2+\
                2.*V[2]*diff_fe(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                2.*V[3]*diff_fi(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                -2.*V[2])
    
    def A3(V, exc_aff=0, inh_aff=0, pure_exc_aff=0): # mu, nu = e,i, then lbd = e then i
        return 1./T*(\
               (TF1(V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)-V[0])*(TF2(V[0]+exc_aff, V[1]+inh_aff)-V[1])+\
                V[2]*diff_fe(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                V[3]*diff_fe(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                V[3]*diff_fi(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                V[4]*diff_fi(TF1, V[0]+exc_aff+pure_exc_aff, V[1]+inh_aff)+\
                -2.*V[3])
    
    def A4(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return 1./T*(\
                1./Ni*TF2(V[0]+exc_aff, V[1]+inh_aff)*(1./T-TF2(V[0]+exc_aff, V[1]+inh_aff))+\
                (TF2(V[0]+exc_aff, V[1]+inh_aff)-V[1])**2+\
                2.*V[3]*diff_fe(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                2.*V[4]*diff_fi(TF2, V[0]+exc_aff, V[1]+inh_aff)+\
                -2.*V[4])
    
    def Diff_OP(V, exc_aff=0, inh_aff=0, pure_exc_aff=0):
        return np.array([A0(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff),\
                         A1(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff),\
                         A2(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff),\
                         A3(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff),\
                         A4(V, exc_aff=exc_aff, inh_aff=inh_aff, pure_exc_aff=pure_exc_aff)])
    return Diff_OP

##### Derivatives taken numerically,
## to be implemented analitically ! not hard...

def diff_fe(TF, fe, fi, df=1e-4):
    return (TF(fe+df/2., fi)-TF(fe-df/2.,fi))/df

def diff_fi(TF, fe, fi, df=1e-4):
    return (TF(fe, fi+df/2.)-TF(fe, fi-df/2.))/df

def diff2_fe_fe(TF, fe, fi, df=1e-4):
    return (diff_fe(TF, fe+df/2., fi)-diff_fe(TF,fe-df/2.,fi))/df

def diff2_fi_fe(TF, fe, fi, df=1e-4):
    return (diff_fi(TF, fe+df/2., fi)-diff_fi(TF,fe-df/2.,fi))/df

def diff2_fe_fi(TF, fe, fi, df=1e-4):
    return (diff_fe(TF, fe, fi+df/2.)-diff_fe(TF,fe, fi-df/2.))/df

def diff2_fi_fi(TF, fe, fi, df=1e-4):
    return (diff_fi(TF, fe, fi+df/2.)-diff_fi(TF,fe, fi-df/2.))/df

def find_fixed_point_first_order(NRN1, NRN2, NTWK,\
                                 Ne=8000, Ni=2000, exc_aff=0.,\
                                 verbose=False):

    TF1, TF2 = load_transfer_functions(NRN1, NRN2, NTWK)
    
    t = np.arange(2000)*1e-4              # time
    
    ### FIRST ORDER ###
    def dX_dt_scalar(X, t=0):
        return build_up_differential_operator_first_order(TF1, TF2, T=5e-3)(X, exc_aff=exc_aff)
    X0 = [1, 10] # need inhibition stronger than excitation
    X = odeint(dX_dt_scalar, X0, t)         # we don't need infodict here
    if verbose:
        print('first order prediction: ', X[-1])
    return X[-1][0], X[-1][1] 

def find_fixed_point(NRN1, NRN2, NTWK, Ne=8000, Ni=2000, exc_aff=0., verbose=False):

    # we start from the first order prediction !!!
    X0 = find_fixed_point_first_order(NRN1, NRN2, NTWK,\
                                      Ne=Ne, Ni=Ni, exc_aff=exc_aff,\
                                      verbose=verbose)
    X0 = [X0[0], X0[1], .5, .5, .5]
    
    TF1, TF2 = load_transfer_functions(NRN1, NRN2, NTWK)
    t = np.arange(200)*1e-4              # time

    ### SECOND ORDER ###
    # def dX_dt_scalar(X, t=0):
    #     return build_up_differential_operator(TF1, TF2,\
    #                                           Ne=Ne, Ni=Ni)(X, exc_aff=exc_aff)
    # X = odeint(dX_dt_scalar, X0, t)         # we don't need infodict here

    # simple euler
    X = X0
    for i in range(len(t)-1):
        X = X + (t[1]-t[0])*build_up_differential_operator(TF1, TF2,Ne=Ne, Ni=Ni)(X, exc_aff=exc_aff)
        last_X = X

    print('Make sure that those two values are similar !!')
    print(X)
    print(last_X)
    
    if verbose:
        print(X)
    if verbose:
        print('first order prediction: ',X[-1])
    
    # return X[-1][0], X[-1][1], np.sqrt(X[-1][2]), np.sqrt(X[-1][3]), np.sqrt(X[-1][4])
    return X[0], X[1], np.sqrt(X[2]), np.sqrt(X[3]), np.sqrt(X[4])

if __name__=='__main__':

    # find_fixed_point('LIF', 'LIF', 'Vogels-Abbott', exc_aff=0., Ne=4000, Ni=1000, verbose=True)
    find_fixed_point('RS-cell', 'FS-cell', 'CONFIG1', exc_aff=4., Ne=8000, Ni=2000, verbose=True)
    

    
