Skip to content

WIP: chempy refactoring #433

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

WIP: chempy refactoring #433

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f92a13f
adding tests for chempy.cpv
ye11owSub Feb 15, 2025
1a02f19
adding tests for chempy models
ye11owSub Feb 20, 2025
6ec0804
refactoring chempy/models.py
ye11owSub Apr 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 58 additions & 68 deletions modules/chempy/cpv.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,43 +21,46 @@
#
# TODO: documentation!

from typing import Sequence
import math
import random
import copy

RSMALL4 = 0.0001

#------------------------------------------------------------------------------
def get_null():
def get_null() -> list[float]:
return [0.0,0.0,0.0]

#------------------------------------------------------------------------------
def get_identity():
def get_identity() -> list[list[float]]:
return [[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]

#------------------------------------------------------------------------------
def distance_sq(v1, v2):
def distance_sq(v1: Sequence[float], v2: Sequence[float]) -> float:
d0 = v2[0] - v1[0]
d1 = v2[1] - v1[1]
d2 = v2[2] - v1[2]
return (d0*d0) + (d1*d1) + (d2*d2)

#------------------------------------------------------------------------------
def distance(v1, v2):
def distance(v1: Sequence[float], v2: Sequence[float]) -> float:
d0 = v2[0] - v1[0]
d1 = v2[1] - v1[1]
d2 = v2[2] - v1[2]
return math.sqrt((d0*d0) + (d1*d1) + (d2*d2))

#------------------------------------------------------------------------------
def length(v):
def length(v: Sequence[float]) -> float:
return math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])

#------------------------------------------------------------------------------
def random_displacement(v,radius):
r_vect = lambda r=random.random:[r()-0.5,r()-0.5,r()-0.5]
def _random_vector() -> list[float]:
r = random.random
return [r()-0.5,r()-0.5,r()-0.5]

def random_displacement(v: Sequence[float], radius: float) -> list[float]:
while 1:
vect = r_vect()
vect = _random_vector()
v_len = length(vect)
if (v_len<=0.5):
break;
Expand All @@ -68,92 +71,85 @@ def random_displacement(v,radius):
return v

#------------------------------------------------------------------------------
def random_sphere(v,radius):
r_vect = lambda r=random.random:[r()-0.5,r()-0.5,r()-0.5]
def random_sphere(v: Sequence[float], radius: float) -> list[float]:
while 1:
vect = r_vect()
vect = _random_vector()
v_len = length(vect)
if (v_len<=0.5) and (v_len!=0.0):
break;
return add(v,scale([vect[0], vect[1], vect[2]],2*radius/v_len))

#------------------------------------------------------------------------------
def random_vector():
r_vect = lambda r=random.random:[r()-0.5,r()-0.5,r()-0.5]
def random_vector() -> list[float]:
while 1:
vect = r_vect()
vect = _random_vector()
if length(vect)<=0.5:
break;
return scale([vect[0], vect[1], vect[2]],2.0)

#------------------------------------------------------------------------------
def add(v1,v2):
def add(v1: Sequence[float], v2: Sequence[float]) -> list[float]:
return [v1[0]+v2[0],v1[1]+v2[1],v1[2]+v2[2]]

#------------------------------------------------------------------------------
def average(v1,v2):
def average(v1: Sequence[float], v2: Sequence[float]) -> list[float]:
return [(v1[0]+v2[0])/2.0,(v1[1]+v2[1])/2.0,(v1[2]+v2[2])/2.0]

#------------------------------------------------------------------------------
def scale(v,factor):
def scale(v: Sequence[float], factor: float) -> list[float]:
return [v[0]*factor,v[1]*factor,v[2]*factor]

#------------------------------------------------------------------------------
def negate(v):
def negate(v: Sequence[float]) -> list[float]:
return [-v[0],-v[1],-v[2]]

#------------------------------------------------------------------------------
def sub(v1,v2):
def reverse(v: Sequence[float]) -> list[float]:
return [ -v[0], -v[1], -v[2] ]

#------------------------------------------------------------------------------
def sub(v1: Sequence[float], v2: Sequence[float]) -> list[float]:
return [v1[0]-v2[0],v1[1]-v2[1],v1[2]-v2[2]]

#------------------------------------------------------------------------------
def dot_product(v1,v2):
def dot_product(v1: Sequence[float], v2: Sequence[float]) -> float:
return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]

