gyral_structure.cost package: defining cost function (full documentation)

An overview of the cost functions shown in full detail here can be found here.

Defines individual terms of the cost function.

Normal arithmetic can be used to add the cost function terms together, e.g. Watson + 0.01 * L2

Terms are separated based on what they affect:

• Density constraints are available in gyral_structure.cost.dens

• Orientation constraints are available in gyral_structure.cost.orient

• Constraints that combine other constraints (e.g., smoothness constraints or adding them) are defined in gyral_structure.cost.meta

• Constraints on the parameters rather than the vector field are defined in gyral_structure.cost.param

• gyral_structure.cost.core contains the base class and the cost optimisation code

gyral_structure.cost.core module

Defines subclasses for all cost function terms

• GenericCost serves as the superclass of all cost functions

• CostFunc serves as any cost function on the vector field rather than the parameters directly (see cost.param for the latter)

• SympyCost puts a constraint on the vector field defined by any sympy equation (which is generally easier than defining the cost function and its derivative in python)

In addition this module defines CostEvaluator, which combines any cost function with a basis function. This class contains the actual optimisation code.

class gyral_structure.cost.core.GenericCost

Bases: object

Superclass of all cost objects.

Used to define how cost functions are added together and how they can be combined with a basis function

symmetric = False

Is the cost function symmetric between the reference field and the computed request (should be True to be used in smoothing)?

quadratic = False

Is the cost function quadratic?

independent_pos = True

Is the cost function evaluated for every requested position independently?

epsilon = 1e-06

step size to numerically evaluate the derivative/hessian

get_evaluator(basis_func, method=())

Combines the cost function with a basis function

Parameters

• basis_func – basis function defining the mapping of parameters to a vector field

• method – algorithm used to __call__ the basis function (overrides the algorithm chosen when defining the basis function if provided)

Returns

object that can __call__ the cost function given a set of parameters (and compute derivatives)

describe_fit(predicted_field)

Describes the fit quality.

Parameters

predicted_field – (npos, ndim) array

Returns

statement of the fit quality

class gyral_structure.cost.core.CostFunc(request)

Bases: gyral_structure.cost.core.GenericCost

__init__(request)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

property ndim

property npos

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((xmin, xmax), (ymin, ymax), (zmin, zmax)) covering the full dataset

save(filename)

Saves the cost function to the filename

Parameters

filename – file to store the basis function & solution (if the extension is .gz, .bz, or .bz2 the result will be zipped)

class gyral_structure.cost.core.SympyCost(request)

Bases: gyral_structure.cost.core.CostFunc

Generates CostFunc classes based on a Sympy equation.

simplify = True

field_symbols = {}

__init__(request)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

property funcstr

vectors()

evaluate_sympy()

get_variables(predicted_field)

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

class gyral_structure.cost.core.CostEvaluator(cost_func, basis_func, method=None)

Bases: object

Evaluates the cost function for a given basis function

Variables

evaluator – actually computes the field at the requested positions

__init__(cost_func, basis_func, method=None)

Evaluates a cost function using the basis function.

Parameters

• cost_func (CostFunc) – Cost function describing all constraints

• basis_func (BasisFunc) – Describes the mapping from parameters to the field

property nparams

property ndim

property npos

property quadratic

cost(params)

Computes the cost for a set of parameters

Parameters

params – (nparams, ) array with the parameters for the basis function

Returns

total cost for the given model parameters

dercost(params)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

params – (nparams, ) array with the parameters for the basis function

Returns

(nparams, ) array defining which direction the parameters should move to INCREASE the cost function

hessp(params, dparams)

Compute the multiplication of the hessian (at params) with dparams.

Parameters

• params – (nparams, ) array indicating where the hessian will be calculated

• dparams – (nparams, ) array indicating in which direction the hessian will be calculated

Returns

(nparams, ) of the local hessian times dparams

describe_fit(params)

Prints a description of the fit quality

Parameters

params – (nparams, ) array with the parameters for the basis function

qp(params)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

cvxopt_solve(p0, ftol=0.001, maxiter=1, L2norm=1e-06)

