from dmipy.data import saved_acquisition_schemes
scheme_hcp = saved_acquisition_schemes.wu_minn_hcp_acquisition_scheme()
from dmipy.signal_models import cylinder_models, gaussian_models, sphere_models
from dmipy.distributions.distribute_models import SD1WatsonDistributed
from dmipy.core.modeling_framework import MultiCompartmentModel
lambda_iso_diff = 3.e-9
lambda_par_diff = 1.7e-9
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
watson_dispersed_bundle = SD1WatsonDistributed(models=[stick,zeppelin])
watson_dispersed_bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp', 'C1Stick_1_lambda_par', 'partial_volume_0')
watson_dispersed_bundle.set_equal_parameter('G2Zeppelin_1_lambda_par', 'C1Stick_1_lambda_par')
watson_dispersed_bundle.set_fixed_parameter('G2Zeppelin_1_lambda_par', lambda_par_diff)
NODDI_mod = MultiCompartmentModel(models=[ball, watson_dispersed_bundle])
NODDI_mod.set_fixed_parameter('G1Ball_1_lambda_iso', lambda_iso_diff)
import numpy as np
Nt = 12
N_models = len(NODDI_mod.models)
def forward_model_matrix(acquisition_scheme, model_dirs):
dir_params = [p for p in NODDI_mod.parameter_names if p.endswith('mu')]
if len(dir_params) != len(model_dirs):
raise ValueError("Length of model_dirs should correspond "
"to the number of directional parameters!")
grid_params = [p for p in NODDI_mod.parameter_names
if not p.endswith('mu') and not p.startswith('partial_volume')]
_amico_grid, _amico_idx = {}, {}
# Compute length of the vector x0
x0_len = 0
for m_idx in range(N_models):
m_atoms = 1
for p in NODDI_mod.models[m_idx].parameter_names:
if NODDI_mod.model_names[m_idx] + p in grid_params:
m_atoms *= Nt
x0_len += m_atoms
for m_idx in range(N_models):
model = NODDI_mod.models[m_idx]
model_name = NODDI_mod.model_names[m_idx]
param_sampling, grid_params_names = [], []
m_atoms = 1
for p in model.parameter_names:
if model_name + p not in grid_params:
continue
grid_params_names.append(model_name + p)
p_range = NODDI_mod.parameter_ranges[model_name + p]
_amico_grid[model_name + p] = np.full(x0_len, np.mean(p_range))
param_sampling.append(np.linspace(p_range[0], p_range[1], Nt, endpoint=True))
m_atoms *= Nt
_amico_idx[model_name] =\
sum([len(_amico_idx[k]) for k in _amico_idx]) + np.arange(m_atoms)
params_mesh = np.meshgrid(*param_sampling)
for p_idx, p in enumerate(grid_params_names):
_amico_grid[p][_amico_idx[model_name]] = np.ravel(params_mesh[p_idx])
_amico_grid['partial_volume_' + str(m_idx)] = np.zeros(x0_len)
_amico_grid['partial_volume_' + str(m_idx)][_amico_idx[model_name]] = 1.
for d_idx, dp in enumerate(dir_params):
_amico_grid[dp] = model_dirs[d_idx]
return NODDI_mod.simulate_signal(acquisition_scheme,_amico_grid).T, _amico_grid, _amico_idx
np.random.seed(123)
n_samples = 100
ods = np.random.uniform(0.03, 0.99, n_samples)
ic_vfs = np.random.uniform(0.1, 0.99, n_samples)
iso_vfs = np.random.uniform(0, 1., n_samples)
theta = np.random.uniform(0, np.pi, n_samples)
phi = np.random.uniform(-np.pi, np.pi, n_samples)
arguments = dict.fromkeys(NODDI_mod.parameter_names)
arguments['SD1WatsonDistributed_1_SD1Watson_1_odi'] = ods
arguments['SD1WatsonDistributed_1_partial_volume_0'] = ic_vfs
arguments['partial_volume_0'] = iso_vfs
arguments['partial_volume_1'] = 1. - iso_vfs
arguments['SD1WatsonDistributed_1_SD1Watson_1_mu'] = np.column_stack((theta, phi))
signals = NODDI_mod.simulate_signal(scheme_hcp, arguments)
print(signals.shape)
from dmipy.optimizers import amico_cvxpy
lambda_1 = [0, 0.0001]
lambda_2 = [0, 0.0001]
x0_vector = np.zeros(145)
amico_opt = amico_cvxpy.AmicoCvxpyOptimizer(NODDI_mod, scheme_hcp, x0_vector=x0_vector,
lambda_1=lambda_1, lambda_2=lambda_2)
for i in range(n_samples):
data = signals[i, :]
mu = [theta[i], phi[i]]
M, grid, idx = forward_model_matrix(scheme_hcp, [mu])
parameters = amico_opt(data, M, grid, idx)
print("estimated:", parameters)
print("ground tr:", iso_vfs[i], 1- iso_vfs[i], ods[i], ic_vfs[i])
print("\n")