## Python Notes

The following are a few notes on the syntax, commands, etc. for python programming Assumes python version 3.x, which differs substantially from the 2.x versions. Python is not backwards-compatible, so 2.x codes don't necessarily work in 3.x.


python prog.py                   # Run python, where prog.py is a python program
python prog.py arg1 arg2         # Run python program with arguments
#  inside the code you must 'import sys', and then
#    str(sys.argv) is 'prog.py'
#    str(sys.argv) is the first argument, etc.

print("It is possible to break long lines of python code",  # breaking long lines
"by preceding the carriage return with a (, \, or \
comma")


---------------Numbers, Statistics, Vectors--------------------

# + - / * do what you'd think. ** is a power, abs is absolute value
# pow(2,7) is the same as 2**7, max and min do what you'd think
# round(2.71828,4) gives 2.7183
# 10%3         # 10mod(3) = 1
# str(10%3)    # '1' returned as a string
# float('3.44') # converts a string to a floating point number
# NOTE: ALL PYTHON inputs from the user are interpreted as strings

# Running it from the command line; most useful as a calculator
python                           # Start python (or python3 for version 3.3)
>>> import numpy as np           # import package numpy
>>> print(np.sin(np.pi/4.))      # prints 0.7071067...; need np.sin as sin is not defined
>>> a=500.3
>>> b=3
>>> a/b                          # prints out 166.7666
>>> _/5.                         # prints 33.3533   _ is the result from the last command
>>> Ctl-D                        # quit

x=[1.3,2.4,8]           # is a list, these behave like strings
x=(1.3,2.4,8)           # is a tuple, these are like lists that can't be changed
x=np.array([1.3,2.4,8]) # is what you want for a vector (1D array)
x.T                     # transpose of x, where x is a numpy array
np.zeros(10, dtype=np.int)  # initialize an integer array of 10 elements to zero
np.ceil([2.87])   # rounds up to 3
np.array([2.87,4.8], dtype=int) # rounds to [2,4]  !!Use this, and not np.int!!
np.floor([2.87,4.8]) # rounds to [2.,4.] but these are reals
np.arange(0,5,0.1) # makes [0, 0.1, 0.2 ... 4.9, 5.0] vector
np.linspace(0,100,5) # makes [0, 25, 50, 75, 100] vector
np.random.rand(5) # five random numbers between 0 and 1
np.random.randn(5) # five random numbers from std normal distribution
np.random.randint(0,99,5) # five random integers between 0 and 99
np.mean(x)       # mean of x
np,median(x)     # median of x
np.var(x)        # variance of x   DONT FORGET THE 'np.'
np.std(x)        # stddev of x
np.log(x)        # ln(x)
np.log10(x)      # log10(x)
np.exp(x)        # exp(x)
np.arctan(x)     # atan(x)
np.nanmedian(x)  # Take the median of values that are not nan.
np.percentile(a,85) # value of the 85th percentile of array a. Very useful for finding 'sigma'
# in the presence of outliers! e.g. sigma=np.percentile(a,83.3)-np.percentile,50)

import time
t0=time.time()  #Current time. Place at beginning and end to time a section of code


---------------Plotting--------------------
import matplotlib.pyplot as plt  #MATLAB-like plots
from matplotlib import pyplot as plt # Equivalent syntax for importing
import matplotlib.ticker as ticker  #for full control over tick mark placement
plt.rcParams['font.size']=14     # make the default labels bigger

