0

tl;dr trying to create plots using for loop in Python 3 but it is resulting in weird inconsistent plots.

I am trying to create plots like Figure 2 of Fox et al. 2015 (https://arxiv.org/pdf/1412.1480.pdf). I have the different ion lines saved in a dictionary with their parameters saved as lists. For example,

DictIon = { 
'''The way this dictionary is defined is as follows: ion name: 
[wavelength, f-value, continuum lower 1, continuum higher 1, continuum lower 2, 
continuum higher 2, subplot x, subplot y, ion name x, ion name y]'''
'C II 1334': [1334.5323,.1278,-400,-200,300,500,0,0,-380,.2],
'Si II 1260': [1260.4221,1.007,-500,-370,300,500,2,0,-380,.15],
'Si II 1193': [1193.2897,0.4991,-500,-200,200,500,3,0,-380,.2]
}

In order to generate the plot, I have written the following code where at first I calculate the continuum of the absorption lines and then use data from a Voigt Profile algorithm that gives me the zval and bval that I use in the Voigt profile section of the code, and finally I combine these two and plot the values in appropriate subplots. The code is:

f, axes = plt.subplots(6, 2, sharex='col', sharey='row',figsize = (15,15))

for ion in DictIon:
    w0 = DictIon[ion][0] #Rest wavelength
    fv = DictIon[ion][1] #Oscillator strengh/f-value

    velocity = (wavelength-Lambda)/Lambda*c

    #Fit the continuum
    low1  = DictIon[ion][2]  #lower bound for continuum fit on the left side
    high1 = DictIon[ion][3]  #upper bound for continuum fit on the left side
    low2  = DictIon[ion][4]  #lower bound for continuum fit on the right side
    high2 = DictIon[ion][5]  #upper bound for continuum fit on the right side
    x1 = velocity[(velocity>=low1) & (velocity<=high1)]
    x2 = velocity[(velocity>=low2) & (velocity<=high2)]
    X = np.append(x1,x2)
    y1 = flux[(velocity>=low1) & (velocity<=high1)]
    y2 = flux[(velocity>=low2) & (velocity<=high2)]
    Y = np.append(y1,y2)
    Z = np.polyfit(X,Y,1)

    #Generate data to plot continuum
    xp = np.linspace(-500,501,len(flux[(velocity>=-500) & (velocity<=500)]))
    p = np.poly1d(Z)

    #Normalize flux
    norm_flux = flux[(velocity>=-500) & (velocity<=500)]/p(xp)

    #Create a line at y=1
    tmp1 = np.linspace(-500,500,10)
    tmp2 = np.full((1,10),1)[0]


    '''Generate Voigt Profile Fits'''
    #Initialize arrays
    vmod = (np.arange(npix+1)-(npix/pixsize))*pixsize #-npix to npix in steps of pix
    fmodraw = np.ndarray((npix+1)); fmodraw.fill(1.0)
    ncom = len(zval)+1

    fitn = 10**(logn); fitne = 10**(logn+elogn)-10**logn
    totcol=np.log10(np.sum(fitn))
    etotcol=np.sqrt(np.sum(elogn**2)) #strictly only true if independent

    #Set up arrays
    sigma=bval/np.sqrt(2.0); tau0=np.ndarray((ncom)); tauv=np.ndarray((ncom,npix))
    find=np.ndarray((ncom, npix)); sfind=np.ndarray((ncom, npix)) #smoothed

    #go from z to velocity
    v0=c*(zval-zmod) / (1.0+zmod); ev0=c*zvale  / (1.0+zmod) 
    bv=bval; ebv=bvale

    #generate models for each comp, where tau is Gaussian with velocity
    for k in range(0, ncom-1):
        tau0[k]=(1.497e-2*(10**logn[k])*(w0/1.0e8)*fv) / (bval[k]*1.0e5)
        for j in range(0, npix-1):
            tauv[k][j]=tau0[k]*np.exp((-1.0*(vmod[j]-v0[k])**2) / (2.0*sigma[k]**2))
            find[k][j]=np.exp(-1.0*tauv[k][j])

    #Transpose
    tauv = tauv.T
    find = find.T

    #Sum over components (pixel by pixel)
    tottauv=np.ndarray((npix+1))
    for j in range(0, npix-1):
        tottauv[j]=tauv[j,:].sum()

    fmodraw=np.exp(-1.0*tottauv)

    #create Gaussian kernel (smoothing function or LSF) 
    #created on 1 km/s grid with default FWHM=20.0 km/s (UVES), integral=1 
    fwhmins=20.0

    sigins=fwhmins/(1.414*1.665); nker=150 #NEED TO FIND NKER
    vt=np.arange(nker)-nker/2 #-75 to +75 in 1 km/s steps
    smfn=(1.0/(sigins*np.sqrt(2.0*np.pi)))*np.exp((-1.0*(vt**2))/(2.0*sigins**2))

    #convolve total and individual comps with LSF
    fmod = np.convolve(fmodraw, smfn, mode='same')

    axes[DictIon[ion][6]][DictIon[ion][7]].axis([-400,400,-.12,1.4])
    axes[DictIon[ion][6]][DictIon[ion][7]].xaxis.set_major_locator(xmajorLocator)
    axes[DictIon[ion][6]][DictIon[ion][7]].xaxis.set_major_formatter(xmajorFormatter)
    axes[DictIon[ion][6]][DictIon[ion][7]].xaxis.set_major_locator(xminorLocator)
        axes[DictIon[ion][6]][DictIon[ion][7]].yaxis.set_major_locator(ymajorLocator)
        axes[DictIon[ion][6]][DictIon[ion][7]].yaxis.set_major_locator(yminorLocator)

    axes[DictIon[ion][6]][DictIon[ion][7]].plot(vmod,fmod,'r', linewidth = 1.5)
    axes[DictIon[ion][6]][DictIon[ion][7]].plot(tmp1,tmp2,'k--')
    axes[DictIon[ion][6]][DictIon[ion][7]].step(velocity[(velocity>=-500) & (velocity<=500)],norm_flux,'k')

Instead of producing plots as Figure 2 of Fox et al. 2015, I am getting situations like this where the code will produce different results at different times when I run it: Fig1 Fig2][2

Bottom line, I have tried to debug this for 3 days now and I am at a loss. I suspect that it may have something to do with how pyplot plots work in for loop and the fact that I am using a dictionary to loop through. Any advice or suggestions would be greatly appreciated. I am using Python 3.

Edit: Data available here: zval,bval values: https://drive.google.com/file/d/0BxZ6b2fEZcGBX2ZTUEdDVHVWS0U/

velocity, flux values:

Si III 1206: https://drive.google.com/file/d/0BxZ6b2fEZcGBQXpmZ01kMDNHdk0/

Si IV 1393: https://drive.google.com/file/d/0BxZ6b2fEZcGBamkxVVA2dUY0Qjg/

3
  • Could you provide a set of sample data? (like velocity, zval, etc.) Commented May 3, 2017 at 3:26
  • I have added links to data in the edit Commented May 3, 2017 at 3:59
  • This could do with being re-asked with a minimum example, i.e., I presume the problem is about using matplotlib in a for loop, not about your complicated data-set etc? Commented May 3, 2017 at 4:59

1 Answer 1

1

Here's a MWE of the operation you want to perform:

import matplotlib.pyplot as plt
import numpy as np

f, axes = plt.subplots(3, 2, sharex='col', sharey='row',figsize=(15,15))
plt.subplots_adjust(hspace=0.)

data_dict = {'C II 1334': [1., 2.],
             'Si II 1260': [2., 3.],
             'Si II 1193': [3., 4.]}

for i, (key, value) in enumerate(data_dict.items()):

    print i, key, value    

    x = np.linspace(0, 100, 10)
    y0 = value[0] * x
    y1 = value[1] * x

    axes[i, 0].plot(x, y0, label=key)
    axes[i, 1].plot(x, y1, label=key)
    axes[i, 0].legend(loc="upper right")
    axes[i, 1].legend(loc="upper right")

plt.legend()
plt.show()

Results in

enter image description here

I don't see any strange behaviour after invoking plt in a for loop over a dict.

I suggest you separate data processing/calculations, i.e. calculations of the scientific quantities of interest, from plotting said data.

Note that the ordering of the dictionary items isn't preserved - better to use a list or an ordered dictionary in my example.

Sign up to request clarification or add additional context in comments.

1 Comment

I took your advice and turned it into a list and it seems to work now. The old code still doesn't work so I am guessing it probably has to do something with the unordered nature of dictionaries. Thanks for the tip!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.