Difference between revisions of "CompSciWeek11"
From Predictive Chemistry
(Created page with "= Class 1 = * Tensor Manipulations - distance distributions in random point set exercise ** [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.distance.pdi…") |
|||
(3 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
+ | = Reading Assignment = |
||
+ | * [http://apps.nrbook.com/empanel/index.html Numerical Recipes] Pages 37-40 |
||
+ | * [http://calnewport.com/blog/2012/10/26/mastering-linear-algebra-in-10-days-astounding-experiments-in-ultra-learning/ Feynman Learning Technique] |
||
+ | ** Try and identify the most challenging missing concept and explain it to someone else! |
||
+ | * Optional Additional reference: Basis and Singular Value lectures from [http://codingthematrix.com/ codingthematrix.com] |
||
+ | |||
= Class 1 = |
= Class 1 = |
||
* Tensor Manipulations - distance distributions in random point set exercise |
* Tensor Manipulations - distance distributions in random point set exercise |
||
** [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.distance.pdist.html pdist] |
** [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.distance.pdist.html pdist] |
||
** [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.histogram.html histogram] |
** [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.histogram.html histogram] |
||
− | * |
+ | * Rodrigues' Rotation Formula - the preferred method when you ''have'' to work with an angle. |
* Matrix Multiplication |
* Matrix Multiplication |
||
** Basic linear algebra operations |
** Basic linear algebra operations |
||
** Timing Strassen (see also Winograd) |
** Timing Strassen (see also Winograd) |
||
+ | |||
+ | == Distribution Function Code == |
||
+ | |||
+ | This illustrates dramatic simplifications that can be obtained by using tensors in numpy. |
||
+ | |||
+ | <source lang="python"> |
||
+ | # Radial distribution function calculator for un-scaled distributions. |
||
+ | # (g(r) is get_rdf() * L**3/N and goes to 1). |
||
+ | from numpy import * |
||
+ | |||
+ | # Get the unnormalized rdf in an isotropic cubic box of side length L |
||
+ | # return unit is particles per L's length unit^3 |
||
+ | # Note that grid methods should be used for very large N. |
||
+ | def get_rdf(x, L, M): |
||
+ | D2 = x - x[:,newaxis] |
||
+ | D2 -= L*floor(D2/L+0.5) # wrap into box |
||
+ | r2 = sqrt(sum(D2*D2, -1)) |
||
+ | s = arange(M+1)*0.5*L/M |
||
+ | h, _ = histogram(r2, s) |
||
+ | h[0] -= len(x) # remove self-distances |
||
+ | # Divide by volume (times an extra (N-1), since we counted |
||
+ | # N*(N-1) distances, but the density scales as N) |
||
+ | h = h.astype(float) |
||
+ | h /= (len(x)-1) * 4*pi/3.0*(s[1:]**3 - s[:-1]**3) |
||
+ | return h |
||
+ | |||
+ | # This should give a uniform RDF of 1000 pt / L**3 (= 1 here). |
||
+ | def test(): |
||
+ | N = 1000 |
||
+ | L = 10.0 |
||
+ | x = random.random((N,3))*L |
||
+ | print get_rdf(x, L, 20) |
||
+ | </source> |
||
+ | |||
+ | Rodrigues' formula code. This is for completeness, so you can generate 3D rotations given in axis-angle notation. It is definitely not for memorization. |
||
+ | |||
+ | <source lang="python"> |
||
+ | # Build rotation matrix to rotate about an arbitrary vector using the right-hand rule. |
||
+ | # Uses Rodrigues' rotation formula (in the quaternion representation). |
||
+ | def build_rotation(u, theta): |
||
+ | n = u |
||
+ | m = sum(n*n) # check that n is normalized |
||
+ | if(m < 1.0-1.0e-10 or m > 1.0+1.0e-10): |
||
+ | n /= sqrt(m) |
||
+ | |||
+ | trans = zeros((3,3), float) |
||
+ | |||
+ | s = sin(theta) |
||
+ | c = cos(theta) |
||
+ | trans[0,0] = c + n[0]*n[0]*(1.0-c) |
||
+ | trans[0,1] = n[0]*n[1]*(1.0-c) - n[2]*s |
||
+ | trans[0,2] = n[1]*s + n[0]*n[2]*(1.0-c) |
||
+ | |||
+ | trans[1,0] = n[2]*s + n[0]*n[1]*(1.0-c) |
||
+ | trans[1,1] = c + n[1]*n[1]*(1.0-c) |
||
+ | trans[1,2] = -n[0]*s + n[1]*n[2]*(1.0-c) |
||
+ | |||
+ | trans[2,0] = -n[1]*s + n[0]*n[2]*(1.0-c) |
||
+ | trans[2,1] = n[0]*s + n[1]*n[2]*(1.0-c) |
||
+ | trans[2,2] = c + n[2]*n[2]*(1.0-c) |
||
+ | |||
+ | return trans |
||
+ | </source> |
||
= Class 2 = |
= Class 2 = |
Latest revision as of 18:12, 5 November 2014
Reading Assignment
- Numerical Recipes Pages 37-40
- Feynman Learning Technique
- Try and identify the most challenging missing concept and explain it to someone else!
- Optional Additional reference: Basis and Singular Value lectures from codingthematrix.com
Class 1
- Tensor Manipulations - distance distributions in random point set exercise
- Rodrigues' Rotation Formula - the preferred method when you have to work with an angle.
- Matrix Multiplication
- Basic linear algebra operations
- Timing Strassen (see also Winograd)
Distribution Function Code
This illustrates dramatic simplifications that can be obtained by using tensors in numpy.
<source lang="python">
- Radial distribution function calculator for un-scaled distributions.
- (g(r) is get_rdf() * L**3/N and goes to 1).
from numpy import *
- Get the unnormalized rdf in an isotropic cubic box of side length L
- return unit is particles per L's length unit^3
- Note that grid methods should be used for very large N.
def get_rdf(x, L, M):
D2 = x - x[:,newaxis] D2 -= L*floor(D2/L+0.5) # wrap into box r2 = sqrt(sum(D2*D2, -1)) s = arange(M+1)*0.5*L/M h, _ = histogram(r2, s) h[0] -= len(x) # remove self-distances # Divide by volume (times an extra (N-1), since we counted # N*(N-1) distances, but the density scales as N) h = h.astype(float) h /= (len(x)-1) * 4*pi/3.0*(s[1:]**3 - s[:-1]**3) return h
- This should give a uniform RDF of 1000 pt / L**3 (= 1 here).
def test():
N = 1000 L = 10.0 x = random.random((N,3))*L print get_rdf(x, L, 20)
</source>
Rodrigues' formula code. This is for completeness, so you can generate 3D rotations given in axis-angle notation. It is definitely not for memorization.
<source lang="python">
- Build rotation matrix to rotate about an arbitrary vector using the right-hand rule.
- Uses Rodrigues' rotation formula (in the quaternion representation).
def build_rotation(u, theta):
n = u m = sum(n*n) # check that n is normalized if(m < 1.0-1.0e-10 or m > 1.0+1.0e-10): n /= sqrt(m)
trans = zeros((3,3), float)
s = sin(theta) c = cos(theta) trans[0,0] = c + n[0]*n[0]*(1.0-c) trans[0,1] = n[0]*n[1]*(1.0-c) - n[2]*s trans[0,2] = n[1]*s + n[0]*n[2]*(1.0-c)
trans[1,0] = n[2]*s + n[0]*n[1]*(1.0-c) trans[1,1] = c + n[1]*n[1]*(1.0-c) trans[1,2] = -n[0]*s + n[1]*n[2]*(1.0-c)
trans[2,0] = -n[1]*s + n[0]*n[2]*(1.0-c) trans[2,1] = n[0]*s + n[1]*n[2]*(1.0-c) trans[2,2] = c + n[2]*n[2]*(1.0-c)
return trans
</source>
Class 2
- Implicit solutions to linear algebraic equations
- Least-squares fitting
- Repetition for point sets - project out and re-fit
- Review projection and orthogonality (see also Gram-Schmidt)
- Solution by factorization + forward / reverse substitution
- Minimization the hard way - nonlinear problems Optimize
- Transcendental problem, 5th order polynomials, etc.
- Geodesic paths via numerical optimization