fig=plt.figure(figsize=(8,15))  #fig is a figure object, change its default shape
ax1=fig.add_subplot(122) #1 row 2 columns, access plot 2  ;  ax1 is a plotting object, a container within fig
fig, [[ax1, ax2], [ax3, ax4]] = plt.subplots(2, 2, figsize=(6, 6)) # another way to set up 4 plotting objects in a 2x2
ax1.plot(data['B-V'],data['mag'], 'bx', markersize=1.3)  #plot blue x's. r-- is red dashed lines, etc., size of points is 1.3
ax1.plot(x,y,lw=2.2)  #Thicker line width
ax1.set_ylim(-2,2) #another way to set y limits
ax1.set_yticks([0,3,6,9,12,15,18])
ax1.tick_params('x',direction='in',length=5,color='b',top='on',
labelrotation=30,width=3,labelsize=10,pad=6,labelcolor='k') #tick mark size, color, angle, location
ax1.xaxis.set_major_locator(ticker.IndexLocator(base=3, offset=2)) # major ticks at 2,5,8,11, etc. other options exist
ax1.set_ylabel('Mag', fontsize=14)  # y-label
ax1.set_xlabel('B-V', fontsize=14)  # x-label
ax1.set_yscale('log')  # log_10 scale on y-axis
ax1.set_title(r'$\sigma_{15}$ Some title')  # A plot title. The r means a raw string (for latex \)
ax1.text(60, .025, r'$\mu=100,\ \sigma=15$')
ax1.hist(x,30,facecolor='b',range=[0.9,1.1]) #histogram of x, 30 bins in [0.9,1.1]
ax1.annotate("",xy=(0.54,0.243),xytext=(0.43,0.243),
ax1.scatter(x,y,c=z,s=sz,vmin=0.8,vmax=1.2) #scatter plot, with colors defined by z, size by sz, over range 0.8-1.2

ax3=ax2.twinx()         #Add a labeled axis on the right that differs from the axis on the left
ax3.set_ylim([0,0.4])
ax3.set_yticks([0.04545,0.0703,0.1005,0.2047,0.3653]) #custom tick locations
ax3.set_yticklabels(["0","1","2","5","10"])           #and custom tick labels
ax3.text(1.14,0.26,r'$\tau_{(^{13}CO)}$',rotation=90,fontsize=16)

plt.subplot(325)  #3 rows, 2 columns of plots; access plot number 5 of these
plt.subplots_adjust(wspace=0.5, hspace=0.5) #extra space between the subplots
plt.subplot2grid((6, 8), (0, 3), rowspan=2, colspan=2)  # More flexible subplots. (6,8) is grid size
# (0,3) is starting location
plt.axis([-0.2,0.8,11,5])  #another way to do limits on x- and y-axes
plt.savefig('output.png') # Save to a file
plt.show()         # Displays the plot
fig.add_subplot(grid[-1:],sharex=ax1) # add subplot at bottom, then plt.xlabel, plt.plot etc to plot there

# colors: (r)ed, (g)reen, (b)lue, (w)hite, blac(k), (c)yan, (m)agenta, (y)ellow,
#  hundreds of others like 'darkblue', 'forestgreen', 'brown', etc.
# linestyles:  : dotted; - solid; -- dashed; -. dot-dashed;  possible to make your own
# point types: o filled circle; d or D: diamonds; s squares; + : plus ; * star ; v and ^ triangle

# Polar plot example plotting (r,theta) and using other columns for point size and labels
fig=plt.figure(figsize=(12,12))
ax=fig.add_subplot(111,projection='polar')   # make a polar plot
# An example using sorted data in a file to control point type
# ax.scatter allows you to use a column to define size and color but not point type. this is a workaround
ax.scatter(np.pi*data['theta'][0:18]/180.,data['r'][0:18],s=data['sz'][0:18],marker='x',color='black')
ax.scatter(np.pi*data['theta'][18:40]/180.,data['r'][18:40],s=data['sz'][18:40],marker='o',color='red')
ax.scatter(np.pi*data['theta'][40:42]/180.,data['r'][40:42],s=data['sz'][40:42],marker='d',color='darkgreen')
ax.scatter(np.pi*data['theta'][42:44]/180.,data['r'][42:44],s=data['sz'][42:44],marker='*',color='black')
ax.set_rmax(3.2)  # radial maximum for a polar plot
ax.set_rticks([0.5,1,1.5,2,2.5,3.0])
ax.set_rlabel_position(313)  # Angle in degrees where the labels run
ax.grid(True)
ax.set_theta_zero_location("S") # Where to put zero degrees on the plot
# Use a label in a column to identify each point
for iv in range(np.size(data['r'])) :
plt.text(np.pi*data['ltheta'][iv]/180.,data['lr'][iv]+0.1,data['name'][iv],fontsize=10)
# Example of a legend
ax.legend(['1-6','6-25','26-50','>50'],fontsize=14,loc='upper right', bbox_to_anchor=(1.15, 1.0))
plt.title('Number of O Stars in Nearby Associations',fontsize=20)
plt.show()


---------------Animations--------------------
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,10))
ax.set_ylim(0,11)
ax.yaxis.grid(False)
ax.xaxis.grid(False)
ax.set_yticklabels([])
yrs=120           #length of simulation in years
tinterval=15./365.25  #time step in years
d=np.int(np.round(yrs/tinterval)) #number of frames in the animation
pjup=11.86     #Jupiter period
psat=29.46     #Saturn period
thear=(2*np.pi)*tinterval #Theta changes per interval
thjup=thear/pjup
thsat=thear/psat
t0e=0.*(np.pi/180.) #Initial thetas
t0j=0.*(np.pi/180.)
t0s=0.*(np.pi/180.)

