Skip to content

Commit

Permalink
added brute force example
Browse files Browse the repository at this point in the history
  • Loading branch information
@rutgerfick
rutgerfick committed Feb 1, 2018
1 parent ea629bf commit 889a545
Showing 1 changed file with 248 additions and 0 deletions.
248 changes: 248 additions & 0 deletions examples/example_brute_force_optimization.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,248 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Brute Force Optimization (brute2fine)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One of the main applications of dmipy is the ability to fit any model combination to measured dMRI data. To find model parameters that appropriately describe the underlying tissue configuration, it is important to have optimization algorithms that are able to find the globally best solution based on the given data.\n",
"\n",
"The most straight-forward implementation of global optimization in dmipy is brute-force. There are two implementations:\n",
"- GlobalBruteOptimizer precomputes a global parameter and accompanying signal grid that is re-used when fitting different voxels. It only returns the nearest precomputed parameter combination based on the sum-squared difference between simulated signal grid and the measured data. For fitting many voxels without any prior initial parameter guess this is the fastest way to get a good intial guess. The result can then be passed to the following Brute2FineOptimizer that can refine the solution to a local minimum using gradient-based methods.\n",
"- Brute2FineOptimizer does a separate brute-force optimization for every voxel and then refines the result to a local minimum using gradient-based methods. This approach is slower, but allows for (possibly partial) voxel-varying initial parameter guesses to be passed to the optimizer. Brute2FineOptimizer then only does brute-force optimization on the parameters that did not get an initial guess, and then refines the solution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To illustrate how these optimizers work we will first set up a simple MultiCompartmentMicrostructureModel and load the hcp acquisition scheme as example"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from dmipy.core.modeling_framework import MultiCompartmentModel\n",
"from dmipy.signal_models import gaussian_models\n",
"from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme\n",
"zeppelin = gaussian_models.G2Zeppelin()\n",
"mc_model = MultiCompartmentModel([zeppelin])\n",
"scheme = wu_minn_hcp_acquisition_scheme()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## GlobalBruteOptimizer using precomputed signal grid"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Normally, these optimizers are only called inside the mc_model.fit() function, but we can also call them separately. The GlobalBruteOptimizer has a very straightforward design: given a model and an acquisition scheme, it will equally sample every model parameter between the minimum and maximum parameter ranges in 'Ns' steps. However, it treats orientation parameter 'mu' differently, where it uses a spherical grid of 'N_sphere_samples' points."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from dmipy.optimizers import brute2fine\n",
"Ns = 5 # number of equal sampling points for standard parameters\n",
"N_sphere_samples = 30 # number of orientations that are sampled for orientation 'mu'\n",
"global_brute_optimizer = brute2fine.GlobalBruteOptimizer(mc_model, scheme,\n",
" Ns=Ns, N_sphere_samples=N_sphere_samples)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see it computed a parameter and signal grid for 750 different parameter combinations"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"parameter grid size: (750, 4)\n",
"signal grid size: (750, 288)\n"
]
}
],
"source": [
"print 'parameter grid size:', global_brute_optimizer.parameter_grid.shape\n",
"print 'signal grid size:', global_brute_optimizer.signal_grid.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Showing mc_model's parameters, we can see that this corresponds to $5 * 5 * 30 = 750$ samples for lambda_par, lambda_perp and mu."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OrderedDict([('G2Zeppelin_1_lambda_perp', 1),\n",
" ('G2Zeppelin_1_mu', 2),\n",
" ('G2Zeppelin_1_lambda_par', 1)])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mc_model.parameter_cardinality"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we just take one of the simulated signal grid points as input for the optimizer, we can see that it will return the corresponding parameter combination."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"actual parameters: [ 1.55000000e-09 1.58994250e+00 -2.91158829e+00 1.00000000e-10]\n",
"optimizer output: [ 1.55000000e-09 1.58994250e+00 -2.91158829e+00 1.00000000e-10]\n"
]
}
],
"source": [
"test_signal = global_brute_optimizer.signal_grid[710]\n",
"test_parameters = global_brute_optimizer.parameter_grid[710]\n",
"\n",
"print 'actual parameters:', test_parameters\n",
"print 'optimizer output: ', global_brute_optimizer(test_signal, parameter_scale_normalization=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case of course we get the perfect result. However, for actual data this optimizer only returns precomputed brute-force parameter combinations and does not refine them further. It must be given to the Brute2FineOptimizer to get a better result."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Brute Force with refining"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Brute2FineOptimizer does not precompute anything, and will do the brute force and refining for every voxel. For this reason it is useful when you want to find specific parameter solutions that are close to some initial parameter guess, but not very efficient when the same global grid can be used for every voxel."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"brute2fine_optimizer = brute2fine.Brute2FineOptimizer(mc_model, scheme)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a simple illustration you can give an initial parameter guess x0 and the optimizer will find a local minimum from that point. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"actual parameters: [ 2.27500000e-09 1.00507992e+00 3.40817407e-01 1.00000000e-10]\n",
"start parameters: [ 1.00000000e-10 8.73750549e-01 -8.94552392e-01 1.00000000e-10]\n",
"optimizer output: [ 2.27499914e-09 1.00508076e+00 3.40817944e-01 1.00000000e-10]\n"
]
}
],
"source": [
"test_signal = global_brute_optimizer.signal_grid[90]\n",
"test_parameters = global_brute_optimizer.parameter_grid[90]\n",
"x0_vector = global_brute_optimizer.parameter_grid[400]\n",
"\n",
"print 'actual parameters:', test_parameters\n",
"print 'start parameters: ', x0_vector\n",
"print 'optimizer output: ', brute2fine_optimizer(test_signal, x0_vector) * mc_model.scales_for_optimization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Parts of the x0_vector can also be set to None, in which case the optimizer will first do brute-force for those parameters and then do gradient descent to a local minimum."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

0 comments on commit 889a545

Please sign in to comment.