Difference between revisions of "CompSciWeek13"

From Predictive Chemistry
Jump to: navigation, search
(→‎Numerical Integration + Binary Output Format)
m (→‎Python Integrator)
 
(9 intermediate revisions by the same user not shown)
Line 8: Line 8:
 
** The relevant numpy methods are x.tofile("file") and x = fromfile("file") -- but save("file", x) and x = load("file") are preferred
  
** The relevant numpy methods are x.tofile("file") and x = fromfile("file") -- but save("file", x) and x = load("file") are preferred
 
** inspecting binary formats with od
  
** inspecting binary formats with od
  +
* [[Media:Integrate.pdf|Lecture Slides]]
   
  +
=== GNU ODE Integrator ===
  +
  +
'''vanderpol.ode'''
  +
  +
<source lang="haskell">
  +
mu = %MU
  +
  +
u = 1.0
  +
y = 0.0
  +
  +
u' = y
  +
y' = mu*(1-u*u)*y - u
  +
  +
print t, u, y, y! every 10 from 50
  +
step 0, 100
  +
</source>
  +
  +
You can do the substitution with sed and run all dynamics with:
  +
  +
for((i=0;i<10;i++)); do
  +
sed -e "s|%MU|$i|" vanderpol.ode >vanderpol-$i.ode
  +
ode -f vanderpol-$i.ode >out-$i.dat </dev/null
  +
done
  +
  +
=== Plotting with Gnuplot ===
  +
  +
<source lang="gnuplot">
  +
set term postscript enh color lw 3 "Times-Roman" 28
  +
set out 'van.eps'
  +
  +
set xlabel "u"
  +
set ylabel "y"
  +
plot 'out-0.dat' u 2:3 w l title "0", \
  +
'out-1.dat' u 2:3 w l title "1", \
  +
'out-4.dat' u 2:3 w l title "{/Symbol m} = 4"
  +
</source>
  +
  +
Run with:
  +
  +
gnuplot plot_van.gplot
  +
  +
=== Python Integrator ===
  +
  +
First, form the heart of the algorithm:
  +
<source lang="python">
  +
def step(N=1):
  +
for i in xrange(N):
  +
v += (mu*(1.0-x**2)*v - x) * dt
  +
x += v * dt
  +
</source>
  +
Note that this will work if x and v are vectors as well as scalars.
  +
That means we can run a bunch of trajectories with different initial starting points
  +
in one fell swoop.
  +
  +
Next, wrap it inside a class:
  +
<source lang="python">
  +
class State:
  +
def __init__(self, x, v, mu, dt):
  +
assert x.shape == v.shape
  +
self.mu = mu
  +
self.dt = dt
  +
self.x = x
  +
self.v = v
  +
def state(self, N=1):
  +
# must change all variables (e.g. "name") to "self.name"
  +
</source>
  +
  +
Next, instantiate and run the class:
  +
  +
<source lang="python">
  +
st = State(5.0, 1.0, 2.0, 0.01)
  +
st.step()
  +
print st.x
  +
st.step()
  +
print st.x
  +
</source>
  +
  +
Next, add a nice printing function to "print the state,"
  +
<source lang="python">
  +
def __str__(self):
  +
s = ""
  +
for x,v in zip(self.x, self.v):
  +
s += " %f %f"%(x, v)
  +
return s
  +
</source>
  +
This prints (x1 v1 x2 v2 ...) in a single, long line.
  +
  +
OK, so now we're ready to run lots of tests. Let's
  +
lift the python program to a command-line program.
  +
Add this to the header:
  +
<source lang="python">
  +
#!/usr/bin/env python
  +
  +
import sys
  +
from numpy import *
  +
from numpy.random import standard_normal
  +
norm = standard_normal
  +
</source>
  +
  +
Add this to the end of the file:
  +
<source lang="python">
  +
def main(argv):
  +
print argv
  +
  +
if __name__=="__main__":
  +
main(sys.argv)
  +
</source>
  +
  +
From the shell, run
  +
<source lang="bash">
  +
chmod +x vanderpol.py
  +
./vanderpol.py a test argument
  +
</source>
  +
  +
Note that you get to provide parameters to your program at the command-line.
  +
Unfortunately, they all have type str. Let's formalize them,
  +
<source lang="python">
  +
def main(argv):
  +
assert len(argv) == 5, "Usage: %s x0 p0 R mu"%(argv[0])
  +