#------------------------------------------------------------------------------
def cross_product(v1,v2):
def cross_product(v1: Sequence[float], v2: Sequence[float]) -> list[float]:
return [(v1[1]*v2[2]) - (v1[2]*v2[1]),
(v1[2]*v2[0]) - (v1[0]*v2[2]),
(v1[0]*v2[1]) - (v1[1]*v2[0])]

#------------------------------------------------------------------------------
def transform(m,v):
def transform(m: Sequence[Sequence[float]], v: Sequence[float]) -> list[float]:
return [m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2],
m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2],
m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2]]

#------------------------------------------------------------------------------
def inverse_transform(m,v):
def inverse_transform(m: Sequence[Sequence[float]], v: Sequence[float]) -> list[float]:
return [m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2],
m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2],
m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2]]

#------------------------------------------------------------------------------
def multiply(m1,m2): # HAVEN'T YET VERIFIED THAT THIS CONFORMS TO STANDARD DEFT
return [[m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0],
m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0],
m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0]],
[m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1],
m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1],
m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1]],
[m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2],
m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2],
m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2]]]
def multiply(m1: Sequence[Sequence[float]], m2: Sequence[Sequence[float]]) -> list[list[float]]:
# HAVEN'T YET VERIFIED THAT THIS CONFORMS TO STANDARD DEFT
# upd: no, it's not(fixed)
Comment on lines +138 to +139
Copy link
Member

@JarrettSJohnson JarrettSJohnson May 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be removed.


#------------------------------------------------------------------------------
def transpose(m1):
return [[m1[0][0],
m1[1][0],
m1[2][0]],
[m1[0][1],
m1[1][1],
m1[2][1]],
[m1[0][2],
m1[1][2],
m1[2][2]]]
return [[m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0],
m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1] + m1[0][2]*m2[2][1],
m1[0][0]*m2[0][2] + m1[0][1]*m2[1][2] + m1[0][2]*m2[2][2]],
[m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0],
m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1] + m1[1][2]*m2[2][1],
m1[1][0]*m2[0][2] + m1[1][1]*m2[1][2] + m1[1][2]*m2[2][2]],
[m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0],
m1[2][0]*m2[0][1] + m1[2][1]*m2[1][1] + m1[2][2]*m2[2][1],
m1[2][0]*m2[0][2] + m1[2][1]*m2[1][2] + m1[2][2]*m2[2][2]]]
Copy link
Contributor Author

@ye11owSub ye11owSub Feb 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Matrix multiplication was not implemented correctly. This and code duplication has also been fixed in this PR


#------------------------------------------------------------------------------
def get_system2(x,y):
def get_system2(x: Sequence[float],y: Sequence[float]) -> list[list[float]]:
z = cross_product(x,y)
z = normalize(z)
y = cross_product(z,x);
Expand All @@ -162,24 +158,24 @@ def get_system2(x,y):
return [x,y,z]

#------------------------------------------------------------------------------
def scale_system(s,factor):
def scale_system(s: Sequence[Sequence[float]], factor: float) -> list[list[float]]:
r = []
for a in s:
r.append([a[0]*factor,a[1]*factor,a[2]*factor])
return r

#------------------------------------------------------------------------------
def transpose(m):
def transpose(m: Sequence[Sequence[float]]) -> list[list[float]]:
return [[m[0][0], m[1][0], m[2][0]],
[m[0][1], m[1][1], m[2][1]],
[m[0][2], m[1][2], m[2][2]]]

#------------------------------------------------------------------------------
def transform_about_point(m,v,p):
def transform_about_point(m: Sequence[Sequence[float]], v: Sequence[float], p: Sequence[float]) -> list[float]:
return add(transform(m,sub(v,p)),p)

#------------------------------------------------------------------------------
def get_angle(v1,v2): # v1,v2 must be unit vectors
def get_angle(v1: Sequence[float], v2: Sequence[float]) -> float: # v1,v2 must be unit vectors
denom = (math.sqrt(((v1[0]*v1[0]) + (v1[1]*v1[1]) + (v1[2]*v1[2]))) *
math.sqrt(((v2[0]*v2[0]) + (v2[1]*v2[1]) + (v2[2]*v2[2]))))
if denom>1e-10:
Expand All @@ -190,7 +186,7 @@ def get_angle(v1,v2): # v1,v2 must be unit vectors
return result

#------------------------------------------------------------------------------
def get_angle_formed_by(p1,p2,p3): # angle formed by three positions in space
def get_angle_formed_by(p1: Sequence[float], p2: Sequence[float], p3: Sequence[float]) -> float: # angle formed by three positions in space