Use cvxopt to exactly solve the quadratic equation.

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((xmin, xmax), (ymin, ymax), (zmin, zmax)) covering the full dataset

fit(params=None, method='L-BFGS-B', options={'maxfun': 50000, 'maxiter': 60000}, tol=1e-08, verbose_every=None)

Fits the basis function to the given cost-function with an added smoothness constraint.

Parameters

• params – initial parameters (default: nearly zero)

• method – name of the optimization algorithm

• options – parameters to pass to the optimization algorithm

• tol – sets the tolerance

• verbose_every – print an update every verbose_every iterations

class gyral_structure.cost.core.VectorCost(request, field=None, norm_length=0, flip=False, mindens=0)

Bases: gyral_structure.cost.core.CostFunc

Computes the offset between the field at the requested point and a reference field

__init__(request, field=None, norm_length=0, flip=False, mindens=0)

Defines a cost function to match the field at the given positions.

The field can represent observed fiber orientations or other constraints (e.g. matching the surface normals). :param request: Defines where the field will be evaluated :param field: (npos, ndim) array with the expected field at the given locations (defaults to 0, i.e. an L2 norm) :param norm_length: Whether to normalize the length of the vectors (0: no normalization; 1: normalize the model field; 2: normalize the target field) :param flip: if True, allow the field to flip sign :param mindens: minimum density to consider

property quadratic

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property symmetric

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

cost(predicted_field, length=None, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field, length=None)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• length – (npos, ) array with the pre-computed expected length of the predicted field (see self.length)

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

length(predicted_field)

Returns the norm that the target field should have to match the observed field

Parameters

predicted_field – (npos, ndim) array with the predicted field at the requested postions

Returns

(npos, ) array

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

gyral_structure.cost.core.read(filename)

Loads a cost function from the given filename

Parameters

filename – file to load the basis function & solution from

gyral_structure.cost.dens module

Defines cost functions that act on the density (i.e. norm) of the vector field, independent of its orientation

• MinInner: Encourages the density of streamlines crossing a surface exceeds a certain value

• TargetInner: Sets a target value for the density of streamlines crossing a surface element

• InnerCost: Sets a target value for the density of streamlines crossing a surface element (same as TargetInner)

• TotalIntersect: Sets a target value for the total number of streamlines crossing a whole surface

• L1Density: L1 norm in the vector field density

• L2Density: L2 norm in the vector field density

• SameDensity: Encourages two fields to have the same vector field density (useful as a smoothness constraint)

class gyral_structure.cost.dens.MinInner(request, field, min_inner=0, normed=False)

Bases: gyral_structure.cost.core.CostFunc

Encourages the inner product to exceed a certain value.

C = 0 if (f cdot F) > min_inner C = min_inner - (f cdot F) if (f cdot F) < min_inner

__init__(request, field, min_inner=0, normed=False)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

property quadratic

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

symmetric = True

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field, tosum=True)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

class gyral_structure.cost.dens.L1Density(request)

Bases: gyral_structure.cost.core.SympyCost

Computes the L1 norm on the fiber density.

This gives a penalty to every streamline proportional to its length

__init__(request)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

class gyral_structure.cost.dens.L2Density(request)

Bases: gyral_structure.cost.core.SympyCost

Computes the L2 norm on the fiber density.

This gives a penalty to every streamline proportional to the streamline length and the surrounding density

__init__(request)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

class gyral_structure.cost.dens.TargetInner(request, field, target_inner=0)

Bases: gyral_structure.cost.core.SympyCost

Encourages a certain number of streamlines to cross the surface.

c = (v - (F cdot f)) ** 2

field_symbols = {'field'}

__init__(request, field, target_inner=0)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

class gyral_structure.cost.dens.SameDensity(request, field, min_norm=1e-20)

Bases: gyral_structure.cost.core.SympyCost

Encourages two fields to have the same density.

c = (F| - |f|)^2 / (|F|^2 + |f|^2 + min_norm)

symmetric = True

field_symbols = {'field'}

__init__(request, field, min_norm=1e-20)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