x0, p0, R, mu = map(float, argv[1:5])
  +
st = State(x0,p0,mu,0.01)
  +
for i in range(500):
  +
st.step(10)
  +
print st
  +
</source>
  +
  +
Finally, we have a new tool that can be run to generate reproducible simulations
  +
of the van der Pol oscillator. Next, we need to consider how to store
  +
the data it generates for later analysis and use.
  +
  +
Here's an example using a 500x2x5 tensor:
  +
<source lang="python">
  +
T = zeros((500,2,5)) # steps by (x or v) by points
  +
for i in range(500):
  +
st.step(10)
  +
T[i,0] = st.x
  +
T[i,1] = st.v
  +
</source>
  +
  +
This can be "flattened" to a 500x10 tensor,
  +
<source lang="python">
  +
T = reshape(T, (500,10)) # or reshape(T, (-1,10)) "figures out" the 500 size
  +
</source>
  +
  +
=== Working with binary ===
 
Example: [http://man7.org/linux/man-pages/man5/elf.5.html ELF] header
  
Example: [http://man7.org/linux/man-pages/man5/elf.5.html ELF] header
   
Line 27: Line 173:
 
|-
  
|-
 
! 2 bytes
  
! 2 bytes
−
! = 1 short
 +
! = 1 short = 1 word
 
|-
  
|-
 
! 4 bytes
  
! 4 bytes
Line 33: Line 179:
 
|-
  
|-
 
! 8 bytes
  
! 8 bytes
−
! = 1 64-bit int = 1 64-bit address = 1 double = 1 float complex
 +
! = 1 64-bit int = 1 double = 1 float complex
 
|-
  
|-
 
! 16 bytes
  
! 16 bytes
Line 39: Line 185:
 
|}
  
|}
   
−
Byte ordering on most modern 64-bit processors is little-endian (Intel, AMD)
 +
