from pylab import *
import pickle
from os.path import exists
from matplotlib.collections import PatchCollection
import scipy.stats
import mytools

def boxoff(ax,horoff='top'):
    ax.spines[horoff].set_visible(False)
    ax.spines['right'].set_visible(False)
    #if horoff == 'top':
    ax.get_xaxis().tick_bottom()
    #else:
    #  ax.get_xaxis().tick_top()
    ax.get_yaxis().tick_left()
    


f,axs = subplots(1,17)
axarr = axs.reshape(prod(axs.shape),).tolist()
for iax in range(0,len(axarr)):
  #for tick in axarr[iax].xaxis.get_major_ticks() + axarr[iax].yaxis.get_major_ticks():
  #  tick.label.set_fontsize(4)
  axarr[iax].tick_params(axis='both', which='major', labelsize=4.7)

#axarr[0].set_position([0.06,0.87,0.15,0.08])
#axarr[1].set_position([0.06,0.75,0.15,0.08])
axarr[0].set_position([0.01,0.8,0.18,0.19])
axarr[1].set_position([0.06,0.67,0.13,0.09])
axarr[2].set_position([0.06,0.48,0.13,0.09])
for iax in range(0,4):
  axarr[3+iax].set_position([0.27+0.195*iax,0.79,0.12,0.15])

for iax in range(0,4):
  axarr[7+iax].set_position([0.27+0.195*iax,0.52,0.12,0.15])

for iax in range(0,2):
  axarr[11+iax].set_position([0.07+0.195*iax,0.2,0.11,0.17])
  axarr[13+iax].set_position([0.07+0.195*2+0.1*iax,0.2,0.03,0.17])
  axarr[15+iax].set_position([0.67+0.195*iax,0.2,0.11,0.17])

axnew = []
for iax in range(0,2):
  axnew.append(f.add_axes([0.07+0.195*iax,0.37,0.11,0.04]))
  for sp in ['left','right','top','bottom']:
    axnew[iax].spines[sp].set_visible(False)
  axnew[iax].patch.set_alpha(0.0)
  
boxoff(axarr[1])
boxoff(axarr[2])

#Plot morphology in axarr[0]
if True:
  unpicklefile = open('morph_Nsyn50.sav', 'rb')
  unpickledlist = pickle.load(unpicklefile,encoding='bytes')
  unpicklefile.close()
  segdata = unpickledlist[:]

  for ipos in range(0,len(segdata)):
    coord1 = segdata[ipos][0]
    coord2 = segdata[ipos][1]
    mystyle = segdata[ipos][2]
    diam = segdata[ipos][3]
    mycol = '#008800'
    if mystyle in ['bx','kx','cx']:
      axarr[0].plot(coord1, coord2,mystyle,color='#00FFFF', linewidth=0.4, ms=1.4, mew=0.7)
    else:
      axarr[0].plot(coord1, coord2,mystyle,color=mycol, linewidth=diam)

  axarr[0].axis('equal')
  axarr[0].plot([100,100,200],[250,150,150],'k-',lw=0.8)
  axarr[0].text(150,140,'100 $\mu$m',fontsize=5,va='top',ha='center')
  axarr[0].text(95,200,'100 $\mu$m',fontsize=5,va='center',ha='right')

  
