Reference¶

Variables¶

mnn.model.G = 0.004302¶: float – Gravitational constant to use when evaluating potential or forces on the models. The value must be changed to match the units required by the user.

Exceptions¶

class mnn.model.MNnError(msg)¶: Miyamoto-Nagai negative exceptions : raised when the models parameters are in invalid ranges or that the user is doing something he should not

Miyamoto-Nagai negative model¶

class mnn.model.MNnModel(diz=1.0)¶

Miyamoto-Nagai negative model. This object is a potential-density pair expansion : it consists of a sum of Miyamoto-Nagai dics allowing

__init__(diz=1.0)¶

Constructor for the summed Miyamoto-Nagai-negative model

Parameters:	diz (float) – Normalization factor applied to all the discs (default = 1.0)

add_disc(axis, a, b, M)¶

Adds a Miyamoto-Nagai negative disc to the model, this disc will be included in the summation process when evaluating quantities with the model.

A disc is a list of three parameters a, b and M. All the parameters of the discs are stored in a flat list with no real separation. This is done so that emcee can be fed the array directly without any transformation.

The model accounts for negative values of a. The constraints on the parameters are the following :

b >= 0

M >= 0

a+b >= 0

Parameters:	axis ({‘x’, ‘y’, ‘z’}) – the normal axis of the plane for the disc. a (float) – disc scale b (float) – disc height M (float) – disc mass
Raises:	`mnn.model.MNnError` – if one of the constraints if not satisfied

Example

Adding a disc lying on the xy plane will be done as follows:

>>> m = MNnModel()
>>> m.add_disc('z', 1.0, 0.1, 10.0)

add_discs(values)¶

Wrapper for the add_disc() method to add multiple MNn discs at the same time.

Parameters:	values (list of 4-tuples) – The parameters of the discs to add. One 4-tuple corresponds to one disc.
Raises:	`mnn.model.MNnError` – if one of the constraints if not satisfied

Example

Adding one disc on the xy place with parameters (1.0, 0.1, 50.0) and one disc on the yz plane with parameters (1.0, 0.5, 10.0) will be done as follows:

>>> m = MNnModel()
>>> m.add_discs([('z', 1.0, 0.1, 50.0), ('x', 1.0, 0.5, 10.0)])

static callback_from_string(quantity)¶

Returns the static function callback associated to a given quantity string.

Returns:	One of the following : `mn_density()`, `mn_potential()`, `mn_force()`
Return type:	A function callback

evaluate_density(x, y, z)¶

Evaluates the summed density over all discs at specific positions