ax.plot(,,'yo',ms=16) #Sun

th=np.arange(0,365,5)*np.pi/180. #plot orbits
t0=0
t1=np.pi*242.7/180.
t2=np.pi*125.4/180.
t3=np.pi*8.15/180.
t4=np.pi*250.9/180.
t5=np.pi*133.6/180.
t6=np.pi*16.3/180.
ax.plot([t0,t0],[0,10.3],'k-.',lw=0.5)
ax.plot([t1,t1],[0,10.3],'k-.',lw=0.5)
ax.plot([t2,t2],[0,10.3],'k-.',lw=0.5)
ax.plot([t3,t3],[0,10.3],'k-.',lw=0.5)
ax.plot([t4,t4],[0,10.3],'k-.',lw=0.5)
ax.plot([t5,t5],[0,10.3],'k-.',lw=0.5)
ax.plot([t6,t6],[0,10.3],'k-.',lw=0.5)
ax.text(0, 10.25, 'C-0')

#animated parts to plot
jupiter, = ax.plot([],[], 'o',mfc='indianred',ms=10)
saturnr, = ax.plot([],[], 'k_',ms=18) #add a ring to saturn as it goes around
saturn, = ax.plot([],[], 'ko',ms=9)
earth, = ax.plot([],[], 'bo',ms=7)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
c1_text = ax.text(242.7*np.pi/180., 10.8, '')
c2_text = ax.text(125.4*np.pi/180., 10.6, '')
c3_text = ax.text(8.2*np.pi/180., 10.25, '')
c4_text = ax.text(250.9*np.pi/180., 10.8, '')
c5_text = ax.text(133.6*np.pi/180., 10.6, '')
c6_text = ax.text(16.3*np.pi/180., 10.25, '')

# initialization function with blank data
def init():
jupiter.set_data([],[],)
saturn.set_data([],[],)
saturnr.set_data([],[],)
earth.set_data([],[],)
time_text.set_text('')
#returns ax.plot commands
return  jupiter, saturn, saturnr, earth, time_text, c1_text, c2_text, c3_text, c4_text, c5_text, c6_text

# animation function.  This is called sequentially
def animate(i):
aa=i*tinterval
if((aa>19.7) & (aa<19.9)):
c1_text.set_text('C-1')
if((aa>39.6) & (aa<39.8)):
c2_text.set_text('C-2')
if((aa>59.4) & (aa<59.6)):
c3_text.set_text('C-3')
if((aa>79.5) & (aa<79.7)):
c4_text.set_text('C-4')
if((aa>99.2) & (aa<99.4)):
c5_text.set_text('C-5')
if(aa>119.0):
c6_text.set_text('C-6')

the = i*thear+t0e
thj = i*thjup+t0j
ths = i*thsat+t0s