axarr[0].spines['top'].set_visible(False)
axarr[0].spines['right'].set_visible(False)
axarr[0].spines['bottom'].set_visible(False)
axarr[0].spines['left'].set_visible(False)
axarr[0].get_xaxis().tick_bottom()
axarr[0].get_yaxis().tick_left()
axarr[0].set_xticks([])
axarr[0].set_yticks([])

    
#1: Effect of RGS concentration on GIRK current deactivation
DATA_RGS7KO = array([[-5.916981132075471, 0.9714640198511166], [-5.615094339622641, 1.063275434243176], [-5.132075471698112, 1.0607940446650124], [-4.920754716981131, 0.4280397022332506], [-4.618867924528301, 0.04466501240694809], [-4.241509433962264, 0.04466501240694809], [-3.8641509433962256, 0.06575682382134018], [-3.50188679245283, 0.09553349875930528], [-3.094339622641509, 0.12406947890818865], [-2.6415094339622636, 0.13771712158808946], [-2.0679245283018863, 0.15632754342431765], [-1.4792452830188667, 0.16749379652605467], [-0.6943396226415093, 0.1736972704714641], [0.12075471698113205, 0.18114143920595538], [0.7547169811320771, 0.18114143920595538], [1.343396226415095, 0.18238213399503722], [1.8867924528301891, 0.1786600496277917], [2.7018867924528305, 0.17990074441687343], [3.3962264150943398, 0.18238213399503722], [4.060377358490566, 0.18362282878411917], [4.467924528301888, 0.152605459057072], [4.981132075471699, 0.15136476426799006], [5.735849056603774, 0.15880893300248144], [6.415094339622643, 0.1786600496277917], [7.1849056603773604, 0.21464019851116634], [7.924528301886793, 0.24441687344913166], [8.588679245283018, 0.2965260545905708], [9.313207547169812, 0.3771712158808934], [9.99245283018868, 0.47518610421836227], [10.656603773584907, 0.5769230769230769], [11.154716981132077, 0.6550868486352357], [11.803773584905661, 0.7444168734491314], [12.347169811320756, 0.8027295285359801], [12.99622641509434, 0.8610421836228288], [13.66037735849057, 0.8858560794044665], [14.354716981132077, 0.9119106699751861], [14.928301886792456, 0.9652605459057072], [15.335849056603772, 0.9937965260545906], [15.788679245283019, 0.966501240694789], [16.196226415094344, 1.001240694789082], [16.588679245283018, 0.9702233250620347], [17.10188679245283, 0.9937965260545906], [17.6, 0.9975186104218362], [17.87169811320755, 1.0285359801488834], [18.4, 1.0459057071960298], [18.867924528301884, 1.0198511166253101]])

DATA_WT = array([[-5.932075471698113, 0.9975186104218362],[-5.766037735849056, 0.9776674937965261],[-5.675471698113207, 1.0483870967741935],[-5.449056603773585, 1.0459057071960298],[-5.313207547169811, 1.0570719602977667],[-5.132075471698112, 1.053349875930521],[-4.950943396226414, 0.7779156327543424],[-4.860377358490566, 0.3498759305210918],[-4.739622641509433, 0.12779156327543428],[-4.437735849056603, 0.04342431761786614],[-4.10566037735849, 0.053349875930521096],[-3.6830188679245275, 0.07444168734491319],[-3.2905660377358483, 0.08312655086848642],[-2.867924528301886, 0.10545905707196046],[-2.4603773584905655, 0.12158808933002496],[-2.0377358490566033, 0.1253101736972706],[-1.6603773584905657, 0.10794044665012414],[-1.222641509433961, 0.14516129032258074],[-0.7849056603773583, 0.14019851116625326],[-0.2566037735849047, 0.1550868486352358],[0.271698113207548, 0.1550868486352358],[0.8603773584905667, 0.16873449131513663],[1.373584905660378, 0.1712158808933003],[1.7811320754716977, 0.15632754342431765],[2.1132075471698126, 0.152605459057072],[2.566037735849056, 0.16749379652605467],[3.0037735849056606, 0.1712158808933003],[3.3962264150943398, 0.17493796526054595],[3.7735849056603765, 0.17741935483870974],[4.090566037735849, 0.19106699751861045],[4.40754716981132, 0.1712158808933003],[4.679245283018869, 0.16377171215880904],[5.011320754716982, 0.1898263027295286],[5.2377358490566035, 0.23449131513647647],[5.509433962264152, 0.3126550868486353],[5.600000000000001, 0.38957816377171217],[5.811320754716981, 0.47146401985111674],[5.962264150943398, 0.5657568238213401],[6.09811320754717, 0.6699751861042184],[6.294339622641511, 0.7419354838709677],[6.475471698113209, 0.8138957816377171],[6.686792452830188, 0.8660049627791564],[6.928301886792454, 0.9181141439205955],[7.230188679245284, 0.9615384615384616],[7.577358490566038, 0.9578163771712159],[7.849056603773585, 0.9851116625310173],[8.19622641509434, 0.9751861042183623],[8.452830188679247, 0.9975186104218362],[8.89056603773585, 1.0347394540942927],[9.116981132075471, 1.0086848635235732],[9.449056603773586, 1.0124069478908189],[9.766037735849059, 0.9937965260545906],[10.067924528301885, 1.0173697270471465],[10.4, 1.0198511166253101],[10.747169811320754, 1.011166253101737],[11.154716981132077, 1.0347394540942927],[11.592452830188678, 1.032258064516129],[12.030188679245283, 1.0186104218362282],[12.37735849056604, 1.011166253101737],[12.58867924528302, 1.0347394540942927]])

