Skip to content
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

Jump to bottom

PEP8 in sims #884 #960

Closed
wants to merge 2 commits into from
Closed

PEP8 in sims #884 #960

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5e6db15
PEP8 in sims #884 Fixes https://github.com/nipy/dipy/issues/884
theaverageguy Mar 2, 2016
aab3b72
Fixed errors
theaverageguy Mar 2, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion dipy/sims/__init__.py
View file
Open in desktop
@@ -1 +1 @@
#init for simulations
# init for simulations
76 changes: 38 additions & 38 deletions dipy/sims/phantom.py
View file
Open in desktop
Expand Up @@ -62,12 +62,12 @@ def add_noise(vol, snr=1.0, S0=None, noise_type='rician'):
return np.reshape(vol_flat, orig_shape)


def diff2eigenvectors(dx,dy,dz):
def diff2eigenvectors(dx, dy, dz):
""" numerical derivatives 2 eigenvectors
"""
basis = np.eye(3)
u = np.array([dx, dy, dz])
u = u/np.linalg.norm(u)
u = u / np.linalg.norm(u)
R = vec2vec_rotmat(basis[:, 0], u)
eig0 = u
eig1 = np.dot(R, basis[:, 1])
Expand Down Expand Up @@ -197,54 +197,54 @@ def orbital_phantom(gtab=None,

if __name__ == "__main__":

## TODO: this can become a nice tutorial for generating phantoms
# TODO: this can become a nice tutorial for generating phantoms

def f(t):
x=np.sin(t)
y=np.cos(t)
#z=np.zeros(t.shape)
z=np.linspace(-1,1,len(x))
return x,y,z
x = np.sin(t)
y = np.cos(t)
# z=np.zeros(t.shape)
z = np.linspace(-1, 1, len(x))
return x, y, z

#helix
vol=orbital_phantom(func=f)
# helix
vol = orbital_phantom(func=f)

def f2(t):
x=np.linspace(-1,1,len(t))
y=np.linspace(-1,1,len(t))
z=np.zeros(x.shape)
return x,y,z
x = np.linspace(-1, 1, len(t))
y = np.linspace(-1, 1, len(t))
z = np.zeros(x.shape)
return x, y, z

#first direction
vol2=orbital_phantom(func=f2)
# first direction
vol2 = orbital_phantom(func=f2)

def f3(t):
x=np.linspace(-1,1,len(t))
y=-np.linspace(-1,1,len(t))
z=np.zeros(x.shape)
return x,y,z
x = np.linspace(-1, 1, len(t))
y = -np.linspace(-1, 1, len(t))
z = np.zeros(x.shape)
return x, y, z

#second direction
vol3=orbital_phantom(func=f3)
#double crossing
vol23=vol2+vol3
# second direction
vol3 = orbital_phantom(func=f3)
# double crossing
vol23 = vol2 + vol3

#"""
# """
def f4(t):
x=np.zeros(t.shape)
y=np.zeros(t.shape)
z=np.linspace(-1,1,len(t))
return x,y,z
x = np.zeros(t.shape)
y = np.zeros(t.shape)
z = np.linspace(-1, 1, len(t))
return x, y, z

#triple crossing
vol4=orbital_phantom(func=f4)
vol234=vol23+vol4
# triple crossing
vol4 = orbital_phantom(func=f4)
vol234 = vol23 + vol4

voln=add_rician_noise(vol234)
voln = add_rician_noise(vol234)

#"""
# """

#r=fvtk.ren()
#fvtk.add(r,fvtk.volume(vol234[...,0]))
#fvtk.show(r)
#vol234n=add_rician_noise(vol234,20)
# r=fvtk.ren()
# fvtk.add(r,fvtk.volume(vol234[...,0]))
# fvtk.show(r)
# vol234n=add_rician_noise(vol234,20)
Copy link
Contributor

Choose a reason for hiding this comment

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

Somebody might want to clean this up at some point, but maybe not on this PR.

2 changes: 1 addition & 1 deletion dipy/sims/tests/__init__.py
View file
Open in desktop
@@ -1,4 +1,4 @@
# Test callable
from numpy.testing import Tester
test = Tester().test
del Tester
del Tester
46 changes: 23 additions & 23 deletions dipy/sims/tests/test_phantom.py
View file
Open in desktop
Expand Up @@ -17,10 +17,10 @@
from dipy.core.gradients import gradient_table


fimg,fbvals,fbvecs=get_data('small_64D')
bvals=np.load(fbvals)
bvecs=np.load(fbvecs)
bvecs[np.isnan(bvecs)]=0
fimg, fbvals, fbvecs = get_data('small_64D')
bvals = np.load(fbvals)
bvecs = np.load(fbvecs)
bvecs[np.isnan(bvecs)] = 0

gtab = gradient_table(bvals, bvecs)

Expand All @@ -29,10 +29,10 @@ def f(t):
"""
Helper function used to define a mapping time => xyz
"""
x = np.linspace(-1,1,len(t))
y = np.linspace(-1,1,len(t))
z = np.linspace(-1,1,len(t))
return x,y,z
x = np.linspace(-1, 1, len(t))
y = np.linspace(-1, 1, len(t))
z = np.linspace(-1, 1, len(t))
return x, y, z


def test_phantom():
Expand All @@ -55,32 +55,32 @@ def test_phantom():
FA[np.isnan(FA)] = 0
# 686 -> expected FA given diffusivities of [1500, 400, 400]
l1, l2, l3 = 1500e-6, 400e-6, 400e-6
expected_fa = (np.sqrt(0.5) * np.sqrt((l1 - l2)**2 + (l2-l3)**2 + (l3-l1)**2 )/np.sqrt(l1**2 + l2**2 + l3**2))
expected_fa = np.sqrt(0.5) * np.sqrt((l1 - l2)**2 + (l2 - l3)**2 + (l3 - l1)**2) / np.sqrt(l1**2 + l2**2 + l3**2))
Copy link
Contributor