th = [the,thj,ths]
earth.set_data(the, Re)
jupiter.set_data(thj,Rj)
saturnr.set_data(ths,Rs)
saturn.set_data(ths,Rs)
time_text.set_text('%.2f years' % aa)
# return ax.plot commands
return jupiter, saturn, saturnr, earth, time_text, c1_text, c2_text, c3_text, c4_text, c5_text, c6_text

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = FuncAnimation(fig, animate, init_func=init, frames=d, interval=20, blit=True)
anim.save('test.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()


---------------Displaying and Altering Images--------------------
plt.figure(figsize=(20,10))  #set the size
plt.axis('off')       # shuts off drawing of x,y axes
plt.imshow(imname[20:42,29:48],origin='lower',vmin=-120,  #imname is np.array, vmin and vmax are scale limits
vmax=20,cmap='jet',interpolation='nearest',aspect='equal') # cmap is colormap, aspect=equal gives square pixels
plt.colorbar()

plt.imshow(blurredz[20:42,29:48],origin='lower',vmin=-120,vmax=20, # OR try using 'extent' to fix the labels
cmap='jet',interpolation='nearest',       # Fits image [x,y] --> np.array [y,x] via astropy fits.open
aspect='equal',extent=(29,48,20,42))      # But imshow displays np.array [y,x] as [x,y] again (I think)
# Note, labels x,y flip for extent.
# If you were to put [2900,4800,2000,4200] for extent it would only
# change the labels

blurredz = gaussian_filter(data3, sigma=1.5, mode='nearest') # smooth an image with a gaussian filter


---------------Strings--------------------

# Assigning, concantating strings
c="hello"
ss=c+" "+c      # ss = "hello hello"

# uppercase and lowercase
ss= "HEllo"
ss.upper()     #would be HELLO
ss.lower()     #would be hello
ss.isupper()   #returns False because HEllo is not all uppercase

# string length
ss="lasjhf"
len(ss)       # equals 6

# substrings and positions in a string, replacing strings
ss=Hi there"
ss[1:4]      # equals 'i th' because the index starts at 0
ss[len(ss)-3:len(ss)-1]  # 'ere' the last three characters
ss[-2]       # second to last character, 'r'
ss[-2:]      # second to last character to the end of the string, 're'
ss[::-1]     #reverses the string, 'ereht iH'
ss.index("h")    # equals 4 in this case
ss.index("her") # also equals 4, the location of the beginning of 'her'
ss.replace("her","eak")  # gives 'Hi teake'

rstr="10:43:24.5800 -59:36:03.300"
ram=float((rstr.split(" ",1)).split(":",3))  # sets ram=43.0

#         Formatted print statements
print("{0:5.3f} {1:3d}".format(var,num))  #the 0: is the first argument (var),
#    1: is the second argument (num)
print('Input filename: ',end="")          # stops the carriage return
print("{0:02.0f}:{1:02.0f}:{2:07.4f}".format(rah,ram,ras)) # does leading zeros for RA, e.g. 09:08:04.233
# Be careful, as this format rounds up or down

print('12'.zfill(5))             # for leading zeros, prints 00012
if(np.mod(ix,1+np.int(ixbig/20))==(np.int(ixbig/20)-1)) :
print("{0:s}Doing the composite median {1:2d}%".format("\r",np.int(100*ix/ixbig)),end="")
# prints out progress of task in a loop, every 5% or so,
# overwrites the 5% in place with \r format


-----I/O into a program, open and write to output files -------

name=input("Input filename: ")                    # read name from terminal
oname="out_"+name               # make output filename
fin= open(name,'r')             # open input filename to read  'a' would be to append
fout= open(oname,'w')           # open output filename to write (overwrite) an entirely new file
for line in fin:                # loop until run out of lines
line=line.strip()            # get rid of invisible characters like \n
columns=line.split()         # split the fields of the line
var1=float(columns)      # columns and columns are strings
var2=float(columns)
fout.write('{0:9.4f} {1:9.4f}\n'.format(var1-var2,var1*var2))
fin.close()
fout.close()

tfile=open('times', "r")
while line:              # read in subsequent nonzero lines
print(line)
tfile.close()

name.write("\n adding this line to file") # for appending to existing files

a,b = float(input("Type in two values: ").split())  # input 2 real values

from astropy.io import ascii
#junk.txt has a line of header: name, B-V, and mag, and then data
data.more()                      #To interactively page through the table
data.info                        #shows what columns are what and their type
data.colnames                    #list column names
data.rename_column('dia','D')    #rename a column
data['newcol'] = [2,4,5,2,3,5]   #add a column named 'newcol'
data.remove_column('newcol')     #remove a column
data['dia'].format = '7.3f'      #changes format of column 'dia'
data['dia','B-V']                #a table with only these two columns
data1=data[np.where((data['long']<20) & (data['dia']>10))] #extract subset
np.where(a>4)  # returns a 'tuple', with the indices where the condition is true
# that's why np.where is inside the [ ] in the above subset example
# most of the time you probably don't want to be doing math on np.where
# because you'll be doing math on the indices, not on the values
np.where(a>4, a, 5*a) # returns 'a' when a>4, and 5*a otherwise, each in a tuple
np.reshape(a,np.product(a.shape))  # converts a 2D array to a 1D array


------- FITS I/O----------
                       #Reading
from astropy.io import fits
hdul=fits.open('test1.fits') # opens the HDU list
data1=np.transpose(hdul.data)   # the data  WARNING: astropy transposes the array
# so x is y unless you include the np.transpose
# do the transpose here and the array behaves like the fits file
hdr1['NAXIS1']  # value of the Fits keyword NAXIS1
ixlim=data1.shape  # x-dimension of array
iylim=data1.shape  # y-dimension of array
print(data1[2:4,5:6]) #prints value of (3,6) and (4,6)

#Writing
hduout=fits.PrimaryHDU(np.transpose(dataout).astype(np.float32)) #puts array dataout into a Fits HDU. Default type is 'double'
# note use of transpose to write out the array
hduout.writeto('test2.fits', overwrite=True) # Writes the HDU to an output file


------- if, for, while statements; looping ----------
flag1=True
flag2 = False
if flag1 and flag2:    # you can also use 'or'.  DONT FORGET the :
print("ok")
elif flag1 and not(flag2): # use of 'not' function
print("maybe ok")
else:
print("not ok")
quit()          #end program if not ok

if num1 >= num2 and num2 != num3:  # For numbers
print(num1,num2)
print("done with loop")

if ((a[i,j] == 0) | (a[i,j] == 2) & (t==8)) :  # note use of ==, |, &
print(a[i,j])

for i in range(0,upper) :  # 0,1,2,...upper-1  i.e., the indexed values for an array of length upper
print(i)

vel=np.zeros(data1.shape)
for i in range(data.shape):
vel[i] = i*15.-150             #set some values to array vel[i]

for iv in list(range(8, 12)) :
print(iv-2.3)

for iv in range(np.size(vs)) :
print(iv-2)

i=0
while (i< len(data['ra'])):
print(i)
i+=1
print("done with loop")

while True:
try:
a=int(input("Input an integer: "))
break                              # break out of what would be an infinite while loop
except ValueError:                  # If they didnt input an integer
a=0
print("You need an integer here")
print(a)

name=np.zeros(length, dtype=np.int)      # initialize a 1D array with zeros
name=np.zeros((nx,ny), dtype=np.float32) # initialize a 2D array with zeros
np.shape(a)                              # gives the size of the array a
np.shape(a)                           # the length of the first axis of the array
cgrid=np.empty((iybig,ixbig,nfiles))     # initializa an array as empty
cgrid[::]=np.NAN                         # ...and then fill it with NANs; useful for masking


-----Lists-------
list=["bobby",23,"billy",55.3]
list         # this is "billy"
list[-1]        # this is 55.3, go from the end
list[2:]        # starts at element 2 and goes to the end:  ["billy",55.3]
list1=["a","b"]
list2=[1,2]
list1.extend(list2) # glues lists together: ["a","b",1,2] and overwrites list1
list2.append(5)     # list2 is now [1,2,5]
list2.index(5)      # returns 2, the index location of element '5'
list2.count(5)      # number of times '5' occurs in the list
list2.insert(1,12)  # list2 is now [1,12,2,5]
list2.remove(12)    # removes first occurrence of 12
list2.pop()         # remove last element from list
list.sort()         # sorts the list
list.reverse()      # reverse sorts the list
list2.clear()       # clears the list
list3=list1.copy()  # copies a list
list4=[[0,3,5],[9,4,2],]  # list4 = 9 in this case

array = (2,4)       # this is a 'tuple' and can never be changed
# so you can't change the 4 to a 5, for example (seems stupid)


---------------Dictionaries--------------------
dateno = {         # this defines the dictionary with a key/value
"jan": "0",        # e.g. mar/59
"feb": "31",
"mar": "59",
"apr": "90",
"may": "120",
"jun": "151",
"jul": "181",
"aug": "212",
"sep": "243",
"oct": "273",
"nov": "304",
"dec": "334",
}
while (True):
a,b = input("Date (e.g. Feb 22): ").split()  #input multiple values on one line
print("screwup")
else:
print(int(aa)+int(b))


---------------Error Handling--------------------

while True:
try:
a=int(input("Input an integer: "))
break                              # break out of the infinite while loop
except ValueError:                  # If they didnt input an integer
print("You need an integer here")
print(a)


---------------Classes --------------------
class exam:              # Define a new class of objects
def __init__(self,question,answer):   # exam objects consist of a question and an answer
self.quest = question   # e1.quest will refer to question part of the exam element e1
self.ans = answer      # e1.ans will refer to answer part of the exam element e1

# make a list of questions
questions = [
"Is the sky blue?\n(a) yes\n(b) no\n\n",
"Does Trump have yellow hair?\n(a) yes\n(b) no\n\n",
"What color are the Rams jerseys?\n(a) purple\n(b) pink\n(c)yellow and blue\n\n",
]

key = [                  #the list 'key' contains exam objects
exam(questions,"a"),
exam(questions,"a"),
exam(questions,"c")
]

score=0
for i in key:            #each i is an exam object
score += 1
print("You got " + str(score) + "/" + str(len(key)) + " Correct")


---------------Subroutines (functions)--------------------
def yrs_since_HS(age):
return(age-18)
print(yrs_since_HS(55)) # gives 37.  Input was age, returned value is age-18


---------------Modules and Importing--------------------
import numpy as np
np.pi                 # 3.14159265...
pip install module
pip uninstall module


---------------Interact with the operating system--------------------

# To execute a command from the OS and print the results to STDOUT, precede with a '!'
!pwd

# To call a command from the OS:
from subprocess import call
call(["/bin/ls","-lt"])  # if successful it simply exits with code '0'

# To call a command from the OS that uses pipes:
import subprocess
p1=subprocess.Popen(["/bin/ls"],stdout=subprocess.PIPE)   #first command
p2=subprocess.Popen(["/usr/bin/awk","END{print NR}"],stdin=p1.stdout,stdout=subprocess.PIPE) #pipe to second command
p1.stdout.close() #apparently so the order of finishing doesn't matter
output=p2.communicate() #set the result from the pipe (a byte string) to 'output'
print(output.decode('ascii')) #print out this byte string as ascii

# To write python commands to a file, and then execute that file in a double loop
# Note if you forget the () in g.close the python code will skip one command
#  and you get a real mess, with only 5 output files instead of 6 and mixed up contents
#  That's because the file remains open, and only closes when it reloops
from subprocess import call
import numpy as np
vs=(20,30)
denz=(10,20,30)
for iv in range(np.size(vs)) :
for id in range(np.size(denz)) :
fout=open("jj",'w')
fout.write('{0:10.2f}  {1:4.2f} VS DEN\n'.format(vs[iv],denz[id]))
fout.close()

g=open("doit",'w')
g.write('cp -f jj test_{0:2d}_{1:4.2f}.txt\n'.format(int(vs[iv]),denz[id]))
g.close()   # Don't forget the ()!!!
call(["/bin/tcsh","doit"])


---------------Jupyter Notebook--------------------

# To start up Jupyter:
jupyter notebook  # will create a .ipynb file (the notebook) in the current directory
# uses  ~/.jupyter/nbconfig/notebook.json for shortcuts
# These are vi-like
# Jupyter has an insert mode and an escape mode like vi
# Don't be afraid to restart the kernel if it gets confused

jupyter nbconvert --to python nbname # Converts nbname.ipynb to nbname.py, a python code file
# you can run via 'python nbname.py'
# Note, jupyter notebooks are really json files, meaning they
# are text files that can be edited