for iax in [1,2]:
  polygon = Polygon(array([[0,0,10,10],[0,1,1,0]]).T)
  p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
  p.set_facecolor('#BBBBBB')
  p.set_edgecolor(None)
  axarr[iax].add_collection(p)
axarr[1].plot(DATA_WT[:,0]-min(DATA_WT[:,0]),1-DATA_WT[:,1],'k-',label='WT')
axarr[1].plot(DATA_RGS7KO[:,0]-min(DATA_WT[:,0]),1-DATA_RGS7KO[:,1],'r-',label='RGS7-KO')


coeffs = [1.0,0.2,0.1,0.05]
cols = ['#000000','#FF66FF','#00BB00','#FFFF00']
labels = ['Wild-type','20% [RGS]','10% [RGS]', '5% [RGS]']
for icoeff in [0,1,2,3]:
  unpicklefile = open('combe_p0_IC_RGSx'+str(coeffs[icoeff])+'_0.033_tstop30000_T1000_dur10000.0_somIsc1.0_max100.0.sav','rb')
  unpickledlist = pickle.load(unpicklefile,encoding='bytes')
  unpicklefile.close()
  minVClamp = unpickledlist[2][-1]
  maxVClamp = unpickledlist[2][8999]
  
  #axarr[2].plot((array(list(range(0,len(unpickledlist[2]))))-5000)/1000,1-(unpickledlist[2] - minVClamp)/(maxVClamp-minVClamp),'k-',color=cols[icoeff])
  axarr[2].plot((array(list(range(0,len(unpickledlist[2])-800))))/1000,(unpickledlist[2][800:] - minVClamp)/(maxVClamp-minVClamp),'k-',color=cols[icoeff],zorder=1+(icoeff==2),label=labels[icoeff])
for iax in [1,2]:
  axarr[iax].set_ylim([-0.01,1.03])
  axarr[iax].set_xlim([-1,30])
  axarr[iax].set_xticks([0,5,10,15,20,25])

axarr[2].set_ylabel('                                                 GIRK current (normalized)',fontsize=4.7)

cols_cases = ['#00BB00','#6666FF','#00BBBB','#BBBB00']
labels_cases = ['Wild-type', '10% [RGS]', '5-fold ligand sensitivity', '10% [RGS] & 5-fold ligand sensitivity']

#iax=1: Effect of RGS concentration on somatic resting membrane potential
#iax=2: Effect of increased GABABR agonist concentration sensitivity on somatic resting membrane potential
#iax=3: Effect of decreased RGS concentration and increased GABABR agonist concentration sensitivity on somatic resting membrane potential
polygon = Polygon(array([[2,200,200,2],[-68.5-(-64.7)+x*0.7 for x in [-1,-1,1,1]]]).T)
p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
p.set_facecolor('#BBBBBB')
p.set_edgecolor(None)
axarr[6].add_collection(p)
concs = [2e-6, 4e-6, 6e-6, 1e-5, 2e-5, 4e-5, 6e-5, 1e-4, 2e-4]
Vrests_all = []
for iax in range(0,4):
  #Vrests = array([[-70.393966,-70.400386],[-70.853872,-70.877969],[-71.280452,-71.330913],[-71.673734,-71.756596],[-72.033041,-72.151405],[-73.404088,-73.700399],[-74.853229,-75.352346],[-75.573096,-76.155422],[-76.000923,-76.617524],[-76.29005,-76.921446],[-76.988728,-77.605759]])
  Vrests = []
  for iconc in range(0,len(concs)):
    unpicklefile = open('combe_restmemb_p0_IC_RGSx'+('0.1' if iax%2 == 1 else '1.0')+'_0.033_tstop3000_T0_somIsc1.0_max'+str(round(1e6*(5 if iax > 1 else 1)*concs[iconc])/1e6)+'.sav','rb')
    unpickledlist = pickle.load(unpicklefile,encoding='bytes')
    unpicklefile.close()
    Vrest = unpickledlist[0][0]
    Vrests.append(Vrest)
    if iax == 0:
      for iax1 in range(0,3):
        axarr[3+iax1].bar(log(concs[iconc])/log(10)-0.02,Vrest,facecolor='#000000',width=0.038)
    else:
      axarr[3+iax-1].bar(log(concs[iconc])/log(10)+0.02,Vrest,facecolor=cols_cases[iax-1],width=0.038)
  Vrests_all.append(Vrests[:])
  if iax > 0:
    axarr[6].semilogx([1e6*x for x in concs], array(Vrests_all[iax])-array(Vrests_all[0]),color=cols_cases[iax-1])