Choose a reason for hiding this comment

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

This line is probably still >80 characters long. A suggestion:

    expected_fa = (np.sqrt(0.5) * np.sqrt((l1 - l2)**2 + 
                   (l2 - l3)**2 + (l3 - l1)**2) / 
                   np.sqrt(l1**2 + l2**2 + l3**2)))


assert_array_almost_equal(FA.max(), expected_fa, decimal=2)
assert_array_almost_equal(FA.max(), expected_fa, decimal = 2)
Copy link
Contributor

Choose a reason for hiding this comment

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

No space needed here.



def test_add_noise():
np.random.seed(1980)

N = 50
S0 = 100
N=50
Copy link
Contributor

Choose a reason for hiding this comment

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

But space would be nice here.

S0=100

options = dict(func=f,
t=np.linspace(0, 2 * np.pi, N),
datashape=(10, 10, 10, len(bvals)),
origin=(5, 5, 5),
scale=(3, 3, 3),
angles=np.linspace(0, 2 * np.pi, 16),
radii=np.linspace(0.2, 2, 6),
S0=S0)
options=dict(func = f,
t = np.linspace(0, 2 * np.pi, N),
datashape = (10, 10, 10, len(bvals)),
origin = (5, 5, 5),
scale = (3, 3, 3),
angles = np.linspace(0, 2 * np.pi, 16),
radii = np.linspace(0.2, 2, 6),
S0 = S0)
Copy link
Contributor

Choose a reason for hiding this comment

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

No spaces around "=" sign when assigning to a dict or when passing arguments to a function, please.


vol = orbital_phantom(gtab, **options)
vol=orbital_phantom(gtab, **options)
Copy link
Contributor

Choose a reason for hiding this comment

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

But space around "=" when used in assignment!


for snr in [10, 20, 30, 50]:
vol_noise = orbital_phantom(gtab, snr=snr, **options)
vol_noise=orbital_phantom(gtab, snr = snr, **options)
Copy link
Contributor

Choose a reason for hiding this comment

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

Same as above.


sigma = S0 / snr
sigma=S0 / snr
Copy link
Contributor

Choose a reason for hiding this comment

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

And same here.


assert_(np.abs(np.var(vol_noise - vol) - sigma ** 2) < 1)

Expand Down
30 changes: 15 additions & 15 deletions dipy/sims/tests/test_voxel.py
View file
Open in desktop
Expand Up @@ -59,7 +59,7 @@ def test_check_directions():
# Testing more than one direction simultaneously
angles = np.array([[90, 0], [30, 0]])
sticks = _check_directions(angles)
ref_vec = [np.sin(np.pi*30/180), 0, np.cos(np.pi*30/180)]
ref_vec = [np.sin(np.pi * 30 / 180), 0, np.cos(np.pi * 30 / 180)]
assert_array_almost_equal(sticks, [[1, 0, 0], ref_vec])
# Testing directions not aligned to planes x = 0, y = 0, or z = 0
the1 = 0
Expand All @@ -68,12 +68,12 @@ def test_check_directions():
phi2 = 45
angles = np.array([(the1, phi1), (the2, phi2)])
sticks = _check_directions(angles)
ref_vec1 = (np.sin(np.pi*the1/180) * np.cos(np.pi*phi1/180),
np.sin(np.pi*the1/180) * np.sin(np.pi*phi1/180),
np.cos(np.pi*the1/180))
ref_vec2 = (np.sin(np.pi*the2/180) * np.cos(np.pi*phi2/180),
np.sin(np.pi*the2/180) * np.sin(np.pi*phi2/180),
np.cos(np.pi*the2/180))
ref_vec1 = (np.sin(np.pi * the1 / 180) * np.cos(np.pi * phi1 / 180),
np.sin(np.pi * the1 / 180) * np.sin(np.pi * phi1 / 180),
np.cos(np.pi * the1 / 180))
ref_vec2 = (np.sin(np.pi * the2 / 180) * np.cos(np.pi * phi2 / 180),
np.sin(np.pi * the2 / 180) * np.sin(np.pi * phi2 / 180),
np.cos(np.pi * the2 / 180))
assert_array_almost_equal(sticks, [ref_vec1, ref_vec2])


