Introduction to Scientific Computing
- Course Numbers CHM 4932-004/6938-010, CRN: 20438/20439
- Credit Hours: 3
- Meeting Dates: January 8 - April 25, 2018
- No Class Jan. 15 or Mar. 12-18
- Meeting Times: Mon. and Wed., 12:30-13:45 in SCA 222
- Midterm Project Presentations: Mon.-Wed., Mar 5 and 7.
- Final Project Presentations: Wed., May 2, 10:00am – 12:00pm (scheduled final date)
- Lab Session & Office Hours: Fri. 12:30-13:45 in IDR 214 (or SCA 222 by announcement)
$ for((i=1;i<=4;i++)); do cal $i 2018; done ____ January 2018 ____ _____ March 2018 ____ | Su Mo Tu We Th Fr Sa | | Su Mo Tu We Th Fr Sa | 1 2 3 4 5 6 | | 1 (2) 3 1 | 7 8 9*10 11 12 13 | 9 | 4 (5) 6 (7) 8 9 10 2 | 14 x 16*17 18 19 20 | | x x 3 | 21*22 23 24 25 26 27 | 10 | 18*19 20 21 22 23 24 4 | 28*29 30 31 | 11 | 25*26 27 28 29 30 31 ___ February 2018 ____ _____ April 2018 ____ | Su Mo Tu We Th Fr Sa | | Su Mo Tu We Th Fr Sa | 1 2 3 | 12 | 1 *2 3 4 5 6 7 5 | 4 *5 6 7 8 9 10 | 13 | 8 *9 10 11 12 13 14 6 | 11*12 13 14 15 16 17 | 14 | 15*16 17 18 19 20 21 7 | 18*19 20 21 22 23 24 | 15 | 22*23 24(25)26 27 28 8 | 25*26 27 28 | | 29 x Week Topic 1 Calculator Form 2 Pseudocode 3 Plotting 1 (Excel) 4 Plotting 2 (Matplotlib and linear spaces) 5 Abstraction 1 (advanced functions and oo method) 6 Abstraction 2 (Booleans and logic) 7 2D transforms and iterated maps 8 Optimization 1 (Newton-Rhapson) 9 midterm project presentations 10 Optimization 2 (MD Energy) 11 ODE and SDE algorithms 12 Statistics 13 Complexity 14 Circe and Linux crash course 15 Compile & install process, ctypes
Mathematical models in the natural sciences take many shapes and forms. The moment they become more complicated than tracking a few variables over time, we have to put away our Excel spreadsheets and make use of a solid foundation in computational science. The course establishes this foundation by introducing data structures and algorithms used in everyday scientific computing using examples in the python scripting language. By the end of the Semester, students will be able to navigate the Linux command-line to go from a mathematical model to a numerical solution. Best practices for working with computer code and visualizing data will also be covered.
- bash (shell)
- python (scripting)
- stack & register-based machine languages: overview of C and Fortran syntax & compiling with gcc, ‘thread safety’
- Shared libraries & language inter-operability (FFI)
- Horner’s algorithm for polynomials
- Newton’s root finding algo.
- Overview of general optimization algo-s (example: linear & nonlin. least-squares)
- Numerical integration (scripted, plus gnu ODE)
- Communication (how HTTP get/put works)
- Data Structures
- Linked Lists
- Trees (file system hierarchies)
- Arrays (dense & sparse) example: vectors & rotations / transpositions, densities and difference operators
- Presentation & visualization
- Flat CSV files (and excel)
- Working with binary data
- Python matplotlib (2D images)
- Binning histograms & weighted averages
- Best Practices
- High-level code design (modularity, interface specifications, unit testing, etc.)
- Code audit: open-source libraries, ex. Gnu Scientific Library, qsort in libc, Boost
- Source code versioning (example: git or mercurial)
- gdb & execution profiling
Grading Scale (out of 400 points possible)
- (280 U / 330 G) Points Possible:
- 9 homework assignments (10 points per homework)
- 14 in-class assignments / quizzes (10 points per week)
- 2 projects / competitions (50 points per project)
- Undergraduates must complete 1 project, graduates must complete 2
- Important Rules
- For maximum benefit, homework should be done on-time. However, homework will be allowed to be re-submitted at any time for full credit, but must present a unique solution to previous work from yourself or others.
- All work must cite co-group members (if any) and the source of any code snippets used online.
Homework points will only be awarded for correct, working code and correct algorithmic analysis, while the homework problems will become progressively tougher. Undergraduates should be able to pass the course with an A if they have correctly answered most of the homework, appear at class prepared and on-time, and made a strong attempt at one of the projects. Graduate students should both perform well on the homework and successfully complete one or two projects. The re-grading for homework allows making up for missed assignments, and encourages returning to core concepts at a later point in the class. Students who have put in the studying effort, completed significant portions of the homework and projects both in- and out of class have thus demonstrated introductory programming proficiency.
- Midterm Project Description Media:Midterm2018.pdf (2 Mar., 2018)
- Midterm Projects:
- Final Project Description (25 Apr., 2018)
Textbooks and Resources
- Beginning Python, Hetland, 2005 (eBook avail. from USF Library)
- Python Algorithms, Hetland, 2010 (eBook avail. from USF Library)
- The Architecture of Open Source Applications
- Python Documentation
- Google Introductory Python Notes
- Tao of Programming
Other (possibly useful) materials:
- The Jargon File
- Chapters 7&8 of Modelling and Simulation, Birta and Arbez, Springer 2007 (eBook avail. from USF Library).
- Chapters 4-6 of Models and Algorithms for Intelligent Data Analysis, Runkler, Springer 2012 (eBook avail. from USF Library).
- Guide to Scientific Computing in C++, Francis and Whiteley, Springer, 2012. (eBook avail. from USF Library)
- Extra info. On the basics:
- Chapters 3 and 7 of Mathematics in Computing, Regan, Springer 2013. (eBook avail. from USF Library)
- Instant feedback graphical programming with Khan Academy
- Local PyBulls Meetings
- Python Libraries Useful for Geospatial Problems
- Python Libraries Useful for Computational Chemistry (very small list at present)
- CodeWars Recommended for building experience with practice problems.
- CodeCombat A very visual way to interact with your python code.
- Exercism A command-line client with a similar learn-through-challenge philosophy.
- repl.it A read-eval-print loop for several computer languages.
- CppReference Documentation with an online evaluator for C++ code.
- Introduction Notes
- Function Syntax and For-Loop Syntax cheat-sheets.
- Plotting Maps
- SIAM Notes on Mathematical Modeling
- The simple, yet amazing harmonic oscillator
Review Notes on Sequences, list comprehensions and arrays
- Say you want to make a plot of the LJ function on the domain x \in [0.9, 3]
- 1. You could use the following two list comprehensions to be specific:
x = [ i*(3-0.9)/100.0 + 0.9 for i in range(101) ] y = [ 4*(xi**-12 - xi**-6) for xi in x]
- 2. You could implement the LJ function as a python function and use MAP,
ix6 = x**-6 # the factoring makes this calculation more efficient return 4*ix6*(ix6 + 1.0)
y = map(LJ, x)
- 3. You could work on arrays, where addition, multiplication, and powers are elementwise:
ax = arange(101)*(3-0.9)/100 + 0.9 ay = LJ(ax)
- Mind the difference in type!
- 4. You could totally cheat, and use numpy's builtin functions.
ax = linspace(0.9, 3.0, 101)
- Exercise: replicate the steps above to test the finite difference computation
- of the derivative against the analytical derivative of LJ.
[ (y[i+1]-y[i])/(x[i+1]-x[i]) for i in ...
import matplotlib.pyplot as plt plt.plot(x, y) plt.show() </source>
Review Notes on Functions with Function arguments (functionals)
- map : (a -> b), [ a ] -> [ b ]
- map: map(str, [1,2,3]) --> [ "1", "2", "3" ]
- ^ function (a -> b)
- map(int, map(str, [1,2,3]) )
- - helpful for parsing
- x = "1 2 3 4 5"
- x.split() # ["1", "2", "3", "4", "5"]
- map(int, x.split()) # [1,2,3,4,5]
str([1,2,3]) # "[1, 2, 3]" (bad) map(str, [1,2,3]) # ["1", "2", "3"] (good)
- filter: (a -> Bool), [a] -> [a]
- Example filter:
- '1' -> no
- ' ' -> no
- 'a-z' -> yes
- 'A-Z' -> yes
def f(x): # str -> Bool
if x >= 'a' and x <= 'z': return True if x >= 'A' and x <= 'Z': return True return False
filter(f, " 1 a b 2 \n")
- reduce: (a,b -> a), [b], a -> a
sum([1,2,3]) # -> 6 def add2(x,y): # x : int, y : int
reduce(add2, [1,2,3], 0)
- For an example where a != b:
def append(x,y): # x : [int], y : int
return x + [y]
reduce(append, [1,2,3], ) # start with an empty list
def nrhap(f, fp, x0, tol = 1e-6):
y = f(x0) iter = 0 # guard against infinity while abs(y) > tol: iter += 1 if iter > 1000: raise ValueError, "Unable to converge!" x0 -= y/fp(x0) y = f(x0) return x0
- solve x^3 - 2 = 0
return x**3 - 2
print( nrhap(f, fp, 1.0)**3 )
- Week 1: Calculator Problems
- Solutions: 
- This solution uses python list comprehensions, and shows that the slowest step is actually unit conversion!
- Week 2: Pseudocode Problems
- We will work these in class. For practice, be sure you have python list ideas down!
- Solutions: 
- Week 3: Plotting 1
- Notes on installing Matplotlib on your own computer
- Challenge: actions vs. functions -- what does the following code do: "[ print(1) for i in range(10) ]" (python3 only)
- Week 4: Plotting 2
- Solutions: 
- Week 5: Read pseudocode and plotting, focusing on types and syntax. Also, re-read the end of chapter summaries for Ch. 3, 5 and 6.
- Week 6: Read the following pieces on first order logic Chapter 1 of Jaynes' Book Hamming's Problem, and pages 71-80 in Chapter 4 of Hetland's Algorithms book (on induction). Just skip directly to chapter 4.
- Solutions 
- Week 7: Review Hetland, Chapters 3-7. We'll do a series of parsing problems and flex our collective induction and recursion muscles.
- Week 8: Minimization Problems!
- Week 10: Molecular Optimization
- Week 11: Dynamical Problems
- Week 12: Inverse Pendulum
- Week 13: No reading, simulate dynamics in the final project, scale up Rule 110 simulation and calculate histograms of local densities.
- Week 15: Elementary Statistics
- Homework will be assigned through Repl.IT
- Homework 6 is described below
- Project presentations will be on Mon.-Wed., Mar 5-7
- Final Projects are due by the end of the last week of classes.
HW6 (hand in on Wednesday, April 18)
After reading through the SIAM guide on mathematical modeling, formulate an interesting question about the chemical behavior of mixed solvent systems. Don't be too constrained by what can and can't be simulated with MD, but do create a well-designed problem statement, list some possible assumptions that could help narrow down the problem scope, and clearly explain any quantitative variables that can be defined. It will help to reference chemical literature, including textbooks from your past courses.
Homework 7 (hand in on Monday, April 23)
Create pseudocode and identify relevant python/numpy/scipy data structures and functions for carrying out each of the following MD data analysis tasks. Be explicit about data types you choose for each variable, as well as the inputs and outputs of each function you define. Your pseudocode should be returned neatly handwritten or typed in outline form and should present a complete, workable solution. 1. How do I determine the best solvent for reaction chemistry? Chemical reactions at an industrial scale are always carried out under optimized conditions. Solvent composition (i.e. the mole fraction of each solvent molecule) are one factor that can be tuned in order to control reaction rate, yield, selectivity, and ease of purification. One factor that can be modeled is the difference in solubility between reactant and product molecules. The solubility product constant for a single molecule, K, is defined as the maximum concentration the molecule can be present in solution before further additions of the molecule `crash out' and form a solid rather than dissolving. Let's test the linear solvation free energy model, K(M) = K0 exp( sum_j n_j(M) u_j ), where K0 and k_j are constants for all molecules (M), the sum runs over the 'j' functional groups present in molecule M, with n_j(M) the number of functional groups of type 'j' in that molecule. Write pseudocode for a "guess and check" method that tries out possible parameters, K0 and u_j, to decide how well those parameters perform. Assume you have a table of K(M) and n_j(M) values available, and can call the function "generate_params" to generate a new parameter guess. 2. How fast is a 'diffusion-limited' reaction? Reaction rates depend on how many times the individual molecules involved in each step of the reaction collide. For a bimolecular reaction like an enzyme-substrate binding event, E + S -> (ES), the reaction rate can be limited by diffusion when binding is very strong (so that every collision leads to "irreversible" binding). MD can't simulate chemical reactions, so the analysis will be limited to determining whether the enzyme-substrate becomes small, indicating binding. Assume you have been given a trajectory of an enzyme in solution with 100 substrate molecules. Define a (distance-based) criteria for determining if a complex has formed when given the positions of 12 atoms on a key residue of the enzyme and all 10 atoms of the substrate. This criteria should be stated in terms of a function, with well-typed inputs and outputs. Then write pseudocode that uses the entire trajectory (and the function) to determine: i) the exact time-point (frame) at which an enzyme-substrate complex forms, and ii) whether that complex ever unbinds during the remainder of the simulation trajectory. For this exercise, you will need to use steps that define and use labels for substrates and/or atoms like "loop through substrates," "label that substrate as 'k'", and "check substrate 'k'", etc. Like adding main characters to a story, it's important to know when to give certain pieces of data individual names and types.
Homework 8 - Due Friday, April 27. (hand in on Wednesday or email by Friday COB).
Revise your midterm project code using what you've learned in this class to make it more structured and maintainable. In some cases, this may simply require adding documentation and in other cases, this may be simplest to accomplish with a complete re-write. The completed code must contain: 1) at least two functions, documented with input and output types and explanations, as well as example input -> output pairs, 2) at least one loop, and 3) at least one new data structure, with a documented type and an example instance that matches the given type. One way to satisfy the last requirement is opting to change the shape of an existing array in the code or to regroup data into an array, list, or dictionary.