#for iax in range(2,5):
for iax in [2]:
  axarr[iax].set_xlabel('time (s)',fontsize=4.7)
for iax in [3,4,5]:
  axarr[iax].set_ylim([-80,-65])
for iax in [3,4,5,7,8,9,15,16]:
  axarr[iax].set_xticks(log(array(concs))/log(10))
  axarr[iax].set_xticklabels([str(x*1e6) for x in concs],rotation=60,ha='right')
for iax in [3,4,5,6,7,8,9,10,15,16]:
  axarr[iax].set_xlabel('basal GABABR activity (A.U.)',fontsize=4.7)
for iax in [3,4,5]:
  axarr[iax].set_ylabel('Resting $V_m$ (mV)  ',fontsize=4.7)
for iax in [6]:
  axarr[iax].set_ylabel('Difference in\nresting $V_m$ (mV)  ',fontsize=4.7)
for iax in [7,8,9]:
  axarr[iax].set_ylabel('Thresh. curr. (nA)  ',fontsize=4.7)
for iax in [10]:
  axarr[iax].set_ylabel('Change in Thresh.\ncurr. (fold change)  ',fontsize=4.7)
  
  
f.savefig("fig5.pdf")


DATA_THRS = [[104.64240903387704,113.11166875784191-104.64240903387704],
             [146.23588456712673,163.5508155583438-146.23588456712673]]
polygon = Polygon(array([[2,200,200,2],[DATA_THRS[1][0]/DATA_THRS[0][0]+x*DATA_THRS[1][1]/DATA_THRS[1][0] for x in [-1,-1,1,1]]]).T)
p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
p.set_facecolor('#BBBBBB')
p.set_edgecolor(None)
axarr[10].add_collection(p)


#3: Effect of decreased RGS concentration and increased GABABR agonist concentration sensitivity on threshold current
for iax in range(0,3):
  mycol = cols_cases[iax]
  thrs = []
  thrs2 = []
  for iconc in range(0,len(concs)):
    filename = 'combe_thr_p0_IC_RGSx1.0_0.033_tstop8000_T3000_DCdur5.0_somIsc1.0_max'+str(concs[iconc])+'.sav'
    if exists(filename):
      print('Loading '+filename)
      unpicklefile = open(filename,'rb');unpickledlist = pickle.load(unpicklefile,encoding='bytes');unpicklefile.close()
      #axarr[7+iax].bar(2.5*iconc+0,unpickledlist[0],facecolor='#000000')
      axarr[7+iax].bar(log(concs[iconc])/log(10)-0.02,unpickledlist[0],facecolor='#000000',width=0.038)
      #axarr[2].bar(log(concs[iconc])/log(10)-0.02,Vrests[iconc,0],facecolor='#000000',width=0.038)
      #axarr[2].bar(log(concs[iconc])/log(10)+0.02,Vrests[iconc,1],facecolor=cols_cases[0],width=0.038)
      thrs.append(unpickledlist[0])
    else:
      print(filename+' does not exist')

    filename = 'combe_thr_p0_IC_RGSx'+str(1.0 if iax==1 else 0.1)+'_0.033_tstop8000_T3000_DCdur5.0_somIsc1.0_max'+str(round(1e6*(5 if iax > 0 else 1)*concs[iconc])/1e6)+'.sav'
    if exists(filename):
      print('Loading '+filename)
      unpicklefile = open(filename,'rb');unpickledlist = pickle.load(unpicklefile,encoding='bytes');unpicklefile.close()
      #axarr[7+iax].bar(2.5*iconc+1,unpickledlist[0],facecolor=mycol)
      axarr[7+iax].bar(log(concs[iconc])/log(10)+0.02,unpickledlist[0],facecolor=mycol,width=0.038)
      thrs2.append(unpickledlist[0])
    else:
      print(filename+' does not exist')
  axarr[10].semilogx([1e6*x for x in concs],[thrs2[i]/thrs[i] for i in range(0,len(concs))],color=cols_cases[iax])
  print('iax = '+str(iax)+'; relative thrs = '+str([thrs2[i]/thrs[i] for i in range(0,len(concs))]))

