<jupyter_start><jupyter_text>Computational Science Best Practices <NAME> CATMS Lunch Talk, 14 Sep 2017  Overview- Organization- Write programs for people- Automate Repetitive Tasks- Keep a record of what you did- Make small changes- Use Version Control, starting yesterday- DRY- Debugging- Optimization / Profiling- Comment on CommentsMuch of this talk is inspired by / taken from https://arxiv.org/pdf/1210.0530v3.pdfIt is supplemented with my opinions, they may or may not be relevant for youBut they probably are  Organization- You're on a computer, use it!- Bad:    - /data        - test.dat        - test2.dat        - test_01_feb.dat        - temp.dat        - data.dat        - data1.dat        - water.dat        - uhh.dat        - minecraft.exe    - /programs        - project1.exe        - project2.exe        - project1_working.exe        - project1_no_actually_working.exe        - All of my hatred:    - /this is a folder        - here is a file lol.xlsx       - My Recommendation    - /Project1        - /src            - /blah.cpp        - /data            - data.json        - /doc        - /papers-presentations            - blah.tex            - blah.ppt        - /images            - blah.eps  Write Programs for people- Meaningful variable names (tt vs times)- Encapsulation (small, digestible chunks of code)- Self-documenting<jupyter_code>from math import sqrt
def f(x, y, z):
    return [x / sqrt(x**2 + y**2 + z**2), y / sqrt(x**2 + y**2 + z**2), z / sqrt(x**2 + y**2 + z**2)

from math import sqrt
def normalize(vec):
    norm = sqrt(sum(i**2 for i in vec))
    return [i / norm for i in vec]

print(f(1, 2, 3))
print(normalize([1, 2, 3]))
def integrals():
    #The total size of the basis
    N = na + nb + nd

    #The starting point for each basis center
    N0a = 0 
    N0b = na 
    N0d = na + nb 

    #Initialize the basis arrays, Tbas is the kinetic energy operator acting on basis functions
    bas  = np.empty([N,Z])
    Tbas = np.empty([N,Z])

    #The basis centered at a 
    for i in range(N0a, N0b):
        for j in range(0,Z):
            bas[i][j]  = (2**(i-N0a)*factorial(i-N0a))**(-0.5) * (wa/pi)**0.25 * np.exp(-0.5*wa*(x[j]-a)**2) * sp.hermite(i-N0a)(wa**0.5*(x[j]-a)) 
            Tbas[i][j] = bas[i][j] * (((i-N0a) + 0.5) * wa - 0.5 * ka * (x[j]-a)**2)

    #The basis centered at b 
    for i in range(N0b, N0d):
        for j in range(0,Z):
            bas[i][j] = (2**(i-N0b)*factorial(i-N0b))**(-0.5) * (wb/pi)**0.25 * np.exp(-0.5*wb*(x[j]-b)**2) * sp.hermite(i-N0b)(wb**0.5*(x[j]-b)) 
            Tbas[i][j] = bas[i][j] * (((i-N0b) + 0.5) * wb - 0.5 * kb * (x[j]-b)**2)

    #The basis centered at d 
    for i in range(N0d, N):
        for j in range(0,Z):
            bas[i][j] = (2**(i-N0d)*factorial(i-N0d))**(-0.5) * (wd/pi)**0.25 * np.exp(-0.5*wd*(x[j]-d)**2) * sp.hermite(i-N0d)(wd**0.5*(x[j]-d)) 
            Tbas[i][j] = bas[i][j] * (((i-N0d) + 0.5) * wd - 0.5 * kd * (x[j]-d)**2)

    #Compute the necessary integrals 
    T  = np.empty([N,N])
    U  = np.empty([N,N])
    H  = np.empty([N,N])
    S  = np.empty([N,N])

    for i in range(N):
        for j in range(N):
            S[i][j] = np.trapz(bas[i]*bas[j],x)
            T[i][j] = np.trapz(bas[i]*Tbas[j],x)
            U[i][j] = np.trapz(bas[i]*bas[j]*V,x)
            H[i][j] = T[i][j] + U[i][j]<jupyter_output><empty_output><jupyter_text>Automate Repetitive Tasks / Keep a recordMy recommendation: use a scripting language and save the scripts!- bash- python- perl (if you're a heathen)- others?Can use for running jobs, data analysis, organizing files based on criteria, writing papers even  Make small changes and use version control, starting yesterdayOptions:    - git     - svn- Keeps track of all the changes to your files throughout history- Keeps folders tidy, and makes sure you never lose progress- Basically essential for collaborative work (multiple people working on different parts of same project)- Easy to remotely backup all necessary filesThis allows you to make small changes, and only include them into your projectwhen you're sure it's working.  No more:    - /project        - file1        - file2            - /project_15_sep_2017        - file1        - file2            - /project_18_jan_2018        -file1         -file2        -file3 Note:We all get free academic github accounts (unlimited private repositories)SCS also has a VCS server  Don't Repeat Yourself (DRY)- Often we need to do the same thing to different objects- Most of the time the common behavior can be separated out and reused- Ensures consistency- Saves time<jupyter_code>def integrals():
    '''stuff...'''
    #The basis centered at a 
    for i in range(N0a, N0b):
        for j in range(0,Z):
            bas[i][j]  = (2**(i-N0a)*factorial(i-N0a))**(-0.5) * (wa/pi)**0.25 * np.exp(-0.5*wa*(x[j]-a)**2) * sp.hermite(i-N0a)(wa**0.5*(x[j]-a)) 
            Tbas[i][j] = bas[i][j] * (((i-N0a) + 0.5) * wa - 0.5 * ka * (x[j]-a)**2)

    #The basis centered at b 
    for i in range(N0b, N0d):
        for j in range(0,Z):
            bas[i][j] = (2**(i-N0b)*factorial(i-N0b))**(-0.5) * (wb/pi)*0.25 * np.exp(-0.5*wb*(x[j]-b)**2) * sp.hermite(i-N0b)(wb**0.5*(x[j]-b)) 
            Tbas[i][j] = bas[i][j] * (((i-N0b) + 0.5) * wb - 0.5 * kb * (x[j]-b)**2)

    #The basis centered at d 
    for i in range(N0d, N):
        for j in range(0,Z):
            bas[i][j] = (2**(i-N0d)*factorial(i-N0d))**(-0.5) * (wd/pi)**0.25 * np.exp(-0.5*wd*(x[j]-d)**2) * sp.hermite(i-N0d)(wd**0.5*(x[j]-d)) 
            Tbas[i][j] = bas[i][j] * (((i-N0d) + 0.5) * wd - 0.5 * kd * (x[j]-d)**2)
            
    '''more stuff...'''
def HO_basis(n, center, freq, x_range):
    x = x_range - center
    norm = sqrt(2**(n) * factorial(n)) * (freq / pi)**(0.25) 
    return norm * np.exp(-0.5 * freq * x**2) * sp.hermite(n)(sqrt(freq) * x)
    
def integrals():
    '''stuff...'''
    #The basis centered at a 
    for i in range(N0a, N0b):
        for j in range(0,Z):
            bas[i][j]  = HO_basis(n=(i - N0a), center=a, freq=wa, x_range=x)
            Tbas[i][j] = diff(0.5 * diff(2, bas[i][j]))

    #The basis centered at b 
    for i in range(N0b, N0d):
        for j in range(0,Z):
            bas[i][j] = HO_basis(n=(i - N0b), center=b, freq=wb, x_range=x)
            Tbas[i][j] = diff(0.5 * diff(2, bas[i][j]))

    #The basis centered at d 
    for i in range(N0d, N):
        for j in range(0,Z):
            bas[i][j] = HO_basis(n=(i - N0d), center=d, freq=wd, x_range=x)
            Tbas[i][j] = diff(0.5 * diff(2, bas[i][j]))<jupyter_output><empty_output><jupyter_text>Debugging / Testing- You will make mistakes, so this is unavoidable- All languages have debuggers- Step through code as you go, checking on all variables, etc- Don't have to know where the problem is!!!!- Visual debuggers exist, and I love them- Stop this:    - code1    - print(1)    - code2    - print(2)    - code3    - print(3)    - code4    - print("done")Should write testing for code to ensure it's functioning as expectedLearn to love assert<jupyter_code>def Kelvin_to_Celsius(temp):
    assert(temp >= 0), 'No temperatures exist below 0K'
    return(temp - 273.15)

print(Kelvin_to_Celsius(73))
print(Kelvin_to_Celsius(-5.2))<jupyter_output>-200.14999999999998<jupyter_text>Premature optimization....- Write correct code first, then make it faster when needed- PROFILE, maybe the slow part is not what you expect    - caveat: Many numerical codes are very obvious (8 fold loop maybe?)  - Recommendation: write code first in a high-level language, then port to lower level if needed<jupyter_code>import cProfile
import numpy as np

N = 50
mat1 = np.asarray([[1.5*i + j for i in range(N)] for j in range(N)])
mat2 = np.asarray([[2.5*i + j for i in range(N)] for j in range(N)])

def func(mat1, mat2):
    exp = np.exp(mat1)
    log = np.log(mat2)
    result = np.dot(exp, log)
    sum_rows = np.sum(result, axis=0)
    normed = sum_rows / np.norm(sum_rows)
    return normed<jupyter_output><empty_output><jupyter_text>Comments- Comments should explain reasoning- If a comment is explaining what a block of code does, should probably be a function- Comments explaining what a line of code is doing are distracting and unnecessary<jupyter_code>def increment_by_one(i):
    """Increments i by one
    
    Uses the plus operator to make i bigger than it used to be
    by the value of one.
    """
    i = i + 1 # increment i by 1
    return i
    
    
def normalize(self, vec):
    """Return a normalized version of given vector.
    
    Vec is not modified in place.
    :param vec: 1D Vector of values along the self.coords.
    :type vec: np.ndarray
    """
    return 1.0 / math.sqrt(self.overlap(np.conjugate(vec), vec))<jupyter_output><empty_output><jupyter_text>A Note on logfiles- Logging your results is important- It's nice to have a logfile you can read- It's usually nicer to have a logfile computers can read- In many cases, you can get both- Why are you writing your own parser?- Standardized formats are more efficient and easier to work withDepending on data type, I like .json files, but HDF5 is good when you need heavy duty. xml probably is good, maybe others.<jupyter_code>import json

mydata = {'n_grid_points': 50,
          'x_min' : -10.0,
          'x_max' : 10.0,
          'integral_value' : 1.36         
         }

s = json.dumps(mydata)
s = s.replace('{', '{\n').replace('}', '\n}').replace(',', ',\n')
print(s)

loaded = json.loads(s)
print(loaded)<jupyter_output>{
"n_grid_points": 50,
 "x_min": -10.0,
 "x_max": 10.0,
 "integral_value": 1.36
}
{'n_grid_points': 50, 'x_min': -10.0, 'x_max': 10.0, 'integral_value': 1.36}
{'fraction_non_alphanumeric': 0.1559680597324484, 'fraction_numerical': 0.03256248055584362, 'mean_word_length': 2.632391713747646, 'pattern_counts': {'":': 4, '<': 24, '<?xml version=': 0, '>': 25, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '468283', 'license': 'MIT', 'n_tokens_mistral': 3736, 'n_tokens_neox': 3495, 'n_words': 1053, 'path': 'doc/Sep_15_2017_Best_Practices.ipynb', 'repo_name': 'CATMS/bootcamp', 'size': '9648'}