Expand Down Expand Up @@ -121,7 +121,7 @@ def test_multi_tensor():
s1 = single_tensor(gtab, 100, mevals[0], mevecs[0], snr=None)
s2 = single_tensor(gtab, 100, mevals[1], mevecs[1], snr=None)

Ssingle = 0.5*s1 + 0.5*s2
Ssingle = 0.5 * s1 + 0.5 * s2

S, sticks = MultiTensor(gtab, mevals, S0=100, angles=[(90, 45), (45, 90)],
fractions=[50, 50], snr=None)
Expand All @@ -144,11 +144,11 @@ def test_snr():


def test_all_tensor_evecs():
e0 = np.array([1/np.sqrt(2), 1/np.sqrt(2), 0])
e0 = np.array([1 / np.sqrt(2), 1 / np.sqrt(2), 0])

# Vectors are returned column-wise!
desired = np.array([[1/np.sqrt(2), 1/np.sqrt(2), 0],
[-1/np.sqrt(2), 1/np.sqrt(2), 0],
desired = np.array([[1 / np.sqrt(2), 1 / np.sqrt(2), 0],
[-1 / np.sqrt(2), 1 / np.sqrt(2), 0],
[0, 0, 1]]).T

assert_array_almost_equal(all_tensor_evecs(e0), desired)
Expand All @@ -166,7 +166,7 @@ def test_kurtosis_elements():
[0.00099, 0, 0], [0.00226, 0.00087, 0.00087]])
angles = [(80, 10), (80, 10), (20, 30), (20, 30)]
fie = 0.49 # intra axonal water fraction
frac = [fie * 50, (1-fie) * 50, fie * 50, (1-fie) * 50]
frac = [fie * 50, (1 - fie) * 50, fie * 50, (1 - fie) * 50]
sticks = _check_directions(angles)
mD = np.zeros((len(frac), 3, 3))
for i in range(len(frac)):
Expand All @@ -176,7 +176,7 @@ def test_kurtosis_elements():
# compute global DT
D = np.zeros((3, 3))
for i in range(len(frac)):
D = D + frac[i]*mD[i]
D = D + frac[i] * mD[i]