for iax in [6,10]:
  axarr[iax].set_xlim([2,200])

#4: Effect of decreased RGS concentration and increased GABABR agonist concentration sensitivity on Ca2+ inputs duriong plasticity-inducing stimuli
Nsyn = 50

iNMDAs_all = []
int_iNMDAs_all = []

for iax in range(0,2):
  iNMDAs = []
  int_iNMDAs = []
  mycol = cols_cases[2] if iax > 0 else '#000000'
  #mycol = cols_cases[iax-1] if iax > 0 else '#000000'
  for iconc in range(0,len(concs)):
    iNMDAs_thisconc = []
    int_iNMDAs_thisconc = []
    for istim in range(0,2):
      iNMDA_samps = []
      int_iNMDA_vals_samps = []
      for myseed in range(1,6):
        filename = 'combe_rhythmic_syn_inputs_inmdarec_carec_p0_IC_RGSx'+str(1.0 if iax==0 else 0.1)+'_0.033_tstop17000_T10000_'+('Ninputs12_Tinputs500.0' if istim==0 else 'Ninputs600_Tinputs10.0')+'_Nsyn'+str(Nsyn)+'_somIsc1.0_max'+str(round(1e6*(5 if iax > 0 else 1)*concs[iconc])/1e6)+'_seed'+str(myseed)+'.sav'
        #filename = 'combe_rhythmic_syn_inputs_inmdarec_p0_IC_RGSx'+str(1.0 if iax%2==0 else 0.1)+'_0.033_tstop17000_T10000_'+('Ninputs12_Tinputs500.0' if istim==0 else 'Ninputs600_Tinputs10.0')+'_Nsyn'+str(Nsyn)+'_somIsc1.0_max'+str(round(1e6*(5 if iax > 1 else 1)*concs[iconc])/1e6)+'_seed1.sav'
        if not exists(filename):
          print(filename+' does not exist')
          continue
        print('Loading '+filename)
        unpicklefile = open(filename,'rb');unpickledlist = pickle.load(unpicklefile,encoding='bytes');unpicklefile.close()
        spikes_all,times,Vsoma,iNMDA,casyn = unpickledlist[0], unpickledlist[1], unpickledlist[2], unpickledlist[3], unpickledlist[4]
        iNMDA = iNMDA*0.1*1000 #Consider only the Ca currents (estimated to be 10% of all NMDA currents; PMID: 7666365. Converted from nA to pA)
        #int_iNMDA_vals = [max(abs(mean(iNMDA,axis=1))[i*500:i*500+10]) for i in range(0,12)] if istim == 0 else [max(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(0,600)]
        #int_iNMDA_vals = [max(abs(mean(iNMDA,axis=1))[i*500:i*500+10]) for i in range(0,12)] if istim == 0 else [max(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(0,12)]
        #int_iNMDA_vals = [mean(abs(mean(iNMDA,axis=1))[i*500:i*500+500]) for i in range(0,12)] if istim == 0 else [mean(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(588,600)]
        #int_iNMDA_vals = [mean(abs(mean(iNMDA,axis=1))[i*500:i*500+500]) for i in range(0,12)] if istim == 0 else [mean(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(0,600)]
        #int_iNMDA_vals = [mean(abs(mean(iNMDA,axis=1))[i*500:i*500+500]) for i in range(1,12)] if istim == 0 else [mean(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(50,600)] #Ignore the first 500 ms (the first stim in LFS, the first 50 stims in HFS)
        #int_iNMDA_vals = [mean(abs(mean(iNMDA,axis=1))[i*500:i*500+500]) for i in range(2,12)] if istim == 0 else [mean(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(100,600)] #Ignore the first 1000 ms (the first 2 stims in LFS, the first 100 stims in HFS)
        int_iNMDA_vals = [mean(abs(mean(iNMDA,axis=1))[i*500:i*500+500]) for i in range(2,12)] if istim == 0 else [mean(abs(mean(iNMDA,axis=1))[i*10:i*10+10]) for i in range(590,600)] #Consider the last 10 in both cases
        iNMDA_samps.append(iNMDA[:])
        int_iNMDA_vals_samps.append(int_iNMDA_vals[:])
        print('Av. iNMDA = '+str(mean(iNMDA))+', av. iNMDA integral = '+str(mean(int_iNMDA_vals))+' for iax = '+str(iax)+", istim="+str(istim)+' ;myseed = '+str(myseed))

      if concs[iconc] == 1e-5:
        if iax == 0:
          axnew[istim].plot(times-10000,Vsoma,'k-',lw=0.3)
        iconc_saved = iconc
        axarr[11+istim].plot(range(0,len(iNMDA)),mean(iNMDA,axis=1),'k-',lw=0.1,color=mycol)

        #axarr[11+istim].plot(range(0,len(casyn)),mean(casyn,axis=1),'k-',lw=0.1,color=mycol)
        axarr[13+istim].bar(2.5*iconc+iax,mean(mean(array(int_iNMDA_vals_samps),axis=1),axis=0),facecolor=mycol)
        #axarr[13+istim].bar(2.5*iconc+iax,mean(int_casyns_vals),facecolor=mycol)
        axarr[13+istim].plot([2.5*iconc+iax]*2,[mean(mean(array(int_iNMDA_vals_samps),axis=1),axis=0)-std(mean(array(int_iNMDA_vals_samps),axis=1),axis=0),mean(mean(array(int_iNMDA_vals_samps),axis=1),axis=0)+std(mean(array(int_iNMDA_vals_samps),axis=1),axis=0)],
                             'k-',lw=0.5)
        #axarr[13+istim].plot([2.5*iconc+iax]*2,[mean(int_casyns_vals)-std(int_casyns_vals),mean(int_casyns_vals)+std(int_casyns_vals)],'k-',lw=0.5)
      int_iNMDAs_thisconc.append(int_iNMDA_vals_samps[:])
      axarr[15+istim].bar(log(concs[iconc])/log(10)-0.02+0.04*iax,mean(int_iNMDAs_thisconc[-1]),facecolor=mycol,width=0.038)
      #axarr[15+istim].bar(log(concs[iconc])/log(10)-0.02+0.04*iax,mean(int_casyns_thisconc[-1]),facecolor=mycol,width=0.038)
      #qwe  

      iNMDAs_thisconc.append(iNMDA[:])
    iNMDAs.append(iNMDAs_thisconc[:])
    int_iNMDAs.append(int_iNMDAs_thisconc[:])
  iNMDAs_all.append(iNMDAs[:])
  int_iNMDAs_all.append(int_iNMDAs[:])

  
