"""\
Create general Channel Proto, pass in name, x and y power, and params

Also, create the library of channels
chan_proto requires alpha and beta params for both activation and inactivation
If channel does not have inactivation, just send in empty Yparam array.
"""

from __future__ import print_function, division
import moose
import numpy as np
import logging

from moose_nerp.prototypes import constants, logutil
from moose_nerp.prototypes.util import NamedList
log = logutil.Logger()

StandardMooseTauInfChannelParams = NamedList('StandardMooseTauInfChannelParams', '''
                                T_rate
                                T_B
                                T_C
                                T_vhalf
                                T_vslope
                                SS_rate
                                SS_B
                                SS_C
                                SS_vhalf
                                SS_vslope''')

TauInfMinChannelParams = NamedList('TauInfMinChannelParams', '''
                                T_min
                                T_vdep
                                T_vhalf
                                T_vslope
                                SS_min
                                SS_vdep
                                SS_vhalf
                                SS_vslope
                                T_power=1''')


AlphaBetaChannelParams = NamedList('AlphaBetaChannelParams', '''
                                A_rate
                                A_B
                                A_C
                                A_vhalf
                                A_vslope
                                B_rate
                                B_B
                                B_C
                                B_vhalf
                                B_vslope''')

ZChannelParams = NamedList('ZChannelParams', 'Kd power tau taumax=0 tau_power=0 cahalf=0')
BKChannelParams=NamedList('BKChannelParams', 'alphabeta K delta')


ChannelSettings = NamedList('ChannelSettings', 'Xpow Ypow Zpow Erev name')

def dist_dependent_cond_equation(cmin, cmax, dhalf, slope):
    '''Returns a function only of distance that specifies a distance
    dependent sigmoid equation for channel conductance.

    cmin: minimal distance-dependent conductance
    cmax: maximal distant-dependent conductance
    dhalf: distance from soma where amplitude of sigmoid is one half max-min
    slope: slope of sigmoid equation. If positive, somatic conductance is higher; If negative, somatic conductance is smaller
    '''

    def equation(distance):
        return cmin+(cmax-cmin)/(1+np.exp((distance-dhalf)/slope))

    print('{}+({}-{})/(1+np.exp((distance-{})/{})))'.format(cmin,cmax,cmin,dhalf,slope))
    return equation

def sigmoid(x,xmin,xmax,xvhalf,xslope):
    return xmin+xmax/(1+np.exp((x-xvhalf)/xslope))
#notice the x-xvhalf in sigmoid, but x+xvhalf used by MOOSE
def quadratic(x,xmin,xmax,xvhalf,xslope):
    tau1 = xmax/(1+np.exp((x-xvhalf)/xslope))
    tau2 = 1/(1+np.exp((x-xvhalf)/-xslope))
    tau_x = xmin+tau1*tau2
    return tau_x

def make_sigmoid_gate(model,params,Gate):
    v = np.linspace(model.VMIN, model.VMAX, model.VDIVS)
    if params.T_power==2:
        print('making quadratic gate', Gate.path)
        tau = quadratic(v,params.T_min,params.T_vdep,params.T_vhalf,params.T_vslope)
    else:
        tau = sigmoid(v,params.T_min,params.T_vdep,params.T_vhalf,params.T_vslope)
    minf = sigmoid(v,params.SS_min,params.SS_vdep,params.SS_vhalf,params.SS_vslope)
    Gate.min = model.VMIN
    Gate.max = model.VMAX
    Gate.divs = model.VDIVS
    Gate.tableA = minf/tau
    Gate.tableB = 1/tau

def interpolate_values_in_table(model, tabA, V_0, l=40):
    '''This function interpolates values in the table
    around tabA[V_0]. '''
    V = np.linspace(model.VMIN, model.VMAX, len(tabA))
    idx =  abs(V-V_0).argmin()
    A_min = tabA[idx-l]
    V_min = V[idx-l]
    A_max = tabA[idx+l]
    V_max = V[idx+l]
    a = (A_max-A_min)/(V_max-V_min)
    b = A_max - a*V_max
    tabA[idx-l:idx+l] = V[idx-l:idx+l]*a+b
    return tabA

def  calc_V0(rate,B,C,vhalf,vslope,Params):
        delta = rate - B * vhalf
        if delta > 1e-10 or delta < -1e-10:
            log.warning("Fixing Singularities. Please verify constraint on Beta: A = B * vhalf {}", Params)
            rate= B * vhalf
        V_0 = vslope*np.log(-C)-vhalf
        return rate,V_0

def fix_singularities(model, Params, Gate):
    #This needs to be extended to work with standardMooseTauInfparams
    if Params.A_C < 0:
        Params.A_rate,V_0=calc_V0(Params.A_rate,Params.A_B,Params.A_C,Params.A_vhalf,Params.A_vslope, Params)
        if model.VMIN < V_0 < model.VMAX:
            #change values in tableA and tableB, because tableB contains sum of alpha and beta
            Gate.tableA = interpolate_values_in_table(model, Gate.tableA, V_0)
            Gate.tableB = interpolate_values_in_table(model, Gate.tableB, V_0)

    if Params.B_C < 0:
        Params.B_rate,V_0=calc_V0(Params.B_rate,Params.B_B,Params.B_C,Params.B_vhalf,Params.B_vslope, Params)
        if model.VMIN < V_0 < model.VMAX:
            #change values in tableB (alpha is stored in tableA)
            Gate.tableB = interpolate_values_in_table(model, Gate.tableB, V_0)