# compute voxel's MD
MD = (D[0][0] + D[1][1] + D[2][2]) / 3
Expand Down Expand Up @@ -208,7 +208,7 @@ def test_kurtosis_elements():
for j in xyz:
for k in xyz:
for l in xyz:
key = (i+1) * (j+1) * (k+1) * (l+1)
key = (i + 1) * (j + 1) * (k + 1) * (l + 1)
assert_almost_equal(kurtosis_element(mD, frac, i, k, j, l),
kt_ref[key])
# Testing optional funtion inputs
Expand Down Expand Up @@ -288,7 +288,7 @@ def test_DKI_crossing_fibers_simulations():
[0.00099, 0, 0], [0.00226, 0.00087, 0.00087]])
angles = [(80, 10), (80, 10), (20, 30), (20, 30)]
fie = 0.49
frac = [fie*50, (1 - fie)*50, fie*50, (1 - fie)*50]
frac = [fie * 50, (1 - fie) * 50, fie * 50, (1 - fie) * 50]
signal, dt, kt = multi_tensor_dki(gtab_2s, mevals, angles=angles,
fractions=frac, snr=None)
# in this simulations dt and kt cannot have zero elements
Expand Down
32 changes: 15 additions & 17 deletions dipy/sims/voxel.py
View file
Open in desktop
Expand Up @@ -176,9 +176,9 @@ def sticks_and_ball(gtab, d=0.0015, S0=100, angles=[(0, 0), (90, 0)],
sticks = _check_directions(angles)

for (i, g) in enumerate(gtab.bvecs[1:]):
S[i + 1] = f0*np.exp(-gtab.bvals[i + 1]*d) + \
np.sum([fractions[j]*np.exp(-gtab.bvals[i + 1]*d*np.dot(s, g)**2)
for (j, s) in enumerate(sticks)])
S[i + 1] = f0 * np.exp(-gtab.bvals[i + 1] * d) + \
np.sum([fractions[j] * np.exp(-gtab.bvals[i + 1] * d * np.dot(s, g)**2)
for (j, s) in enumerate(sticks)])

S[i + 1] = S0 * S[i + 1]

Expand Down Expand Up @@ -298,16 +298,15 @@ def multi_tensor(gtab, mevals, S0=100, angles=[(0, 0), (90, 0)],
sticks = _check_directions(angles)

for i in range(len(fractions)):
S = S + fractions[i] * single_tensor(gtab, S0=S0, evals=mevals[i],
evecs=all_tensor_evecs(
sticks[i]), snr=None)
S = S + fractions[i] * single_tensor(gtab, S0=S0, evals=mevals[i],
evecs=all_tensor_evecs(
Copy link
Contributor

Choose a reason for hiding this comment

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

put the function and its param on the same line, it's a bit hard to read like that

sticks[i]), snr=None)

return add_noise(S, snr, S0), sticks


def multi_tensor_dki(gtab, mevals, S0=100, angles=[(90., 0.), (90., 0.)],
fractions=[50, 50], snr=20):

r""" Simulate the diffusion-weight signal, diffusion and kurtosis tensors
based on the DKI model

Expand Down Expand Up @@ -386,7 +385,7 @@ def multi_tensor_dki(gtab, mevals, S0=100, angles=[(90., 0.), (90., 0.)],
# compute voxel's DT
DT = np.zeros((3, 3))
for i in range(len(fractions)):
DT = DT + fractions[i]*D_comps[i]
DT = DT + fractions[i] * D_comps[i]
dt = np.array([DT[0][0], DT[0][1], DT[1][1], DT[0][2], DT[1][2], DT[2][2]])

# compute voxel's MD
Expand Down Expand Up @@ -462,7 +461,7 @@ def kurtosis_element(D_comps, frac, ind_i, ind_j, ind_k, ind_l, DT=None,
if DT is None:
DT = np.zeros((3, 3))
for i in range(len(frac)):
DT = DT + frac[i]*D_comps[i]
DT = DT + frac[i] * D_comps[i]

if MD is None:
MD = (DT[0][0] + DT[1][1] + DT[2][2]) / 3
Expand All @@ -471,13 +470,13 @@ def kurtosis_element(D_comps, frac, ind_i, ind_j, ind_k, ind_l, DT=None,

for f in range(len(frac)):
wijkl = wijkl + frac[f] * (
D_comps[f][ind_i][ind_j]*D_comps[f][ind_k][ind_l] +
D_comps[f][ind_i][ind_k]*D_comps[f][ind_j][ind_l] +
D_comps[f][ind_i][ind_l]*D_comps[f][ind_j][ind_k])
D_comps[f][ind_i][ind_j] * D_comps[f][ind_k][ind_l] +
D_comps[f][ind_i][ind_k] * D_comps[f][ind_j][ind_l] +
D_comps[f][ind_i][ind_l] * D_comps[f][ind_j][ind_k])

wijkl = (wijkl - DT[ind_i][ind_j]*DT[ind_k][ind_l] -
DT[ind_i][ind_k]*DT[ind_j][ind_l] -
DT[ind_i][ind_l]*DT[ind_j][ind_k]) / (MD**2)
wijkl = (wijkl - DT[ind_i][ind_j] * DT[ind_k][ind_l] -
DT[ind_i][ind_k] * DT[ind_j][ind_l] -
DT[ind_i][ind_l] * DT[ind_j][ind_k]) / (MD**2)

return wijkl

Expand Down Expand Up @@ -523,7 +522,7 @@ def DKI_signal(gtab, dt, kt, S0=150, snr=None):

# define vector of DKI parameters
MD = (dt[0] + dt[2] + dt[5]) / 3
X = np.concatenate((dt, kt*MD*MD, np.array([np.log(S0)])), axis=0)
X = np.concatenate((dt, kt * MD * MD, np.array([np.log(S0)])), axis=0)

# Compute signals based on the DKI model
S = np.exp(dot(A, X))
Expand All @@ -533,7 +532,6 @@ def DKI_signal(gtab, dt, kt, S0=150, snr=None):
return S



def single_tensor_odf(r, evals=None, evecs=None):
""" Simulated ODF with a single tensor.

Expand Down