#for istim in [0,1]:
#  pval = scipy.stats.ranksums(mean(iNMDAs_all[0][iconc_saved][istim],axis=0),mean(iNMDAs_all[1][iconc_saved][istim],axis=0))[1]
#  print('istim = '+str(istim)+'; Ns = '+str(iNMDAs_all[0][iconc_saved][istim].shape[1])+','+str(iNMDAs_all[1][iconc_saved][istim].shape[1])+': pval = '+str(pval)+'; means = '+str(mean(iNMDAs_all[0][iconc_saved][istim]))+','+str(mean(iNMDAs_all[1][iconc_saved][istim]))+' (iNMDA)')
#  pval = scipy.stats.ranksums(mean(casyns_all[0][iconc_saved][istim],axis=0),mean(casyns_all[1][iconc_saved][istim],axis=0))[1]
#  print('istim = '+str(istim)+'; Ns = '+str(casyns_all[0][iconc_saved][istim].shape[1])+','+str(casyns_all[1][iconc_saved][istim].shape[1])+': pval = '+str(pval)+'; means = '+str(mean(casyns_all[0][iconc_saved][istim]))+','+str(mean(casyns_all[1][iconc_saved][istim]))+' (cai)')
#  #if pval < 0.01:
#  #  axarr[12+istim].plot([2.5*iconc_saved+x for x in [0,0,1,1]],[(6e-4+6e-5*x) if istim == 0 else (0.0018+0.00018*x) for x in [0,1,1,0]], 'k-',lw=0.5)
#  #  axarr[12+istim].plot(2.5*iconc_saved+0.5, 6e-4+6e-5*2 if istim == 0 else 0.0018+0.00018*2 , 'k*', lw=0.5,ms=3,mew=0.5)