class gyral_structure.cost.dens.InnerCost(request, field, inner, flip=False)

Bases: gyral_structure.cost.core.CostFunc

Set the inner product between the provided and predicted field.

Can be used to constrain the number of streamlines passing through a normal by setting field to the surface normal and inner to the value. Set flip to True so that the streamlines can pass through in either direction.

quadratic = True

symmetric = True

__init__(request, field, inner, flip=False)

Creates a new constraint on the inner product.

Parameters

• request – Defines where the field will be evaluated

• field – (N, ndim) array defining the field elements

• inner – (N, ) array with the expected inner product

• flip – If True compares the absolute of the inner product rather than the value itself

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

class gyral_structure.cost.dens.TotalIntersect(request, normals, nstream, size=None, flip=False)

Bases: gyral_structure.cost.core.CostFunc

Constrains the total number of streamlines crossing a surface.

quadratic = True

independent_pos = False

__init__(request, normals, nstream, size=None, flip=False)

Creates a new constraint on the total number of streamlines crossing a surface

Parameters

• request – Defines the evaluation of the field across a surface

• normals – (N, ndim) array defining the surface normals

• nstream – total number of streamlines crossing a surface

• size – (N, ) array with the size of each requested element (defaults to request.size())

• flip – Set to True if the nstream can be positive or negative

cost(predicted_field)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

gyral_structure.cost.meta module

Defines cost functions that take other cost functions as input

Most of these can be created by adding/subtracting/multiplying cost functions, except for:

• ExponentAdd: add cost function terms exponentially

• Smoothness: compares the vector field in neighbouring voxels

• MinCost: use the minimum of several cost functions for fitting (useful when fitting to crossing fibers)

class gyral_structure.cost.meta.ExponentAdd(sub_cost)

Bases: gyral_structure.cost.core.CostFunc