# based on code submitted by Paul Sherwood
r1 = distance(p1,p2)
Expand All @@ -207,37 +203,33 @@ def get_angle_formed_by(p1,p2,p3): # angle formed by three positions in space
return theta;

#------------------------------------------------------------------------------
def project(v,n):
def project(v: Sequence[float], n: Sequence[float]) -> list[float]:
dot = v[0]*n[0] + v[1]*n[1] + v[2]*n[2]
return [ dot * n[0], dot * n[1], dot * n[2] ]

#------------------------------------------------------------------------------
def remove_component(v, n):
def remove_component(v: Sequence[float], n: Sequence[float]) -> list[float]:
dot = v[0]*n[0] + v[1]*n[1] + v[2]*n[2]
return [v[0] - dot * n[0], v[1] - dot * n[1], v[2] - dot * n[2]]

#------------------------------------------------------------------------------
def normalize(v):
def normalize(v: Sequence[float]) -> list[float]:
vlen = math.sqrt((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]))
if vlen>RSMALL4:
return [v[0]/vlen,v[1]/vlen,v[2]/vlen]
else:
return get_null()

#------------------------------------------------------------------------------
def reverse(v):
return [ -v[0], -v[1], -v[2] ]

#------------------------------------------------------------------------------
def normalize_failsafe(v):
def normalize_failsafe(v: Sequence[float]) -> list[float]:
vlen = math.sqrt((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]))
if vlen>RSMALL4:
return [v[0]/vlen,v[1]/vlen,v[2]/vlen]
else:
return [1.0,0.0,0.0]

#------------------------------------------------------------------------------
def rotation_matrix(angle,axis):
def rotation_matrix(angle: float, axis: Sequence[float]) -> list[list[float]]:

x=axis[0]
y=axis[1]
Expand Down Expand Up @@ -271,7 +263,7 @@ def rotation_matrix(angle,axis):
[ (one_c * zx) - ys, (one_c * yz) + xs, (one_c * zz) + c ]]

#------------------------------------------------------------------------------
def transform_array(rot_mtx,vec_array):
def transform_array(rot_mtx: Sequence[Sequence[float]], vec_array: Sequence[Sequence[float]]) -> list[list[float]]:

'''transform_array( matrix, vector_array ) -> vector_array

Expand All @@ -280,7 +272,7 @@ def transform_array(rot_mtx,vec_array):
return [transform(rot_mtx, x) for x in vec_array]

#------------------------------------------------------------------------------
def translate_array(trans_vec,vec_array):
def translate_array(trans_vec: Sequence[float], vec_array: Sequence[Sequence[float]]) -> list[list[float]]:

'''translate_array(trans_vec,vec_array) -> vec_array

Expand All @@ -291,17 +283,17 @@ def translate_array(trans_vec,vec_array):
return [add(trans_vec, x) for x in vec_array]

#------------------------------------------------------------------------------
def fit_apply(fit_result,vec_array):
def fit_apply(fit_result: tuple[Sequence[float], Sequence[float], Sequence[Sequence[float]], float], vec_array: Sequence[Sequence[float]]) -> list[list[float]]:
'''fit_apply(fir_result,vec_array) -> vec_array

Applies a fit result to an array of vectors
'''

t1, mt2, m = fit_result[:3]
return [add(t1, transform(m, add(mt2, x))) for x in vec_array]

#------------------------------------------------------------------------------
def fit(target_array, source_array):
def fit(target_array: Sequence[Sequence[float]], source_array: Sequence[Sequence[float]]) -> tuple[list[float], list[float], list[list[float]], float]:

'''fit(target_array, source_array) -> (t1, t2, rot_mtx, rmsd) [fit_result]

Expand Down Expand Up @@ -370,7 +362,6 @@ def fit(target_array, source_array):
iters = 0
while (iters < maxiter):
iters = iters + 1
ix = (iters-1)%ndim
iy = iters%ndim
iz = (iters+1)%ndim
sig = corr_mtx[iz][iy] - corr_mtx[iy][iz]
Expand Down Expand Up @@ -403,5 +394,4 @@ def fit(target_array, source_array):
return(t1, t2, rot_mtx, rmsd)

# Too many iterations; something wrong.
print ("Error: Too many iterations in RMS fit.")
raise ValueError
raise ValueError("Error: Too many iterations in RMS fit.")
Loading
Loading