for istim in [0,1]:
  data1 = mean(array(int_iNMDAs_all[0][iconc_saved][istim]),axis=1)
  data2 = mean(array(int_iNMDAs_all[1][iconc_saved][istim]),axis=1)
  pval = scipy.stats.ranksums(data1,data2)[1]
  print('istim = '+str(istim)+'; int Ns = '+str(len(int_iNMDAs_all[0][iconc_saved][istim]))+','+str(len(int_iNMDAs_all[1][iconc_saved][istim]))+': pval = '+str(pval)+'; means = '+str(mean(int_iNMDAs_all[0][iconc_saved][istim]))+','+str(mean(int_iNMDAs_all[1][iconc_saved][istim]))+' (iNMDA)')
  if pval < 0.01:
    axarr[13+istim].plot([2.5*iconc_saved+x for x in [0,0,1,1]],[(0.01+1e-3*x) if istim == 0 else (0.18+0.018*x) for x in [0,1,1,0]], 'k-',lw=0.5)
    axarr[13+istim].plot(2.5*iconc_saved+0.5, 0.01+1e-3*2 if istim == 0 else 0.18+0.018*2 , 'k*', lw=0.5,ms=3,mew=0.5)

for iax in [11,12]:
  axarr[iax].set_xlim([-50,1500])
  axarr[iax].set_xlabel('$t$ (ms)',fontsize=4.7)
for iax in [0,1]:
  axnew[iax].set_xlim([-50,1500])
  axnew[iax].set_ylim([-80,-20])
  axnew[iax].set_xticks([])
  axnew[iax].set_yticks([])
  axnew[iax].plot([-150,-150],[-75,-25],'k-',lw=0.7,clip_on=False)
for iax in [13,14]:
  axarr[iax].set_xticks([])
#for iax in [10]: #,12,14]:
#  axarr[iax].set_ylim([-0.0007,0])
#for iax in [12,14]: #,12,14]:
#  axarr[iax].set_ylim([-0.00008,0])
#for iax in [11]: #,13,15]:
#  axarr[iax].set_ylim([-0.0015,0])
#for iax in [13,15]: #,13,15]:
#  axarr[iax].set_ylim([-0.0012,0])

for iax in [11,12]:
  axarr[iax].set_ylabel('$I_{\mathrm{NMDA,Ca}^{2+}}$ (pA)',fontsize=4.7)
for iax in [13,14,15,16]:
  axarr[iax].set_ylabel('Av. $I_{\mathrm{NMDA,Ca}^{2+}}$ (pA)',fontsize=4.7)
  #axarr[iax].set_ylabel('[Ca$^{2+}$ (mM)',fontsize=5)
for iax in [0,1,3,4,5,6,7,8,9,10,11,12,13,14,15,16]:
  pos = axarr[iax].get_position()
  f.text(pos.x0 - 0.03 - 0.02 - 0.015*(iax in [3,4,5,7,8,9]) - 0.015*(iax in [11,12,13,16]) + 0.045*(iax==0) - 0.008*(iax==13), pos.y1 - 0.02 + 0.005*(iax in [7,8,9,10]) - 0.01*(iax==0), chr(ord('A')+iax-(iax>1)), fontsize=11)

axarr[14].set_yticks([0,0.05,0.1])
axarr[15].set_yticks([0,0.003,0.006])
axarr[16].set_yticks([0,0.05,0.1])
#axarr[0].legend(fontsize=5)
#axarr[1].legend(fontsize=5)

legs = []
#legs.append(mytools.mylegend(f,[0.04,0.95,0.15,0.035],['k-','r-'],['WT','RGS7-KO'],nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=['#000000','#FF0000'],dashes=[],linewidths=[],myfontsize=4.7))
#legs.append(mytools.mylegend(f,[0.04,0.73,0.15,0.05],['k-','k-','k-','k-'],labels,nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=cols,dashes=[],linewidths=[],myfontsize=4.7))
legs.append(mytools.mylegend(f,[0.04,0.763,0.15,0.035],['k-','r-'],['WT','RGS7-KO'],nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=['#000000','#FF0000'],dashes=[],linewidths=[],myfontsize=4.7))
legs.append(mytools.mylegend(f,[0.04,0.573,0.15,0.05],['k-','k-','k-','k-'],labels,nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=cols,dashes=[],linewidths=[],myfontsize=4.7))
legs.append(mytools.mylegend(f,[0.3,0.945,0.65,0.05],['k-','k-','k-','k-'],labels_cases,nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=['#000000']+cols_cases[0:3],dashes=[],linewidths=[],myfontsize=4.7))
for leg in legs:
  for q in ['left','right','top','bottom']:
    leg.spines[q].set_visible(False)

refDCdur = 100.0
if len(sys.argv) > 1:
  refDCdur = float(sys.argv[1])
  
f.savefig("fig5.pdf")