Python in Scientific Computing: Why and How?

Sage Shaw - Oct 26^th, 2023

Graduate Student Seminar

Who are you?

Would you rather...

Who are you?

Do you use Python? v3.11

Do you work with sequences? (yes you do)

Are your codes proceedural?

Functional?

Object oriented?

Who are you?

Have you written a class?

An __iter__ method?

A generator comprehension?

Who are you?

How many keywords do you recognize (use)?

False      await      else       import     pass
None       break      except     in         raise
True       class      finally    is         return
and        continue   for        lambda     try
as         def        from       nonlocal   while
assert     del        global     not        with
async      elif       if         or         yield
match      case

Disclaimers

• Working code is good.
• Python is slow^*.
• This talk is about generators.

Loops

Types of loops

• for loop
• while loop
• do while loop
• foreach loop
• recursion?

Loops are a kind of control flow.

Often defined as a primitive recursive function.

What kind of loops are these?

n = 0
while not done(n):
    print(n)
    n += 1

from itertools import count as nonnegative_integers
for n in nonnegative_integers():
    if done(n): break
    print(n)

Python "for loops"

 # What you write
for n in range(3):
    #  loop body
    print(n)

 # What happens behind the scenes
generator_object = range(3).__iter__()
while True:
    try:
        n = generator_object.__next__()
    except StopIteration:
        break
    #  loop body
    print(n)

Why does Python work this way?

Python loops are flexible

for x, y in zip(x_values, y_values): 

for alpha, beta in product(alphas, betas):

for time, u, v, w in solver.solution_generator():

for index, state in enumerate(markov_chain):

Generators

Look like functions.

def example():
	yield 1
	yield 2
	yield 3

But returns generator objects.

def example():
	yield 1
	yield 2
	yield 3

gen = example()
type(gen)  # class 'generator'
next(gen)  # 1
next(gen)  # 2
next(gen)  # 3
next(gen)  # StopIteration exception

The Single Responsibility Principle

A unit of code should be responsible for one thing.

(First in the SOLID principles.)

Example: Rootfinding

Newton's Method

Newton's Method

How many things is this function doing?

def newton_old(x0, f, df):
    seq = [x0]
    x = x0
    for _ in range(10):
        x = x - f(x) / df(x)
        seq.append(x)
    return seq

Let's separate the pieces.

The update step.

def newton_update(x, f, df):
    return x - f(x)/df(x)
		

Let's separate the pieces.

The iteration.

def iterate(x, func):
    while True:
        yield x
        x = func(x)

Let's separate the pieces.

Truncate to length 10 and store in a list.

def take10(seq):
    return list(islice(seq, 0, 10))
        

And combine them.

Make the infinte Newton sequence.

def newton_seq(x, f, df):
    update = partial(newton_update, f=f, df=df)
    yield from iterate(x, update)
        

And combine them.

Truncate and make into a list.

def newton_same_as_old(x0, f, df):
    return take10(newton_seq(x0, f, df))
        

The Open/Closed principle.

Code should be open for extension, but closed to modification.

(Second in the SOLID principles.)

How would you change the old code to do...

Multidimesional Newton's?

Broyden's method?

To use an $f$ tolarance?

To use a Cauchy tolarance?

To have a maximum length?

1D secant method?

Questions

For you...

• Did I convince you generators are useful?
• Are you going to explore other Python features?
• Do you want to help me build a Python community?

Questions for me?

Slides: shawsa.github.io