\(C = log(\sum_i(exp(c_i))\) for all \(c_i\) in sub_cost

Each element of the requested cost function is added together separately

__init__(sub_cost)

New cost function combining all functions in sub_cost

Parameters

sub_cost (List[CostFunc]) – cost functions to be combined. All cost functions should have the same request

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

property request

describe_fit(predicted_field)

Describes the fit quality.

Parameters

predicted_field – (npos, ndim) array

Returns

statement of the fit quality

class gyral_structure.cost.meta.MinCost(sub_cost)

Bases: gyral_structure.cost.core.CostFunc

Returns the minimum for several cost functions

Useful to fit to crossing fiber data

__init__(sub_cost)

Computes the minimum for every field requested over the provided cost functions

Parameters

sub_cost – list of cost functions with the same request

Raises

ValueError – if the cost functions have different requests

property request

arg_min(predicted_field)

Finds the closest vector to the predicted field.

Parameters

predicted_field – (npos, ndim) array with the field predicted at the requested positions

Returns

tuple with 2 elements:

• (npos, ) integer array with the indices of the minimum cost function (between 0 and number of cost functions -1)

• (npos, ) array with the minimum cost at each requested position

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

describe_fit(predicted_field)

Describes the fit quality.

Parameters

predicted_field – (npos, ndim) array

Returns

statement of the fit quality

class gyral_structure.cost.meta.Smoothness(request, sub_cost, index1, index2)

Bases: gyral_structure.cost.core.CostFunc

Enforces a smoothness contraint.

Applies a cost function to the fields found in neighbouring voxels

__init__(request, sub_cost, index1, index2)

Creates a comparison between the field values at index1 and index2 according to sub_cost.

Parameters

• request (FieldRequest) – Where the field values will be calculated

• sub_cost (CostFunc) – Cost function to compare the field at index1 and index2

• index1 – index array

• index2 – index array

classmethod from_mask(sub_cost, mask, affine=None, jacobian=False, nsplit=3)

Creates a cost function that evaluates the sub_cost in neighbouring voxels in the mask

Parameters

• sub_cost – cost function that expects two fields which will be compared

• mask – (Nx, Ny, Nz) array; the smoothness constraint will be applied wherever the mask in non-zero

• affine – (4, 4) array with the mapping from voxel to mm space

• jacobian – if True evaluates the smoothness constraint on the first derivative of the field rather than the field itself

• nsplit – number of individual elements to split the voxels in

Returns

cost function that compares the field in neighbouring voxels

classmethod from_request(sub_cost, request, max_dist)

Uses sub_cost to compare the fields in any requested positions within min_dist from each other

Parameters

• sub_cost (CostFunc) – cost function that expects two field which will be compared

• request (FieldRequest) – defines where the field will be evaluated

• max_dist – all pairs of requested fields, whose centers are within max_dist from each other will be compared

Returns

cost function that compares the field in neighbouring requested positions

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

class gyral_structure.cost.meta.WeightedCost(sub_cost, weight)

Bases: gyral_structure.cost.core.CostFunc

Returns the cost from the sub_cost multiplied by a float stored in weight

__init__(sub_cost, weight)

Returns a cost-function that is a multiplication of the sub-cost by weight.

Parameters

• sub_cost (CostFunc:) – basis cost function

• weight – number with which the cost function is multiplied.

property request

property quadratic

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

mult_iter(iter)

Helper function to iteratively multiply the result of sub_cost with the given weight.

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

describe_fit(predicted_field)

Describes the fit quality.

Parameters

predicted_field – (npos, ndim) array

Returns

statement of the fit quality

class gyral_structure.cost.meta.PowerCost(sub_cost, power, tosum=False)

Bases: gyral_structure.cost.core.CostFunc

Raises the cost function result in sub_cost to a given power

__init__(sub_cost, power, tosum=False)

Returns the cost-function to the nth power.

Parameters

• sub_cost (CostFunc) – basis cost function

• power – number to which the sub_cost is raised

• tosum – if True sum the cost function before raising it to the nth power

property request

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

class gyral_structure.cost.meta.SumCost(constraints)

Bases: gyral_structure.cost.core.CostFunc

One or more constraints.

__init__(constraints)

Base functionality cost-function.

Parameters

request (FieldRequest) – Defines where the field will be evaluated

property request

property quadratic

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property ndim

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

qp(predicted_field)

Returns the (sparse) matrix \(d\) and vector \(q\) needed to rephrase the problem as a quadratic programming problem.

For a field F, minimize \((1/2) * F' P F + F' q\) where element \(i\), dimension \(d\) is encoded in the vector \(F\) as \(F[i + d * n]\)

Parameters

predicted_field – Returns the quadratic field at the requested position

Returns

tuple with (P, q):

• P: (npos * ndim, npos * ndim) matrix

• q: (npos * ndim, ) vector

add_constraint(new_constraint)

Extends the cost function by adding an additional constraint

describe_fit(predicted_field)

Describes the fit quality.

Parameters

predicted_field – (npos, ndim) array

Returns

statement of the fit quality

gyral_structure.cost.orient module

Defines terms for the cost function that constrain the fibre orientation independent of fiber density

• VonMises: encourages alignment with the direction of a reference field (disallows flipping)

• Watson: encourages alignment with the orientation of a reference field (allows for flipping)

• BinghamDistribution: encourages alignment with a reference field with an ellipsoidal uncertainty (allows for flipping)

class gyral_structure.cost.orient.Watson(request, field, kappa=1.0, min_norm=1e-20, mult_term='None')

Bases: gyral_structure.cost.core.SympyCost

Computes the (unnormalized) -log(p) for the Watson distribution.

field_symbols = {'field'}

__init__(request, field, kappa=1.0, min_norm=1e-20, mult_term='None')

Returns the Watson distribution.

Parameters

• request – where the field will be calculated

• field – fiber orientations to compare the observed field with

• kappa – width of the Watson distribution

• min_norm – correction term to stop the cost function to go to NaN if the requested field is zero

• mult_term –

  how each field-request comparison is weighted

  • ’predicted’: weight with the density of the model field

  • ’field’: weight with the density of the comparison field

  • ’None’: all comparisons have equal weight

  • ’both’: weight with the sum of the model and comparison field densities

property symmetric

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

class gyral_structure.cost.orient.VonMises(request, field, kappa=1.0, min_norm=1e-20, mult_term='None')

Bases: gyral_structure.cost.core.SympyCost

Computes the (unnormalized) -log(p) for the von Mises distribution.

field_symbols = {'field'}

__init__(request, field, kappa=1.0, min_norm=1e-20, mult_term='None')

Returns the Von Mises distribution.

Parameters

• request – where the field will be calculated

• field – fiber orientations to compare the observed field with

• kappa – width of the von Mises distribution

• min_norm – correction term to stop the cost function to go to NaN if the requested field is zero

• mult_term –

  how each field-request comparison is weighted

  • ’predicted’: weight with the density of the model field

  • ’field’: weight with the density of the comparison field

  • ’None’: all comparisons have equal weight

  • ’both’: weight with the sum of the model and comparison field densities

property symmetric

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

class gyral_structure.cost.orient.BinghamDistribution(request, theta, phi, psi, ka, kb, maxk=100000.0, flipx=True)

Bases: gyral_structure.cost.core.CostFunc

Computes the (unnormalized) -log(p) for the Bingham distribution

__init__(request, theta, phi, psi, ka, kb, maxk=100000.0, flipx=True)

Returns the Bingham distribution.

Parameters

• request – where the field will be calculated

• theta – elevation of the primary eigenvector

• phi – orientation of the primary eigenvector in the x-y plane

• psi – determines the orientation of the secondary eigenvector

• ka – inverse disperison along the primary eigenvector

• kb – inverse disperison along the secondary eigenvector

• maxk – maximum value for the inverse dispersion

• flipx – if True flips the x-axis

evaluate_B()

Evaluates the matrix with which the observed field orientation should be multiplied

This method is automatically called whenever maxk is set.

property maxk

Maximum value for the inverse dispersion

static rotation_matrix(angle, around_axis)

Creates a rotation matrix around the given axis

Parameters

• angle – (N, ) array with the angular rotation in radians

• around_axis – around which axis the matrix should rotate

:return (N, 3, 3) array for the rotation around the N angles

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

classmethod read_wb_bingham(filename, voxel_size=0, fmin=0.05, maxk=100000.0, flipx=True)

Reads the Bingham distribution from the CIFTI file produced by the wb_command -estimate-fiber-binghams command.

NOTE: the model orientations are in mm-space not voxel-space like in bedpostx

Parameters

• filename – CIFTI file with the Bingham distributions

• voxel_size – size over which each voxel is averaged

• fmin – minimum mean f to be included in the analysis

• maxk – maximum ka/kb values (can be very high in the results from -estimate-fiber-binghams)

• flipx – if True flip the x-coordinates of the dyad (useful for a negative affine)

class gyral_structure.cost.orient.BinghamDistributionSympy(request, theta, phi, psi, ka, kb)

Bases: gyral_structure.cost.core.SympyCost

Alternative implementation of the Bingham distribution using sympy

Used for comparison with BinghamDistribution (which is faster)

__init__(request, theta, phi, psi, ka, kb)

Returns the Bingham distribution.

Parameters

• request – where the field will be calculated

• theta – elevation of the primary eigenvector

• phi – orientation of the primary eigenvector in the x-y plane

• psi – determines the orientation of the secondary eigenvector

• ka – accuracy along the primary eigenvector

• kb – accuracy along the secondary eigenvector

class gyral_structure.cost.orient.NormedInner(request, field, min_norm=1e-10)

Bases: gyral_structure.cost.core.CostFunc

Encourages alignment with an observed field.

For a model field \(F\) and observed field \(f\) minimizes: \(c = - (F \cdot f)^2 / ((F \cdot F) (f \cdot f))\)

__init__(request, field, min_norm=1e-10)

Creates a new constrained on the normalized inner product

Parameters

• request – Defines where the model field will be evaluated

• field – (N, ndim) array that the model field will be compared against

symmetric = True

cost(predicted_field, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• predicted_field – (npos, ndim) array

• tosum – if False return the cost evaluated for each requested position rather than the full cost

Returns

float, which should be minimized (or (npos, ) array if tosum is False)

dercost(predicted_field)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

predicted_field – (npos, ndim) array

Returns

(npos, ndim) defining which direction the predicted field should move to INCREASE the cost function

hessp(predicted_field, pfield)

Compute the multiplication of the hessian (at predicted_field) with pfield.

Parameters

• predicted_field – (npos, ndim) array indicating where the hessian will be calculated

• pfield – (npos, ndim) array indicating in which direction the hessian will be calculated

Returns

(npos, ndim) of the local hessian times pfield

gyral_structure.cost.param module

Terms of the cost function that take the model parameters rather than the vector field as input

• ParamCost: superclass of all such cost functions

• L2Params: L2 constraint on the parameters

• ParamVecMult: minimize the multiplication of the parameters with a vector

• ParamMatMult: minimize \(p A p - q p\) for some matrix \(A\) and vector \(q\) (for parameters \(p\))

• ParamTargetVec: minimize \(|A p - q|\) for some matrix \(A\) and vector \(q\) (for parameters \(p\))

class gyral_structure.cost.param.ParamCost

Bases: gyral_structure.cost.core.GenericCost

Cost function based on parameters rather than the vector field.

request = ParamRequest

class gyral_structure.cost.param.L2Params

Bases: gyral_structure.cost.param.ParamCost

Computes the L2 norm on the raw parameters rather than the vector field.

cost(params, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• params – (nparams, ) array

• tosum – if False return the cost evaluated for each parameter rather than the full cost

Returns

float, which should be minimized (or (nparams, ) array if tosum is False)

dercost(params)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

params – (nparams, ) array

Returns

(nparams, ) defining which direction the predicted field should move to INCREASE the cost function

class gyral_structure.cost.param.ParamVecMult(vector)

Bases: gyral_structure.cost.param.ParamCost

Minimizes the multiplication of the parameters with some vector of equal length

__init__(vector)

Initialize self. See help(type(self)) for accurate signature.

cost(params, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• params – (nparams, ) array

• tosum – if False return the cost evaluated for each parameter rather than the full cost

Returns

float, which should be minimized (or (nparams, ) array if tosum is False)

dercost(params)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

params – (nparams, ) array

Returns

(nparams, ) defining which direction the predicted field should move to INCREASE the cost function

class gyral_structure.cost.param.ParamMatMult(matrix, vector=None)

Bases: gyral_structure.cost.param.ParamCost

Minimizes the result of matrix multiplication with the parameters

For parameters \(p\) the cost function will return: \(p' \cdot self.matrix \cdot p + self.vector \cdot p\)

__init__(matrix, vector=None)

Creates a new ParamMatMult

Parameters

• matrix – matrix with which the parameter vector will be multiplied twice

• vector – vector with which the parameter vector will be multiplied (defaults to all zeros)

cost(params, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• params – (nparams, ) array

• tosum – if False return the cost evaluated for each parameter rather than the full cost

Returns

float, which should be minimized (or (nparams, ) array if tosum is False)

dercost(params)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

params – (nparams, ) array

Returns

(nparams, ) defining which direction the predicted field should move to INCREASE the cost function

class gyral_structure.cost.param.ParamTargetVec(matrix, vector)

Bases: gyral_structure.cost.param.ParamCost

Minimizes the Euclidean offset between the multiplication of a matrix with the parameters and a target vector

minimizes \(|self.matrix \cdot parameters - self.vector|_2\)

__init__(matrix, vector)

Creates a new ParamTargetVec

Parameters

• matrix – (nvec, nparams) matrix mapping the parameters to a (nvec, ) array

• vector – (nvec, ) target vector

cost(params, tosum=True)

Computes the cost-function given the predicted field at the provided positions.

Parameters

• params – (nparams, ) array

• tosum – if False return the cost evaluated for each vector element rather than the full cost

Returns

float, which should be minimized (or (nparams, ) array if tosum is False)

dercost(params)

Computes the gradient of the cost-function with respect to the predicted field at the provided positions.

Parameters

params – (nparams, ) array

Returns

(nparams, ) defining which direction the predicted field should move to INCREASE the cost function