Parameters:	x, y, z (float or Nx1 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed density over all discs at position `(x, y, z)`.

Note

If x, y and z are numpy arrays, then the return value is a Nx1 vector of the evaluated potential at every point (x[i], y[i], z[i])

evaluate_density_vec(x)¶

Returns the summed density of all the discs at specific points.

Parameters:	x (Nx3 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed density over all discs at every position in vector `x`.

evaluate_force(x, y, z)¶

Evaluates the summed force over all discs at specific positions

Parameters:	x, y, z (float or Nx1 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed force over all discs at position `(x, y, z)`.

Note

If x, y and z are numpy arrays, then the return value is a Nx3 vector of the evaluated potential at every point (x[i], y[i], z[i])

evaluate_force_vec(x)¶

Returns the summed force of all the discs at specific points.

Parameters:	x (Nx3 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed force over all discs at every position in vector `x`.

evaluate_potential(x, y, z)¶

Evaluates the summed potential over all discs at specific positions

Parameters:	x, y, z (float or Nx1 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed potential over all discs at position `(x, y, z)`.

Note

If x, y and z are numpy arrays, then the return value is a Nx1 value of the potential evaluated at every point (x[i], y[i], z[i])

evaluate_potential_vec(x)¶

Returns the summed potential of all the discs at specific points.

Parameters:	x (Nx3 numpy array) – Cartesian coordinates of the point(s) to evaluate
Returns:	The summed potential over all discs at every position in vector `x`.

generate_dataset_meshgrid(xmin, xmax, nx, quantity='density')¶

Generates a numpy meshgrid of data from the model

Parameters:

xmin (3-tuple of floats) – The low bound of the box
xmax (3-tuple of floats) – The high bound of the box
nx (3-tuple of floats) – Number of points in every direction
quantity ({‘density’, ‘potential’, ‘force’}) – Type of quantity to fill the box with (default=’density’)

Returns:

A 4-tuple containing

vx, vy, vz (N vector of floats): The x, y and z coordinates of each point of the mesh
res (N vector of floats): The values of the summed quantity over all discs at each point of the mesh

Raises:

MemoryError – If the array is too big
mnn.model.MNnError – If the quantity parameter does not correspond to anything known

get_model()¶

Copies the discs currently stored and returns them as a list of 4-tuples [(axis1, a1, b1, M1), (axis2, a2, b2, …), … ]

Returns:	A list of 4-tuples (axis, a, b, M).

Example

>>> m = MNnModel()
>>> m.add_discs([('z', 1.0, 0.1, 50.0), ('x', 1.0, 0.5, 10.0)])
>>> m.get_model()
[('z', 1.0, 0.1, 50.0), ('x', 1.0, 0.5, 10.0)]

static get_tangent_coordinates(x, y, z, axis)¶

Returns the tangent and normal coordinates used in mn_force() from a set of cartesian coordinates and an axis. The correspondence between axis and tangent coordinates are the following :

axis	t1	t2	n
x	y	z	x
y	x	z	y
z	x	y	z

Parameters:

x (float or numpy-array) – x coordinate of the points to convert
y (float or numpy-array) – y coordinate of the points to convert
z (float or numpy-array) – z coordinate of the points to convert
axis ({‘x’, ‘y’, ‘z’}) – the normal axis of the disc

Returns:

t1 (float or numpy-array): The first tangential coordinate for a disc aligned on axis
t2 (float or numpy-array): The second tangential coordinate for a disc aligned on axis
n (float or numpy-array): The normal component for a disc aligned on axis

Return type:

A tuple containing three coordinates

is_positive_definite(max_range=None)¶

Returns true if the sum of the discs are positive definite.

The methods tests along every axis if the minimum of density is positive. If it is not the case then the model should NOT be used since we cannot ensure positive density everywhere.

Parameters:	max_range (a float or None) – Maximum range to evaluate, if None the maximum scale radius will be taken. (default = None)
Returns:	A boolean indicating if the model is positive definite.

static mn_density(r, z, a, b, M)¶

Evaluates the density of a single Miyamoto-Nagai negative disc (a, b, M) at polar coordinates (r, z).

Parameters:	r (float) – radius of the point where the density is evaluated z (float) – height of the point where the density is evaluated a (float) – disc scale b (float) – disc height M (float) – disc mass
Returns:	the density (scaled to the model) at (r, z)
Return type:	float

Note

This method does not check the validity of the constraints b>=0, M>=0, a+b>=0

static mn_force(t1, t2, n, a, b, M, axis)¶

Evaluates the force of a single Miyamoto-Nagai negative disc (a, b, M) at a set of tangent/radial coordinates.

Parameters:	t1 (float) – first tangent coordinate of the point where the density is evaluated t2 (float) – second tangent coordinate of the point where the density is evaluated n (float) – height of the point where the density is evaluated a (float) – disc scale b (float) – disc height Mo (float) – disc mass axis ({‘x’, ‘y’, ‘z’}) – the normal axis of the disc
Returns:	the force applied at point (r, z) relative to the disc in cartesian coordinates.
Return type:	numpy array

Note

This method does not check the validity of the constraints b>=0, M>=0, a+b>=0

Note

This method relies on user-specified value for the gravitational constant. This value can be overriden by setting the value mnn.model.G.

Note

The tangent coordinates allow us to abstract the orientation of the disc to sum everything up for the model. Although it might seem a bit heavy here, it is done to simplify the summation process for the model. Since we require a vector as output we can’t use anymore the “simple” cylindrical coordinates.

The correspondence between axis and tangent coordinates are given in the definition of get_tangent_coordinates()

static mn_potential(r, z, a, b, M)¶

Evaluates the potential of a single Miyamoto-Nagai negative disc (a, b, M) at polar coordinates (r, z).

Parameters:	r (float) – radius of the point where the density is evaluated z (float) – height of the point where the density is evaluated a (float) – disc scale b (float) – disc height Mo (float) – disc mass
Returns:	the potential (scaled to the model) at (r, z)
Return type:	float

Note

This method does not check the validity of the constraints b>=0, M>=0, a+b>=0

Note

This method relies on user-specified value for the gravitational constant. This value can be overriden by setting the value mnn.model.G.

Miyamoto-Nagai negative fitter¶

class mnn.fitter.MNnFitter(n_walkers=100, n_steps=1000, n_threads=1, random_seed=123, fit_type='density', check_positive_definite=False, cdp_range=None, allow_negative_mass=False, verbose=False)¶

Miyamoto-Nagai negative fitter.

This class is used to fit the parameters of a Miyamoto-Nagai negative model to data. (with a predefined number of discs) to a datafile.

__init__(n_walkers=100, n_steps=1000, n_threads=1, random_seed=123, fit_type='density', check_positive_definite=False, cdp_range=None, allow_negative_mass=False, verbose=False)¶

Constructor for the Miyamoto-Nagai negative fitter. The fitting is based on emcee.

Parameters:

n_walkers (int) – How many parallel walkers emcee will use to fit the data (default=100).
n_step (int) – The number of steps every walker should perform before stopping (default=1000).
n_threads (int) – Number of threads used to fit the data (default=1).
random_seed (int) – The random seed used for the fitting (default=123).
fit_type ({‘density’, ‘potential’}) – What type of data is fitted (default=’density’).
check_positive_definite (bool) – Should the algorithm check if every walker is positive definite at every step ?
cdp_range ({float, None}) – Maximum range to which check positive definiteness. If none, the criterion will be tested, for each axis on 10*max_scale_radius
allow_negative_mass (bool) – Allow the fitter to use models with negative masses (default=False)
verbose (bool) – Should the program output additional information (default=False).

Note

Using check_positive_definite=True might guarantee that the density will be always positive. But this takes a toll on the computation. We advise to fit the data with check_positive_definite=False. If the result is not definite positive, then switch this flag on and re-do the fitting.

corner_plot(model=None)¶

Computes the corner plot of the fitted data.

Note

If this method fails it might mean the fitting has not properly converged yet.

Parameters:	model (tuple) – A flattened model or None. If None no truth value will be displayed on the plot.
Returns:	The corner plot object.

fit_data(burnin=100, x0=None, x0_range=0.01, plot_freq=0, plot_ids=[])¶

Runs emcee to fit the model to the data.

Fills the mnn.fitter.sampler object with the putative models and returns the burned-in data. The walkers are initialized randomly around position x0 with a maximum dispersion of x0_range. This ball is the initial set of solutions and should be centered on the initial guess of what the parameters are.

Parameters:

burnin (int) – The number of timesteps to remove from every walker after the end (default=100).
x0 (numpy array) – The initial guess for the solution (default=None). If None, then x0 is determined randomly.
x0_range (float) – The radius of the inital guess walker ball. Can be either a single scalar or a tuple of size 3*n_discs (default=1e-2)
plot_freq (int) – The frequency at which the system outputs control plot (default=0). If 0, then the system does not plot anything until the end.
plot_ids (array) – The id of the discs to plot during the control plots (default=[]). If empty array, then every disc is plotted.

Returns:

A tuple containing

samples (numpy array): A 2D numpy array holding every parameter value for every walker after timestep burnin
lnprobability (numpy array): The samplers pointer to the matrix value of the log likelihood produced by each walker at every timestep after burnin

Raises:

MNnError – If the user tries to fit the data without having called load_data() before.

Note

The plots are outputted in the folder where the script is executed, in the file current_state.png.

get_residuals(model)¶

Computes the residual between the data and the model you provide as input

Parameters:	model (numpy array) – The Ndiscs*3 parameter values of the model you want to compute the residuals on.
Returns:	A numpy array storing the residual value for every point of the data.
Raises:	`MNnError` – If the user tries to compute the residual without having called `load_data()` before.

load_data(filename)¶

Loads the data that will be fitted to the model.

The data should be in an ascii file with four columns tab or space separated : X Y Z quantity

Parameters:	filename (string) – The filename to open.

loglikelihood(discs)¶

Computes the log likelihood of a given model

Parameters:	discs (tuple) – the list of parameters for the model stored in a flat-tuple (a1, b1, M1, a2, b2, …)
Returns:	The loglikelihood of the model given in parameter

make_model(model)¶

Takes a flattened model as parameter and returns a mnn.model.MNnModel object.

Parameters:	model (a tuple or a numpy object) – The flattened model
Returns:	A `mnn.model.MNnModel` instance corresponding to the flattened model

maximum_likelihood()¶

Computation of the maximum likelihood for a given model and stores them in MNnFitter.model

Parameters:	samples – The nd-array given by the fit_data routine
Returns:	The parameters corresponding to the maximized log likelihood

plot_disc_walkers(id_discs=None)¶

Plotting the walkers on each parameter of a certain disc.

Parameters:	id_disc (int of list) – the ids of the disc parameters you want to plot. If None, all the discs are plotted
Returns:	The matplotlib figure object. You can either plot it or save it.

set_model_type(nx=0, ny=0, nz=1)¶

Defines the type of Miyamoto-nagai negative model that will be fitted

This method allows to set the number of discs available to put in the model along each plane.

Parameters:	nx (int) – Number of discs on the yz plane (default=0). ny (int) – Number of discs on the xz plane (default=0). nz (int) – Number of discs on the xy plane (default=1).