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