def make_gate(params,model,gate):
    if isinstance(params,AlphaBetaChannelParams):
        gate.setupAlpha(params + [model.VDIVS, model.VMIN, model.VMAX])
        fix_singularities(model, params, gate)
    elif isinstance(params,StandardMooseTauInfChannelParams):
        gate.setupTau(params + [model.VDIVS, model.VMIN, model.VMAX])
        fix_singularities(model, params, gate)
    elif isinstance(params,TauInfMinChannelParams):
        make_sigmoid_gate(model,params,gate)

def chan_proto(model, chanpath, params):
    log.info("{}: {}", chanpath, params)
    chan = moose.HHChannel(chanpath)

    chan.Xpower = params.channel.Xpow
    if params.channel.Xpow > 0:
        xGate = moose.HHGate(chan.path + '/gateX')
        make_gate(params.X,model,xGate)

    chan.Ypower = params.channel.Ypow
    if params.channel.Ypow > 0:
        yGate = moose.HHGate(chan.path + '/gateY')
        make_gate(params.Y,model,yGate)

    if params.channel.Zpow > 0:
        chan.Zpower = params.channel.Zpow
        zGate = moose.HHGate(chan.path + '/gateZ')
        if params.Z.__class__==ZChannelParams:
            #
            ca_array = np.linspace(model.CAMIN, model.CAMAX, model.CADIVS)
            zGate.min = model.CAMIN
            zGate.max = model.CAMAX
            caterm = (ca_array/params.Z.Kd) ** params.Z.power
            inf_z = caterm / (1 + caterm)
            if params.Z.taumax>0:
                tauterm=(ca_array/params.Z.cahalf)**params.Z.tau_power
                taumax_z=(params.Z.taumax-params.Z.tau)/(1+tauterm)
                taumin_z= params.Z.tau * np.ones(len(ca_array))
                tau_z = taumin_z+taumax_z
            else:
                tau_z = params.Z.tau * np.ones(len(ca_array))
            #
            zGate.tableA = inf_z / tau_z
            zGate.tableB = 1 / tau_z
            chan.useConcentration = True
        else:
            chan.useConcentration = False
            make_gate(params.Z,model,zGate)

    chan.Ek = params.channel.Erev
    chan.tick=-1
    return chan

def BKchan_proto(model, chanpath, params):
    ZFbyRT= 2 * constants.Faraday / (constants.R * constants.celsius_to_kelvin(model.Temp))
    v_array = np.linspace(model.VMIN, model.VMAX, model.VDIVS)
    ca_array = np.linspace(model.CAMIN, model.CAMAX, model.CADIVS)
    if model.VDIVS<=5 and model.CADIVS<=5:
        log.info("{}, {}", v_array, ca_array)
    gatingMatrix = []
    for i,pars in enumerate(params.X):
        Vdepgating=pars.K*np.exp(pars.delta*ZFbyRT*v_array)
        if i == 0:
            gatingMatrix.append(pars.alphabeta*ca_array[None,:]/(ca_array[None,:]+pars.K*Vdepgating[:,None]))
            #table.tableVector2D=gatingMatrix
        else:
            gatingMatrix.append(pars.alphabeta/(1+ca_array[None,:]/pars.K*Vdepgating[:,None]))
            gatingMatrix[i] += gatingMatrix[0]
            #table.tableVector2D=gatingMatrix


    chan = moose.HHChannel2D(chanpath)
    chan.Xpower = params.channel.Xpow
    chan.Ek=params.channel.Erev
    chan.Xindex="VOLT_C1_INDEX"
    xGate = moose.HHGate2D(chan.path + '/gateX')
    xGate.xminA=xGate.xminB=model.VMIN
    xGate.xmaxA=xGate.xmaxB=model.VMAX
    xGate.xdivsA=xGate.xdivsB=model.VDIVS
    xGate.yminA=xGate.yminB=model.CAMIN
    xGate.ymaxA=xGate.ymaxB=model.CAMAX
    xGate.ydivsA=xGate.ydivsB=model.CADIVS
    xGate.tableA=gatingMatrix[0]
    xGate.tableB=gatingMatrix[1]
    if log.isEnabledFor(logging.INFO):
        for ii in np.arange(0, model.VDIVS,1200):
            log.info("{} V={}", chan.path,model.VMIN+ii*(model.VMAX-model.VMIN)/(model.VDIVS-1))
            for jj in np.arange(0,model.CADIVS,2000):
                log.info("    Ca={} A,B={},{}",
                         model.CAMIN+jj*(model.CAMAX-model.CAMIN)/(model.CADIVS-1),
                         xGate.tableA[ii][jj], xGate.tableB[ii][jj])
    chan.tick=-1
    return chan
#Channels (model.py) includes channel function name in the dictionary


TypicalOneD = NamedList('TypicalOneD',
                             '''channel X Y Z=[] calciumPermeable=False calciumDependent=False''')
TwoD             = NamedList('TwoD',
                             '''channel X        calciumPermeable=False calciumDependent=False''')

_FUNCTIONS = {
    TypicalOneD: chan_proto,
    TwoD: BKchan_proto,
}

def make_channel(model, chanpath, params):
    func = _FUNCTIONS[params.__class__]
    return func(model, chanpath, params)

def chanlib(model):
    if not moose.exists('/library'):
        moose.Neutral('/library')
    #Adding all the channels to the library.
    chan = [make_channel(model, '/library/'+key, value) for key, value in model.Channels.items()]
    if model.ghkYN:
        ghk = moose.GHK('/library/ghk')
        ghk.T = model.Temp
        ghk.Cout = model.ConcOut
        ghk.valency=2
