import matplotlib
matplotlib.use('Agg')
import numpy
from pylab import *
import mytools
import pickle
import sys
import time
from os.path import exists
coeffCoeffs = [[0.25,0],[0.125,0],[0.5,0],[0.5,1.0/3],[0.5,2.0/3],[0.5,1.0],[-0.25,0],[-0.125,0],[-0.5,0]]
import mutation_stuff
MT = mutation_stuff.getMT()
defVals = mutation_stuff.getdefvals()
geneNames = mutation_stuff.getgenenames()
keyList = defVals.keys()
for idefval in range(0,len(keyList)):
if type(defVals[keyList[idefval]]) is not list:
defVals[keyList[idefval]] = [defVals[keyList[idefval]], defVals[keyList[idefval]]] #make the dictionary values [somatic, apical]
updatedVars = ['somatic','apical','basal'] # the possible classes of segments that defVals may apply to
whichDefVal = [0,1,0] # use the defVal[0] for somatic and basal segments and defVal[1] for apical segments
variants = [[0,5,1],[2,4,7],[1,2,13],[3,1,0],[5,0,0],[8,3,0],[12,2,0],[13,4,0]] #from drawallmeangains (maxCountersAll)
oscamp = 0.25
## These are the frequencies used in the article figure 2:
#oscfreqs = [0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 7.5, 10.0, 15.0]
# This is a smaller set of frequencies that captures most of the shape of the response curve:
oscfreqs = [0.5, 0.625, 0.75, 0.875, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 7.5, 10.0, 15.0]
gsyn = 1.07
gNoise = 1.07
fts_control = []
ft_df = 0.00001
if exists('spectrum_freq1.0_0.sav'): #Do this just to get the fs_control
unpicklefile = open('spectrum_freq1.0_0.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fs_control = unpickledlist[0]
print "Loading control..."
amps_control = []
amps_control_std = []
if exists('spectra_fig2_0.sav'):
unpicklefile = open('spectra_fig2_0.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fts = unpickledlist[0]
fts_std = unpickledlist[1]
amps_control = unpickledlist[2]
amps_control_std = unpickledlist[3]
else:
for iosc in range(0,len(oscfreqs)):
oscfreq = oscfreqs[iosc]
if exists('spectrum_freq'+str(oscfreq)+'_0.sav'):
unpicklefile = open('spectrum_freq'+str(oscfreq)+'_0.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fs_control = unpickledlist[0]
ft_fs = unpickledlist[1]
fts = [mean([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs_control))]
fts_std = [std([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs_control))]
iosc_fs = next(i for i,x in enumerate(fs_control) if x>oscfreq/1000)
amps_control.append(fts[iosc_fs])
amps_control_std.append(fts_std[iosc_fs])
else:
print 'spectrum_freq'+str(oscfreq)+'_0.sav not found!'
print "Loading done"
picklelist = [fts, fts_std, amps_control, amps_control_std]
file = open('spectra_fig2_0.sav','w')
pickle.dump(picklelist,file)
file.close()
styles = ['g-','g-','g-','g-','g-','g-','g-','g-','g-']
cols = ['#666666','#012345','#cc00aa','#bbaa00','#ee6600','#ff0000', '#00aaaa','#772277','#00cc00']
col_control = '#2222ff'
counter = -1
close("all")
f, axarr = plt.subplots(1, len(variants)+1)
lenvarper2 = 4
for ix in range(0,lenvarper2):
for iy in range(0,2):
axarr[ix+lenvarper2*iy].set_position([0.08+0.176*ix, 0.1+0.44*(1-iy), 0.176, 0.37])
axarr[8].set_position([0.08+0.176*4, 0.1+0.44*(1-1), 0.176, 0.37])
timesAll = []
VsomaAll = []
spikeFreqsAll = []
for igene in range(0,len(MT)):
for imut in range(0,len(MT[igene])):
nVals = len(MT[igene][imut])*[0]
thesemutvars = []
for imutvar in range(0,len(MT[igene][imut])):
thesemutvars.append(MT[igene][imut][imutvar][0])
if type(MT[igene][imut][imutvar][1]) is int or type(MT[igene][imut][imutvar][1]) is float:
MT[igene][imut][imutvar][1] = [MT[igene][imut][imutvar][1]]
nVals[imutvar] = len(MT[igene][imut][imutvar][1])
cumprodnVals = cumprod(nVals)
allmutvars = cumprodnVals[len(MT[igene][imut])-1]*[thesemutvars]
allmutvals = []
for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
allmutvals.append([0]*len(thesemutvars))
for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
for imutvar in range(0,len(MT[igene][imut])):
if imutvar==0:
allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][iallmutval%nVals[imutvar]]
else:
allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][(iallmutval/cumprodnVals[imutvar-1])%nVals[imutvar]]
for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
mutval = allmutvals[iallmutval]
thisCoeff = 1
mutText = ""
for imutvar in range(0,len(MT[igene][imut])):
if imutvar > 0 and imutvar%2==0:
mutText = mutText+"\n"
mutvars = allmutvars[iallmutval][imutvar]
mutvals = allmutvals[iallmutval][imutvar]
if type(mutvars) is str:
mutvars = [mutvars]
mutText = mutText + str(mutvars) + ": "
for kmutvar in range(0,len(mutvars)):
mutvar = mutvars[kmutvar]
if mutvar.find('offm') > -1 or mutvar.find('offh') > -1 or mutvar.find('ehcn') > -1:
newVal = [x+mutvals*thisCoeff for x in defVals[mutvar]]
if mutvals >= 0 and kmutvar==0:
mutText = mutText + "+" + str(mutvals) +" mV"
elif kmutvar==0:
mutText = mutText + str(mutvals) +" mV"
else:
newVal = [x*(mutvals**thisCoeff) for x in defVals[mutvar]]
if kmutvar==0:
mutText = mutText + "*" + str(mutvals)
if kmutvar < len(mutvars)-1:
mutText = mutText + ", "
if igene == 0 and imut == 0 and iallmutval == 0:
iters = [-1, 0, 2, 6, 8]
else:
iters = [0, 2, 6, 8]
doSkip = True
ivar = -1
for iiter in range(0,len(iters)):
iter = iters[iiter]
counter = counter+1
for ivar2 in range(0,len(variants)):
if igene == variants[ivar2][0] and imut == variants[ivar2][1] and iallmutval == variants[ivar2][2]:
ivar = ivar2
doSkip = False
break
if doSkip:
continue
if exists('spectra_fig2_'+str(counter)+'.sav'):
unpicklefile = open('spectra_fig2_'+str(counter)+'.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fts = unpickledlist[0]
fts_std = unpickledlist[1]
amps = unpickledlist[2]
amps_std = unpickledlist[3]
else:
amps = []
amps_std = []
for iosc in range(0,len(oscfreqs)):
oscfreq = oscfreqs[iosc]
iosc_fs = next(i for i,x in enumerate(fs_control) if x>oscfreq/1000)
if exists('spectrum_freq'+str(oscfreq)+'_'+str(counter)+'.sav'):
unpicklefile = open('spectrum_freq'+str(oscfreq)+'_'+str(counter)+'.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fs = unpickledlist[0]
ft_fs = unpickledlist[1]
fts = [mean([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs))]
fts_std = [std([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs))]
amps.append(fts[iosc_fs])
amps_std.append(fts_std[iosc_fs])
else:
print 'spectrum_freq'+str(oscfreq)+'_'+str(counter)+'.sav does not exist'
amps.append([])
amps_std.append([])
picklelist = [fts, fts_std, amps, amps_std]
file = open('spectra_fig2_'+str(counter)+'.sav','w')
pickle.dump(picklelist,file)
file.close()
if len(amps) > 0:
if not any([type(amps[iosc]) is list for iosc in range(0,len(amps))]):
axarr[ivar].semilogx(oscfreqs, [amps[iosc] for iosc in range(0,len(amps))], 'b-', color=cols[iter])
else:
myvec = zeros([len(amps),])
for iosc in range(0,len(amps)):
if type(amps[iosc]) is not list:
myvec[iosc] = amps[iosc]
else:
myvec[iosc] = nan
axarr[ivar].semilogx(oscfreqs, myvec, 'b.-', color=cols[iter])
print "Some frequencies not found in counter="+str(counter)+", iter="+str(iter)+"!"
else:
print "No frequencies not found in counter="+str(counter)+", iter="+str(iter)+"!"
print "igene="+str(igene)+", imut="+str(imut)+", iallmutval="+str(iallmutval)+", ivar="+str(ivar)+", counter="+str(counter)+", iter="+str(iter)
if doSkip:
continue
print mutText
axarr[ivar].semilogx(oscfreqs, amps_control, 'b-', color=col_control)
axarr[ivar].set_title(geneNames[igene])
axarr[ivar].set_xlim([0.4,15])
axarr[ivar].set_ylim([0,7.5e7])
axarr[ivar].set_yticks([0, 3e7, 6e7])
if ivar < lenvarper2:
axarr[ivar].set_xticklabels(['']*len(axarr[ivar].get_xticks()))
elif ivar == 6:
axarr[ivar].set_xlabel('Input frequency $f$ (Hz)')
if ivar % lenvarper2 > 0:
axarr[ivar].set_yticklabels(['', '', ''])
elif ivar == 0:
axarr[ivar].set_ylabel('Power of the frequency component corresponding $f$ ')
f.savefig("fig2c.eps")
iters = [0, 2, 6, 8]
combmutIDnums = [0, 1, 2, 3]
ivar = 8
amps_thismut = []
for iiter in range(0,4):
iter = iters[iiter]
if exists('spectra_fig2_comb'+str(combmutIDnums[iiter])+'.sav'):
unpicklefile = open('spectra_fig2_comb'+str(combmutIDnums[iiter])+'.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fts = unpickledlist[0]
fts_std = unpickledlist[1]
amps = unpickledlist[2]
amps_std = unpickledlist[3]
else:
amps = []
amps_std = []
for iosc in range(0,len(oscfreqs)):
oscfreq = oscfreqs[iosc]
iosc_fs = next(i for i,x in enumerate(fs_control) if x>oscfreq/1000)
if exists('spectrum_freq'+str(oscfreq)+'_comb'+str(combmutIDnums[iiter])+'.sav'):
unpicklefile = open('spectrum_freq'+str(oscfreq)+'_comb'+str(combmutIDnums[iiter])+'.sav','r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
fs = unpickledlist[0]
ft_fs = unpickledlist[1]
fts = [mean([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs))]
fts_std = [std([abs(ft_fs[isamp][ifreq])**2 for isamp in range(0,len(ft_fs))]) for ifreq in range(0,len(fs))]
amps.append(fts[iosc_fs])
amps_std.append(fts_std[iosc_fs])
else:
print 'spectrum_freq'+str(oscfreq)+'_comb'+str(combmutIDnums[iiter])+'.sav does not exist'
amps.append(nan)
amps_std.append(nan)
picklelist = [fts, fts_std, amps, amps_std]
file = open('spectra_fig2_comb'+str(combmutIDnums[iiter])+'.sav','w')
pickle.dump(picklelist,file)
file.close()
if len(amps) > 0:
if not any([type(amps[iosc]) is list for iosc in range(0,len(amps))]):
axarr[ivar].semilogx(oscfreqs, [amps[iosc] for iosc in range(0,len(amps))], 'b-', color=cols[iter])
else:
myvec = zeros([len(amps),])
for iosc in range(0,len(amps)):
if type(amps[iosc]) is not list:
myvec[iosc] = amps[iosc]
else:
myvec[iosc] = nan
axarr[ivar].semilogx(oscfreqs, myvec, 'b.-', color=cols[iter])
print "Some frequencies not found in counter="+str(counter)+", iter="+str(iter)+"!"
else:
print "No frequencies not found in counter="+str(counter)+", iter="+str(iter)+"!"
amps_thismut.append(amps[:])
axarr[ivar].semilogx(oscfreqs, amps_control, 'b-', color=col_control)
axarr[ivar].set_ylim([0,8e7])
axarr[ivar].set_title("Combination")
axarr[ivar].set_xlim([0.4,15])
axarr[ivar].set_ylim([0,7.5e7])
axarr[ivar].set_yticks([0, 3e7, 6e7])
axarr[ivar].set_yticklabels(['', '', ''])
for ivar in range(0,9):
t = axarr[ivar].yaxis.get_offset_text()
if type(t) is matplotlib.text.Text:
if t.get_text()[0:2] == '1e':
t.set_text('$\\times 10^{'+t.get_text()[2:]+'}$')
t.set_position((t.get_position()[0]-0.15,t.get_position()[1]))
f.text(0.01, 0.9, 'C', fontsize=31)
f.savefig("fig2c.eps")