Byte ordering on most modern 64-bit processors is little-endian (Intel, AMD). For more info, see [http://betterexplained.com/articles/understanding-big-and-little-endian-byte-order/].
   
−
in base 10, this would mean the 4-digit representation of the number 123 is
 +
In an imaginary system where base 10 numbers are encoded by 4-bit nibbles, a 4-digit representation of the number 123 is
   
 
3 = 0011
  
3 = 0011
Line 52: Line 198:
 
0011 0010 0001 0000
  
0011 0010 0001 0000
   
−
Since we have to use 4 digits, but the little end (least significant byte) goes first.
  
  +
Since we have to use 4 digits, the most significant digit is set to all zeros, and goes last. The little end (least significant byte) goes first. Notice that within the 4-bit nibbles, the least-significant bit is still on the right. So the number is encoded using this left-to right (within a digit) and bottom to top (digits within a word) order.
  +
  +
Real binary encoding is just like the above, except that each 'digit' is a byte, so every 8 bits are in order, but the bytes are reversed.
   
 
= Basis Functions =
  
= Basis Functions =
Line 58: Line 206:
 
* Representing the differentiation operator
  
* Representing the differentiation operator
 
* Solving PDEs using the implicit method
  
* Solving PDEs using the implicit method
  +
* [[Media:Basis.pdf|Lecture Slides]]

Latest revision as of 12:00, 6 April 2016

Numerical Integration + Binary Output Format

  • Van der Pol Oscillator: <math>\ddot u - \mu (1-u^2) \dot u + u = 0</math>
    • Implementations in OpenOffice Spreadsheet, GNU ODE, and Python
    • Perturbed initial conditions, Lyapunov exponents, and large deviation principles
  • Working well with others: text
  • Working with machines: binary
    • This is preferable when you have lots and lots of data
    • The relevant numpy methods are x.tofile("file") and x = fromfile("file") -- but save("file", x) and x = load("file") are preferred
    • inspecting binary formats with od
  • Lecture Slides

GNU ODE Integrator

vanderpol.ode

<source lang="haskell"> mu = %MU

u = 1.0 y = 0.0

u' = y y' = mu*(1-u*u)*y - u

print t, u, y, y! every 10 from 50 step 0, 100 </source>

You can do the substitution with sed and run all dynamics with: 

 for((i=0;i<10;i++)); do
   sed -e "s|%MU|$i|" vanderpol.ode >vanderpol-$i.ode
   ode -f vanderpol-$i.ode >out-$i.dat </dev/null
 done

Plotting with Gnuplot

<source lang="gnuplot"> set term postscript enh color lw 3 "Times-Roman" 28 set out 'van.eps'

set xlabel "u" set ylabel "y" plot 'out-0.dat' u 2:3 w l title "0", \ 

    'out-1.dat' u 2:3 w l title "1", \
    'out-4.dat' u 2:3 w l title "{/Symbol m} = 4"

</source>

Run with: 

gnuplot plot_van.gplot

Python Integrator

First, form the heart of the algorithm: <source lang="python"> def step(N=1): 

   for i in xrange(N):
       v += (mu*(1.0-x**2)*v - x) * dt
       x += v * dt

</source> Note that this will work if x and v are vectors as well as scalars. That means we can run a bunch of trajectories with different initial starting points in one fell swoop.

Next, wrap it inside a class: <source lang="python"> class State: 

   def __init__(self, x, v, mu, dt):
       assert x.shape == v.shape
       self.mu = mu
       self.dt = dt
       self.x = x
       self.v = v
   def state(self, N=1):
  1. must change all variables (e.g. "name") to "self.name"

</source>

Next, instantiate and run the class:

<source lang="python"> st = State(5.0, 1.0, 2.0, 0.01) st.step() print st.x st.step() print st.x </source>

Next, add a nice printing function to "print the state," <source lang="python"> 

   def __str__(self):
       s = ""
       for x,v in zip(self.x, self.v):
           s += " %f %f"%(x, v)
       return s

</source> This prints (x1 v1 x2 v2 ...) in a single, long line.

OK, so now we're ready to run lots of tests. Let's lift the python program to a command-line program. Add this to the header: <source lang="python">

  1. !/usr/bin/env python

import sys from numpy import * from numpy.random import standard_normal norm = standard_normal </source>

Add this to the end of the file: <source lang="python"> def main(argv): 

   print argv

if __name__=="__main__": 

   main(sys.argv)

</source>

From the shell, run <source lang="bash"> chmod +x vanderpol.py ./vanderpol.py a test argument </source>

Note that you get to provide parameters to your program at the command-line. Unfortunately, they all have type str. Let's formalize them, <source lang="python"> def main(argv): 

   assert len(argv) == 5, "Usage: %s x0 p0 R mu"%(argv[0])
   x0, p0, R, mu = map(float, argv[1:5])
   st = State(x0,p0,mu,0.01)
   for i in range(500):
       st.step(10)
       print st

</source>

Finally, we have a new tool that can be run to generate reproducible simulations of the van der Pol oscillator. Next, we need to consider how to store the data it generates for later analysis and use.

Here's an example using a 500x2x5 tensor: <source lang="python"> T = zeros((500,2,5)) # steps by (x or v) by points for i in range(500): 

   st.step(10)
   T[i,0] = st.x
   T[i,1] = st.v

</source>

This can be "flattened" to a 500x10 tensor, <source lang="python"> T = reshape(T, (500,10)) # or reshape(T, (-1,10)) "figures out" the 500 size </source>

Working with binary

Example: ELF header 

od -t x1 -t c -N 8 /bin/bash

Binary comes in lots of units


8 bits = 1 byte
4 bits = 1 nibble
1 byte = 1 ascii character
2 bytes = 1 short = 1 word
4 bytes = 1 32-bit int = 1 float
8 bytes = 1 64-bit int = 1 double = 1 float complex
16 bytes = 1 long double = 1 double complex

Byte ordering on most modern 64-bit processors is little-endian (Intel, AMD). For more info, see [1].

In an imaginary system where base 10 numbers are encoded by 4-bit nibbles, a 4-digit representation of the number 123 is 

3 = 0011
2 = 0010
1 = 0001
0 = 0000

The binary together would read: 

0011 0010 0001 0000

Since we have to use 4 digits, the most significant digit is set to all zeros, and goes last. The little end (least significant byte) goes first. Notice that within the 4-bit nibbles, the least-significant bit is still on the right. So the number is encoded using this left-to right (within a digit) and bottom to top (digits within a word) order.

Real binary encoding is just like the above, except that each 'digit' is a byte, so every 8 bits are in order, but the bytes are reversed.

Basis Functions

  • Constructing B-splines - tensor method
  • Representing the differentiation operator
  • Solving PDEs using the implicit method
  • Lecture Slides
Retrieved from "https://www.predictivestatmech.org/w/index.php?title=CompSciWeek13&oldid=504"