kdsource Python API

Python module for creating and optimizing KDSource objects.

KDSource’s are particle sources for Monte Carlo radiation transport simulations. The full distribution and documentation can be found in the project GitHub page:

https://github.com/KDSource/KDSource

A KDSource object is based on a particle list in MCPL format (see https://mctools.github.io/mcpl/), on which a Kernel Density Estimation (KDE) is applied, using KDEpy library (see https://github.com/tommyod/KDEpy).

With kdsource Python library you can create, optimize, and export KDSource objects as XML files. These files can be later used as distributional sources in other Monte Carlo simulations, using the tools available in the full distribution.

If you use KDSource tools in your work, please add the following reference:

Abbate, O. I., Schmidt, N. S., Prieto, Z. M., Robledo, J. I., Dawidowski, J., Márquez, A. A., & Márquez Damián, J. I. KDSource, a tool for the generation of Monte Carlo particle sources using kernel density estimation [Computer software]. https://github.com/KDSource/KDSource

kdsource.kdsource module

Module for KDSource object

kdsource.kdsource.load(xmlfilename, N=-1)

Load KDSource from XML parameters file. After building the KDSource object, a fit is performed to load the particle list, leaving the KDSource ready to be evaluated. The bandwidths given in the XML file are not modified.

Parameters

xmlfilename: str: Name of XML parameters file
N: int: Number of particles to use for fitting. Default is -1, meaning that all particles in the list will be used.

Returns

kdsource: KDSource object

class kdsource.kdsource.KDSource(plist, geom, bw='silv', J=1.0)

Bases: object

__init__(plist, geom, bw='silv', J=1.0)

Object representing Kernel Density Estimation (KDE) sources.

A KDSource object is based on a particle list in MCPL format, wrapped in a PList object. It also includes a Geometry object, which defines variables (energy, position, direction) treatment. Having these basic components, KDSource can fit (optimize) a KDE model on the particle list, creating a distributional source.

After optimization, analysis, plotting, etc., a KDSource can be exported in XML format, for later use in Python or other KDSource tools.

Parameters

plist: PList object: The PList wrapping the MCPL file containing the particle list.
geom: Geometry object: The Geometry defining particle variables treatment.
bw: float, str or numpy.ndarray, optional: Bandwidth or bandwidth selection method for KDE. If a float is passed, it is the bandwidth for all particles. If a numpy.ndarray is passed, it is the bandwidth for each particle. If a string is passed it is the bandwidth selection method. See bw_methods for available methods.
J: float, optional: The source total current, in [1/s]. If set, the density plots will have the correct units.

fit(N=-1, skip=0, scaling=None, **kwargs)

Fit KDE model to particle list.

A number of particles are loaded from MCPL file to build the KDE model. If the KDSource object has a bandwidth selection method, it will be applied to optimize the bandwidth.

Parameters

N: int: Number of particles to use for fitting. The real number of particles used may be lower if end of particle list is reached or there are particles with zero weight. -1 to use all particles.
skip: int: Number of particles to skip in the list before starting to read.
scaling: array-like, optional: Scaling to be applied to each variable. This means each particle variable will be divided by the corresponding scaling element before applying KDE. By default, the standard deviation of each variable is used.
**kwargs: optional: Parameters to be passed to bandwidth selection method. Refer corresponding method for docs (see bw_methods for method names).

evaluate(parts)

Evaluate density and statistic error in a set of points.

Variables and scaling treatment is so that the evaluated density units will be the inverse of the product of metrics volunits parameter, times [s-1] (e.g.: [s-1 MeV-1 cm-2 sr-1]).

Parameters

parts: array-like: Array of particles where evaluate density. Must have shape (obs, 8), with columns being [ekin, x, y, z, dx, dy, dz, t].

Returns

[evals, errs]: list of arrays: Array of evaluated densities and estimated statistic errors.

save(xmlfilename=None, bwfile=None, adjust_N=True)

Save KDSource object to XML parameters file.

Parameters

xmlfilename: str: Name of XML parameters file. By default it is set as: [MCPL file name]_source.xml
bwfile: str or file object: File name or file object to write bandwidths in binary float32 format, if a variable bandwidth is used. A file object with append mode can be used to save several KDSource’s in same parameters file, useful when particle list is so big it can’t be loaded at once.
adjust_N: bool: If True, before saving the following factor is applied on bandwidth: bw_silv(dim, N_tot) / bw_silv(dim, N) where N_tot is the total number of particles in the MCPL file, and N is the one passed to fit method. This way the bandwidth optimized for a subset of the particle list can be adapted to the full list.

Returns

xmlfilename: str: Name of the written XML file.

plot_point(var, grid, part0, **kwargs)

1D plot of the full multivariate distribution.

The current densities plotted are evaluated with the evaluate method, along a grid of particles, varying only one variable.

Parameters

var: int or str

Variable to be plotted, or its index. Names and indices of variables can be found in varnames list.

grid: array-like

1D array with values of var variable.

part0: array-like

Particle defining where to evaluate the rest of variables. Must have shape (7,), with variables ordered as in varnames list. Value of var variable will be ignored.

**kwargs: optional

Additional parameters for plotting options:

xscale: ‘linear’ or ‘log’: Scale for x axis. Default: ‘linear’
yscale: ‘linear’ or ‘log’. Default: ‘log’: Scale for y axis.
fact: float: Factor to apply on all densities. Default: 1
label: str: Line label in plot legend.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

plot_integr(var, grid, vec0=None, vec1=None, **kwargs)

1D plot of the univariate distribution along one variable.

Densities are integrated over the specified range of the other variables, and evaluated on a newly created 1D KDE model. Units are the inverse of the units of the parametrized chosen variable (see geometry units parameter), times [s-1] (e.g.: [s-1 cm-1], [s-1 deg-1]).

Parameters

var: int or str

Variable to be plotted, or its index. Names and indices of variables can be found in varnames parameter of geometry. Note that this is a parametrized variable (e.g.: lethargy, theta, etc.).

grid: array-like

1D array with values of var variable.

vec0: array-like

Lower limits of integration range for other variables. Must be of shape (geom.dim,), with variables ordered as in varnames parameter of the geometry. Value of var variable will be ignored.

vec1: array-like

Upper limits of integration range for other variables. Same considerations as for vec0.

**kwargs: optional

Additional parameters for plotting options:

xscale: ‘linear’ or ‘log’: Scale for x axis. Default: ‘linear’
yscale: ‘linear’ or ‘log’. Default: ‘log’: Scale for y axis.
fact: float: Factor to apply on all densities. Default: 1
label: str: Line label in plot legend.
adjust_bw: str: If True, a Silverman factor will be applied to adjust the bandwidth to the dimension and number of samples used in plot. Default: False.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

plot_E(grid_E, vec0=None, vec1=None, **kwargs)

1D plot of the energy spectrum.

Same as plot(“ekin”, grid_E), but applying jacobian of energy parametrization (if any). When using Lethargy metric, plot method will give plots with units [s-1] (lethargy is dimesionless), while this method will produce the expected spectrum with units [s-1 MeV-1].

Parameters

grid_E: array-like

1D array with values of ekin.

vec0: array-like

Lower limits of integration range for other variables. Must be of shape (geom.dim,), with variables ordered as in varnames parameter of the geometry. Value of energy variable will be ignored.

vec1: array-like

Upper limits of integration range for other variables. Same considerations as for vec0.

**kwargs: optional

Additional parameters for plotting options:

fact: float: Factor to apply on all densities. Default: 1
label: str: Line label in plot legend.
adjust_bw: str: If True, a Silverman factor will be applied to adjust the bandwidth to the dimension and number of samples used in plot. Default: False.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

plot_t(grid_t, vec0=None, vec1=None, **kwargs)

1D plot of the time distribution.

Same as plot(“t”, grid_t), but applying jacobian of time parametrization (if any). When using Decade metric, plot method will give plots with units [s-1] (decade is dimesionless), while this method will produce the expected spectrum with units [s-2].

Parameters

grid_t: array-like

1D array with values of t.

vec0: array-like

Lower limits of integration range for other variables. Must be of shape (geom.dim,), with variables ordered as in varnames parameter of the geometry. Value of energy variable will be ignored.

vec1: array-like

Upper limits of integration range for other variables. Same considerations as for vec0.

**kwargs: optional

Additional parameters for plotting options:

fact: float: Factor to apply on all densities. Default: 1
label: str: Line label in plot legend.
adjust_bw: str: If True, a Silverman factor will be applied to adjust the bandwidth to the dimension and number of samples used in plot. Default: False.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

plot2D_point(vrs, grids, part0, **kwargs)

2D plot of the full multivariate distribution.

The current densities plotted are evaluated with the evaluate method, along a grid of particles, varying only two variable.

Parameters

vrs: list

List of the 2 variables to be plotted, or its indices. Names and indices of variables can be found in varnames list.

grids: list

List of the 2 1D arrays with values of vrs variables.

part0: array-like

Particle defining where to evaluate the rest of variables. Must have shape (7,), with variables ordered as in varnames list. Value of vrs variables will be ignored.

**kwargs: optional

Additional parameters for plotting options:

scale: ‘linear’ or ‘log’: Scale for color map. Default: ‘log’
fact: float: Factor to apply on all densities. Default: 1
title: str: Plot title.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

__dict__ = mappingproxy({'__module__': 'kdsource.kdsource', '__init__': <function KDSource.__init__>, 'fit': <function KDSource.fit>, 'evaluate': <function KDSource.evaluate>, 'save': <function KDSource.save>, 'plot_point': <function KDSource.plot_point>, 'plot_integr': <function KDSource.plot_integr>, 'plot_E': <function KDSource.plot_E>, 'plot_t': <function KDSource.plot_t>, 'plot2D_point': <function KDSource.plot2D_point>, 'plot2D_integr': <function KDSource.plot2D_integr>, '__dict__': <attribute '__dict__' of 'KDSource' objects>, '__weakref__': <attribute '__weakref__' of 'KDSource' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.kdsource'

__weakref__: list of weak references to the object (if defined)

plot2D_integr(vrs, grids, vec0=None, vec1=None, **kwargs)

2D plot of the bivariate distribution along two variables.

Densities are integrated over the specified range of the other variables, and evaluated on a newly created 2D KDE model. Units are the inverse of the units of the parametrized chosen variables (see geometry units parameter), times [s-1] (e.g.: [s-1 cm-2], [s-1 deg-2]).

Parameters

vrs: list

List of the 2 variables to be plotted, or its indices. Names and indices of variables can be found in varnames parameter of geometry. Note that these are parametrized variables (e.g.: lethargy, theta, etc.).

grids: list

1D array with values of vrs variable.

vec0: array-like

Lower limits of integration range for other variables. Must be of shape (geom.dim,), with variables ordered as in varnames parameter of the geometry. Value of vrs variable will be ignored.

vec1: array-like

Upper limits of integration range for other variables. Same considerations as for vec0.

**kwargs: optional

Additional parameters for plotting options:

scale: ‘linear’ or ‘log’: Scale for color map. Default: ‘log’
fact: float: Factor to apply on all densities. Default: 1
title: str: Plot title.
adjust_bw: str: If True, a Silverman factor will be applied to adjust the bandwidth to the dimension and number of samples used in plot. Default: False.

Returns

[fig, [scores, errs]]:: Figure object, and evaluated densities and statistic errors.

kdsource.plist module

Module for PList object

kdsource.plist.convert2mcpl(filename, readformat)

Convert particle list with MCPL-compatible format to MCPL format.

If there is a file with same name as the file to convert, and MCPL format, it will be assumed that particle list has already been converted.

Conversion is executed with subprocess.run, calling the corresponding conversion executable. For this to work, KDSource binaries must be in PATH.

Parameters

filename: str

Name of file with particle list to convert.

readformat: str

Particle list format. Valid formats are:

‘ssw’: MCNP binary surface source.
‘phits’: PHITS particle list.
‘ptrac’: MCNP text surface source.
‘stock’: TRIPOLI-4 stored particles.
‘ssv’: Text file with space-separated values, with format resulting from ‘mcpltool –text’ command.

Returns

mcplname: str: Name of the created or existing MCPL file.

kdsource.plist.join2mcpl(filenames, readformat)

Merge particle lists with MCPL-compatible format into MCPL file.

Each file is converted to MCPL with convert2mcpl method, and then merged with mcpltool command.

Merged MCPL file name is constructed as the intersection of filenames, if any, or set as ‘merged.mcpl.gz’ otherwise.

Parameters

filenames: list: Names of files with particle lists to convert.
readformat: str: Particle list format for all files. See convert2mcpl for valid formats.

Returns

mergedname: str: Name of the created MCPL file.

kdsource.plist.savessv(pt, parts, ws, outfile)

Save particle list to Space-Separated Values file.

This function is equivalent to convert2ascii from mcpl, and results in the same format as ‘mcpltool –text’ command. Generated SSV file can be converted to MCPL format with convert2mcpl.

Parameters

pt: str: Particle type. See pt2pdg for available particle types.
parts: array-like: Array of particles. Must have shape (N, 8), with columns ordered as in varnames list.
ws: array-like: Particle weights.
outfile: str: Name of SSV file. If it exist, its content will be overwritten.

kdsource.plist.appendssv(pt, parts, ws, outfile)

Append particle list to Space-Separated Values file.

This function uses the same format as savessv, but appends particles without writing a header.

Parameters

pt: str: Particle type. See pt2pdg for available particle types.
parts: array-like: Array of particles. Must have shape (N, 8), with columns ordered as in varnames list.
ws: array-like: Particle weights.
outfile: str: Name of SSV file. If it exist, its content will be overwritten.

class kdsource.plist.PList(filename, readformat='mcpl', pt='n', trasl=None, rot=None, switch_x2z=False, set_params=True)

Bases: object

__init__(filename, readformat='mcpl', pt='n', trasl=None, rot=None, switch_x2z=False, set_params=True)

Object defining particle list. It is a wrapper for a MCPL file.

It also allows performing a spatial transformation right after reading particles (during get method), which is useful for using particle lists in simulations with different reference systems.

Parameters

filename: str or list: File or list of files containing the particle list, in any MCPL-compatible format.
readformat: str: Particle list format. Can be ‘mcpl’ or any valid value of convert2mcpl readformat parameter.
pt: str: Particle type. See pt2pdg for available particle types.
trasl: array-like, optional: Spatial translation. Default is no translation.
rot: numpy.ndarray or scipy.spatial.transform.Rotation, optional: Spatial rotation. Can be a scipy Rotation object, or any array-like which can be used to generate a scipy Rotation object (rotation vector, quaternion or rotation matrix). Rotation is applied after translation.
switch_x2z: bool: If true, the following permutation is applied after translation and rotation:

(x, y, z) -> (y, z, x)
set_params: bool: Whether to set I (sum of weights) and p2 (sum of squared weights) after construction.

set_params()

Set parameters particle list parameters.

The following parameters are set:: sum_seights: Sum of weights p2: Sum of square weights N: Total number of valid particles in list N_eff = I^2/p2: Effective number of particles

get(N=-1, skip=0)

Get particles from MCPL file.

After loading particles, translation, rotation and switch_x2z transformations will be applied (if any).

Parameters

N: int: Number of particles to read. The real number N’ of particles read may be lower if end of particle list is reached or there are particles with zero weight. -1 to read all particles.
skip: int: Number of particles to skip in the list before starting to read.

Returns

[parts,ws]: list: Array of particles and weights. Particles array will have shape (N’, 8), with columns ordered as in varnames list.

save(pltree): Save PList parameters into XML tree.

__dict__ = mappingproxy({'__module__': 'kdsource.plist', '__init__': <function PList.__init__>, 'set_params': <function PList.set_params>, 'get': <function PList.get>, 'save': <function PList.save>, 'load': <staticmethod object>, '__dict__': <attribute '__dict__' of 'PList' objects>, '__weakref__': <attribute '__weakref__' of 'PList' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.plist'

__weakref__: list of weak references to the object (if defined)

static load(pltree): Load parameters from XML tree and build PList.

kdsource.geom module

Module for Geometry and Metric objects

class kdsource.geom.Metric(partvars, varnames, units, volunits)

Bases: object

__init__(partvars, varnames, units, volunits)

Abstract object defining metrics for a subset of variables.

The main function of a Metric is definig a transformation from a subset of particle variables to certain parametrized variables, which can be more suitable for applying KDE.

See geom._metrics for available metrics.

Parameters

partvars: list of int: Indices of the particle variables to parametrize (see varnames global list).
varnames: list of str: Names of parametrized variables.
units: list of str: Units of each parametrized variable.
volunits: str: Units of the product of variables to parametrize.

transform(parts): Transform particle variables to parametrized variables.

inverse_transform(vecs): Transform parametrized variables to particle variables.

jac(parts): Jacobian of transformation.

mean(parts=None, vecs=None, weights=None)

Mean of particle variables.

Mean is computed in parametrized space, and transformed back to particle variables.

Parameters

parts: array-like, optional: Array of particle variables.
vecs: array-like, optional: Array of parametrized variables. If set, overrides parts.
weights: array-like, optional: Array of particle statistic weights.

std(parts=None, vecs=None, weights=None)

Standard deviation of particle variables.

Standard deviation is computed in parametrized space, and transformed back to particle variables.

Parameters

parts: array-like, optional: Array of particle variables.
vecs: array-like, optional: Array of parametrized variables. If set, overrides parts.
weights: array-like, optional: Array of particle statistic weights.

save(mtree): Save Metric parameters into XML tree.

static load(mtree): Load parameters from XML tree and build Metric.

__dict__ = mappingproxy({'__module__': 'kdsource.geom', '__init__': <function Metric.__init__>, 'transform': <function Metric.transform>, 'inverse_transform': <function Metric.inverse_transform>, 'jac': <function Metric.jac>, 'mean': <function Metric.mean>, 'std': <function Metric.std>, 'save': <function Metric.save>, 'load': <staticmethod object>, '__dict__': <attribute '__dict__' of 'Metric' objects>, '__weakref__': <attribute '__weakref__' of 'Metric' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.geom'

__weakref__: list of weak references to the object (if defined)

class kdsource.geom.Geometry(metrics, trasl=None, rot=None)

Bases: Metric

__init__(metrics, trasl=None, rot=None)

Object defining particle variables treatment (metrics).

The main function of a Geometry is definig a transformation from particle variables to certain parametrized variables, which can be more suitable for applying KDE.

Parameters

metrics: list: List of metrics for each subset of variables. They must cover all particle variables (energy, position, direction and time).
trasl: array-like, optional: Spatial traslation for source. Default is no traslation.
rot: numpy.ndarray or scipy.spatial.transform.Rotation, optional: Spatial rotation for source. Can be a scipy Rotation object, or any array-like which can be used to generate a scipy Rotation object (rotation vector, quaternion or rotation matrix). Rotation is applied after traslation in transform, and before traslation in inverse_transform. Default is no rotation.

transform(parts): Transform particle variables to parametrized variables.

inverse_transform(vecs): Transform parametrized variables to particle variables.

jac(parts): Jacobian of transformation.

mean(parts=None, vecs=None, weights=None)

Mean of particle variables.

Mean is computed in parametrized space, and transformed back to particle variables.

Parameters

parts: array-like, optional: Array of particle variables.
vecs: array-like, optional: Array of parametrized variables. If set, overrides parts.
weights: array-like, optional: Array of particle statistic weights.

std(parts=None, vecs=None, weights=None)

Standard deviation of particle variables.

Standard deviation is computed in parametrized space, and transformed back to particle variables.

Parameters

parts: array-like, optional: Array of particle variables.
vecs: array-like, optional: Array of parametrized variables. If set, overrides parts.
weights: array-like, optional: Array of particle statistic weights.

save(gtree): Save Geometry parameters into XML tree.

static load(gtree): Load parameters from XML tree and build Geometry.

__module__ = 'kdsource.geom'

class kdsource.geom.Energy

Bases: Metric

__init__(): Simple metric for energy, with no transformation.

load(): Build Energy.

__module__ = 'kdsource.geom'

class kdsource.geom.Lethargy(E0=10)

Bases: Metric

__init__(E0=10)

Lethargy metric for energy.

Lethargy is defined as:: u = log(E0 / ekin)

Parameters

E0: float: Reference energy. Typically, it is the highest energy in the system.

transform(ekins): Transform energy to lethargy.

inverse_transform(us): Transform lethargy to energy.

jac(ekins): Jacobian of lethargy transformation.

save(mtree): Save Lethargy parameters into XML tree.

static load(mtree): Load parameters from XML tree and build Lethargy.

__module__ = 'kdsource.geom'

class kdsource.geom.Time

Bases: Metric

__init__(): Simple metric for time, with no transformation.

load(): Build Time.

__module__ = 'kdsource.geom'

class kdsource.geom.Decade

Bases: Metric

__init__()

Decade metric for time.

Decade is defined as:: d = log(t)

transform(ts): Transform time to decade.

inverse_transform(ds): Transform deacde to time.

jac(ts): Jacobian of lethargy transformation.

static load(mtree): Build Decade.

__module__ = 'kdsource.geom'

class kdsource.geom.Vol(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, zmin=-inf, zmax=inf)

Bases: Metric

__init__(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, zmin=-inf, zmax=inf)

Simple metric for 3D position, with no transformation.

Each spatial variable (x, y, z) is delimited between a min and max value. By default these are -infinity and infinity, respectively. All positions in the particle list should be inside these limits.

Axis system of reference can be changed by means of the ‘trasl’ and ‘rot’ arguments of the Geometry object.

save(mtree): Save Vol parameters into XML tree.

static load(mtree): Load parameters from XML tree and build Vol.

__module__ = 'kdsource.geom'

class kdsource.geom.SurfXY(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, z=0)

Bases: Metric

__init__(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, z=0)

Simple metric for 2D position, with no transformation.

Spatial variables x and y are delimited between a min and max value. By default these are -infinity and infinity, respectively. All positions in the particle list should be inside these limits.

z has the fixed value given as argument.

Axis system of reference can be changed by means of the ‘trasl’ and ‘rot’ arguments of the Geometry object.

transform(poss): Transform volume position (x,y,z) to flat position (x,y).

inverse_transform(poss): Transform flat position (x,y) to volume position (x,y,z).

save(mtree): Save SurfXY parameters into XML tree.

static load(mtree): Load parameters from XML tree and build SurfXY.

__module__ = 'kdsource.geom'

class kdsource.geom.Guide(xwidth, yheight, zmax=inf, rcurv=None)

Bases: Metric

__init__(xwidth, yheight, zmax=inf, rcurv=None)

Guide metric for position and direction.

Position is parametrized with following variables:

z: distance along guide, following curvature (if any).

t: transversal direction along mirrors, starting at (x+,y-) corner, towards (x+,y+) corner.

Direction is parametrized with following variables:

mu: cosine of angle between particle direction and mirror normal.

phi: azimuthal angle, starting from z direction, in [deg].

Axis system of reference can be changed by means of the ‘trasl’ and ‘rot’ arguments of the Geometry object.

Parameters

xwidth: float

Guide width.

yheight: float

Guide height.

zmax: float

Guide length.

rcurv: float

Curvature radius, defined as follows. Default is no curvature.

rcurv > 0 for curvature towards negative x

rcurv < 0 for curvature towards negative x

rcurv = 0 or rcurv = infinity for no curvature

transform(posdirs): Transform position and direction to guide variables.

inverse_transform(posdirs): Transform guide variables to position and direction.

save(mtree): Save Guide parameters into XML tree.

static load(mtree): Load parameters from XML tree and build Guide.

__module__ = 'kdsource.geom'

class kdsource.geom.Isotrop(keep_xdir=False, keep_ydir=False, keep_zdir=False)

Bases: Metric

__init__(keep_xdir=False, keep_ydir=False, keep_zdir=False)

Simple metric for direction, with no transformation.

Distance is measured as the euclidean distance between 3D unitary direction vectors.

Axis orientation can be changed by means of the rot’ argument of the Geometry object.

Parameters

keep_xdir: bool: If True, when using the source for sampling new particles, perturbation will not change dirx sign.
keep_ydir: bool: If True, when using the source for sampling new particles, perturbation will not change diry sign.
keep_zdir: bool: If True, when using the source for sampling new particles, perturbation will not change dirz sign.

mean(dirs=None, vecs=None, weights=None)

Mean of directions.

Mean is computed as the euclidean mean of direction vectors, normalized to 1.

Parameters

dirs: array-like, optional: Array of directions.
vecs: array-like, optional: Array of parametrized directions. If set, overrides dirs.
weights: array-like, optional: Array of particle statistic weights.

std(dirs=None, vecs=None, weights=None)

Standard deviation of directions.

Standard deviation is computed as the euclidean standard deviation of direction vectors minus its mean (computed with mean method).

Parameters

dirs: array-like, optional: Array of directions.
vecs: array-like, optional: Array of parametrized directions. If set, overrides dirs.
weights: array-like, optional: Array of particle statistic weights.

save(mtree): Save Guide parameters into XML tree.

static load(mtree): Load parameters from XML tree and build Isotrop.

__module__ = 'kdsource.geom'

class kdsource.geom.Polar

Bases: Metric

__init__()

Polar parametrization for direction.

Polar angles are defined as follows:: theta: angle between direction and z, in [deg]. phi: azimuthal angle, starting from x direction, in [deg].

Axis orientation can be changed by means of the rot’ argument of the Geometry object.

transform(dirs): Transform directions to polar angles.

inverse_transform(tps): Transform polar angles to directions.

jac(dirs): Jacobian of polar transformation.

static load(mtree): Build Polar.

__module__ = 'kdsource.geom'

class kdsource.geom.PolarMu

Bases: Metric

__init__()

Polar parametrization for direction, with mu = cos(theta).

Polar parameters are defined as follows:: mu: cosine of angle between direction and z, in [deg]. phi: azimuthal angle, starting from x direction, in [deg].

Axis orientation can be changed by means of the rot’ argument of the Geometry object.

transform(dirs): Transform directions to polar parameters.

inverse_transform(tps): Transform polar parameters to directions.

static load(mtree): Build PolarMu.

__module__ = 'kdsource.geom'

kdsource.geom.GeomFlat(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, z=0, E0=10, keep_zdir=True, trasl=None, rot=None)

Build flat neutron source.

Energy metric is Lethargy, position metric is SurfXY and direction metric is Isotrop.

See Metric’s and Geometry constructors for parameters docs.

kdsource.geom.GeomFlatTemp(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, z=0, E0=10, keep_zdir=True, trasl=None, rot=None)

Build flat neutron source with time treatment.

Energy metric is Lethargy, position metric is SurfXY, direction metric is Isotrop, and time metric is Decade.

See Metric’s and Geometry constructors for parameters docs.

kdsource.geom.GeomGuide(xwidth, yheight, zmax=inf, rcurv=None, E0=10, trasl=None, rot=None)

Build neutron source for leaks thru guide mirrors.

Energy metric is Lethargy, position and direction metric is Guide.

See Metric’s and Geometry constructors for parameters docs.

kdsource.geom.GeomActiv(xmin=-inf, xmax=inf, ymin=-inf, ymax=inf, zmin=-inf, zmax=inf, trasl=None, rot=None)

Build photon volumetric activation source.

Energy metric is Energy, position metric is Vol and direction metric is Isotrop.

See Metric’s and Geometry constructors for parameters docs.

kdsource.kde module

Module for bandwidth selection methods

Eventually, all the content of this module should be merged into KDEpy.

kdsource.kde.bw_silv(dim, N_eff)

Silverman’s Rule.

Estimates optimal bandwidth based on dimension and effective number of samples.

Parameters

dim: int

Dimension of the KDE model.

N_eff: float

Effective number of samples, defined as:: N_eff = sum(weights)**2 / sum(weights**2)

Returns

bw_opt: float: Bandwidth resulting from Silverman’s Rule.

kdsource.kde.bw_knn(data, weights=None, K_eff=100, k=None, batch_size=10000)

K Nearest Neighbors method.

For each sample, computes its bandwidth as the distance to the k-th neighbor, within each batch.

Parameters

data: array-like: Array of samples. Must have shape (obs, dim).
weights: array-like, optional: Array of sample weights. By default all weights are 1.
K_eff: float, optional: Effective k for all dataset. It is used to compute the k so that the estimated number of neighbors in all dataset is K_eff.
k: int, optional: Number of neighbors within each batch. If set, overrides the value estimated with K_eff.
batch_size: int, optional: Batch size for KNN search.

Returns

bw_opt: numpy.array: Array of KNN bandwidths.

kdsource.kde.bw_mlcv(data, weights=None, n_splits=10, seed=None, grid=None, show=True)

Maximum Likelihood Cross-Validation bandwidth.

Builds a grid of bandwidths, and computes the cross-validation likelihood score for each. Chooses the bandwidth that maximizes the score, or raises exception if maximum score is found in beginning or end of bandwidth grid. Also plots the scores over the bandwidth grid.

Bandwidths grid is computed with seed and grid, so that:: bw_grid[i] = seed * grid[i]

Parameters

data: array-like: Array of samples. Must have shape (obs, dim).
weights: array-like, optional: Array of sample weights. By default all weights are 1.
n_splits: int, optional: Number of folds for cross-validation. Must be at least 2.
seed: float or array-like: Seed bandwidth for computing bandwidth grid. Can be constant or adaptive. It is recommended to use the output of other bandwidth selection method. By default it is computed with Silverman’s Rule.
grid: array-like: Grid of factors for computing bandwidth grid. Default: numpy.logspace(-1, 1, 20)
show: bool: Whether to show the scores plot.

Returns

bw_opt: float or numpy.array: Optimal bandwidth, with same shape as seed.

kdsource.kde.optimize_bw(bw_method, vecs, ws=None, weightfun=None, maskfun=None, **kwargs)

Optimize bandwidth with given method.

Parameters

bw_method: str: Bandwidth selection method. See bw_methods for available methods.
vecs: array-like: Array of samples. Must have shape (obs, dim).
ws: array-like, optional: Array of sample weights. By default all weights are 1.
weightfun: function, optional: Weighting function. If set, ws will be multiplied by weightfun(vecs).
maskfun: function, optional: Masking function. If set, vecs and ws will be replaced by vecs[maskfun(vecs)] and ws[maskfun(vecs)], respectively.
**kwargs: optional: Parameters to be passed to bandwidth selection method.

Returns

bw_opt: float or array-like: Output of the selected bandwidth selection method.

kdsource.stats module

Module for statistic analysis

kdsource.stats.apply_weight_mask(vecs, ws, weightfun=None, maskfun=None)

Apply weighting and masking functions to particle list.

Parameters

vecs: array-like: Array of particles. Must have shape (N, dim).
ws: array-like, optional: Array of particle weights.
weightfun: function, optional: Weighting function. If set, ws will be multiplied by weightfun(vecs).
maskfun: function, optional: Masking function. If set, vecs and ws will be replaced by vecs[maskfun(vecs)] and ws[maskfun(vecs)], respectively.

Returns

[vecs, ws]: list: Particle list and weights after applying weighting and masking functions.

kdsource.stats.convergence(vecs, ws, param, fracmin=0.1, steps=10, plot=True)

Compute statistical parameter over subsets of particle list.

A number of subsets of particle list are built, with linearly growing size and the last one being the full particle list. For each subset the specified statistic parameter is computed.

Parameters

vecs: array-like: Array of particles. Must have shape (N, dim).
ws: array-like, optional: Array of particle weights.
param: function: Function that computes the statistic parameter. Must have signature:

param(vecs, ws) -> [param, err]
fracmin: float, optional: Size of the first subset, as a fraction of the full particle list size.
steps: int, optional: Number of subsets to build and compute statistic parameter.
plot: bool, optional: Whether to show the plot of statistic parameter vs. subset size.

Returns

[Ns, params, errs]: list: Subsets sizes, statistic parameter values, and errors.

kdsource.stats.mean_weight(vecs, ws): Compute mean weight and its statistic error.

kdsource.stats.mean(vecs, ws, var): Compute mean particle and its statistic error.

kdsource.stats.std(vecs, ws, var): Compute particles standard deviation and its statistic error.

class kdsource.stats.Stats(vecs, ws, weightfun=None, maskfun=None)

Bases: object

__init__(vecs, ws, weightfun=None, maskfun=None)

Object for statistic analysis of particle list.

This class methods show the variation and convergence of some statistic parameters as the number of particles in the list grows. This can be used to determine whether the list size is enough or not, based on some user-defined criterion.

Weighting and masking functions allow using an importance function based on particle parameters, and selecting a region of the phase-space to analyze.

Parameters

vecs: array-like: Array of particles. Must have shape (N, dim).
ws: array-like, optional: Array of particle weights.
weightfun: function, optional: Weighting function. If set, ws will be multiplied by weightfun(vecs).
maskfun: function, optional: Masking function. If set, vecs and ws will be replaced by vecs[maskfun(vecs)] and ws[maskfun(vecs)], respectively.

mean_weight(fracmin=0.1, steps=10, plot=True)

Compute convergence of mean weight.

Parameters

fracmin: float, optional: Size of the first subset, as a fraction of the full particle list size.
steps: int, optional: Number of subsets to build and compute mean weight.
plot: bool, optional: Whether plot mean weight vs. subset size.

Returns

[Ns, params, errs]: list: Subsets sizes, mean weight values, and errors.

mean(var, varname=None, fracmin=0.1, steps=10, plot=True)

Compute convergence of a variable mean.

Parameters

var: int: Index of variable to compute mean.
varname: str, optional: Selected variable name for display in plot.
fracmin: float, optional: Size of the first subset, as a fraction of the full particle list size.
steps: int, optional: Number of subsets to build and compute variable mean.
plot: bool, optional: Whether plot variable mean vs. subset size.

Returns

[Ns, params, errs]: list: Subsets sizes, variable mean values, and errors.

std(var, varname=None, fracmin=0.1, steps=10, plot=True)

Compute convergence of a variable standard deviation.

Parameters

var: int: Index of variable to compute standard deviation.
varname: str, optional: Selected variable name for display in plot.
fracmin: float, optional: Size of the first subset, as a fraction of the full particle list size.
steps: int, optional: Number of subsets to build and compute variable standard deviation.
plot: bool, optional: Whether plot variable standard deviation vs. subset size.

Returns

[Ns, params, errs]: list: Subsets sizes, variable standard deviation values, and errors.

__dict__ = mappingproxy({'__module__': 'kdsource.stats', '__init__': <function Stats.__init__>, 'mean_weight': <function Stats.mean_weight>, 'mean': <function Stats.mean>, 'std': <function Stats.std>, '__dict__': <attribute '__dict__' of 'Stats' objects>, '__weakref__': <attribute '__weakref__' of 'Stats' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.stats'

__weakref__: list of weak references to the object (if defined)

kdsource.summary module

Module for collecting McStas and TRIPOLI-4 simulation outputs

kdsource.summary.read_bashoutput(bashoutput, mccode)

Process bash output of McStas or TRIPOLI-4 simulation.

Searches for number of produced particles (with KDSource) and simulation times.

Parameters

bashoutput: str: Name of file containing bash output of simulation.
mccode: “McStas” or “TRIPOLI”: MC code used for simulation.

class kdsource.summary.Summary(mccode, folder, bashoutput='bash.out', n_detectors=[], p_detectors=[], t4output=None, tallies=[])

Bases: object

__init__(mccode, folder, bashoutput='bash.out', n_detectors=[], p_detectors=[], t4output=None, tallies=[])

Object representing summary of MC simulation.

After a Monte Carlo simulation, this class helps collecting the results, such as recorded particle lists or tallies, and simulation times. It also computes integral magnitudes of said results.

Parameters

mccode: “McStas” or “TRIPOLI”: MC code used for simulation.
folder: str: Directory containing simulation output files.
bashoutput: str: Name of file containing bash output of simulation.
n_detectors: list, optional: Neutron lists recorded.
p_detectors: list, optional: Photon lists recorded.
t4outup: str, optional: Name of TRIPOLI output file. Ignore if mccode is not TRIPOLI.
tallies: list, optional: Tally names recorded in TRIPOLI. Ignore if mccode is not TRIPOLI.

compute(): Load results and compute integral magnitudes.

save(filename)

Save results in text file, able to be imported to spreadsheet.

Parameters

filename: str: Name of file to save results.

__dict__ = mappingproxy({'__module__': 'kdsource.summary', '__init__': <function Summary.__init__>, 'compute': <function Summary.compute>, 'save': <function Summary.save>, '__dict__': <attribute '__dict__' of 'Summary' objects>, '__weakref__': <attribute '__weakref__' of 'Summary' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.summary'

__weakref__: list of weak references to the object (if defined)

kdsource.tally module

Module for reading tallies from MC simulations

kdsource.tally.read_spectrum(spectrum=None)

Read and load decay spectrum from CSV file.

Decay spectrum files can be downloaded from:

https://www-nds.iaea.org/relnsd/vcharthtml/VChartHTML.html

Select nuclide, and in ‘Decay Radiation’ tab download Gamma table as CSV.

Parameters

spectrum: str: Name of CSV file with decay spectrum. Energy and intensity values must be in first and third column, respectively. If None, empty spectrum is returned.

Returns

[Es, ws]: list: Energy and intensity values.

class kdsource.tally.T4Tally(outputfile, tallyname, spectrum=None, geomplot=None, J=1.0)

Bases: object

varnames = ['x', 'y', 'z']

varmap = {'x': 0, 'y': 1, 'z': 2}

units = ['cm', 'cm', 'cm']

__init__(outputfile, tallyname, spectrum=None, geomplot=None, J=1.0)

Object for reading TRIPOLI-4 3D tallies.

Tally must have mesh of type EXTENDED_MESH, and FRAME CARTESIAN.

This object has two main uses:

Reading activation tallies and converting them to gamma lists, which can then be used to generate KDSource object.
Reading and plotting dose maps.

Parameters

outputfile: str: Name of output file generated by TRIPOLI-4.
talyname: str: Name of the tally to read.
spectrum: str, optional: Name of file with decay spectrum. Needed only for converting activation tally to gamma list.
geomplot: str, optional: Name of file with geometry graph, to add contour lines to plot. This file must be generated with GRAPH command, on the same region than the tally.
J: float, optional: Intensity of the source of the TRIPOLI-4 simulation, in [1/s].

save_tracks(filename=None)

Save activation tally as gamma list.

For each tally cell, one gamma per decay energy is generated in its center, with random direction. Each gamma weight is the product of the cell tally value and the intensity of the energy value, normalized to have mean weight of 1 (or close).

Parameters

filename: str: Name of the particle list to generate. Will have MCPL format, but the .mcpl suffix is not needed.

Returns

trackfile: str, optional: Name of generated MCPL file. By default the tally name is used.

plot(var, cells=None, **kwargs)

1D plot of tally.

Parameters

var: str or int

Variable to be plotted, or its index. Names and indices of variables can be found in cls.varnames.

cells: list or None

Indices of fixed cells for other variables. If None, tally is averaged over other variables.

**kwargs: optional

Additional parameters for plotting options:

xscale: ‘linear’ or ‘log’
Scale for x axis. Default: ‘linear’
yscale: ‘linear’ or ‘log’. Default: ‘log’
Scale for y axis.
fact: float
Factor to apply on all densities. Default: 1
label: string
Label for plot legend.

Returns

[fig, [scores, errs]]:: Figure object, and plotted tally values and statistic errors.

plot2D(vrs, cell=None, geomplot=False, levelcurves=None, **kwargs)

2D plot of tally.

Parameters

vrs: list

Variables to be plotted, or its indices. Names and indices of variables can be found in cls.varnames.

cell: int or None

Index of fixed cell for other variable. If None, tally is averaged over other variable.

geomplot: bool

Whether to plot geometry contours.

levelcurves: list or None

If a list, gives the tally values to plot level curves.

**kwargs: optional

Additional parameters for plotting options:

scale: ‘linear’ or ‘log’: Scale for color map. Default: ‘log’
fact: float: Factor to apply on all densities. Default: 1

Returns

[fig, [scores, errs]]:: Figure object, and plotted tally values and statistic errors.

__dict__ = mappingproxy({'__module__': 'kdsource.tally', 'varnames': ['x', 'y', 'z'], 'varmap': {'x': 0, 'y': 1, 'z': 2}, 'units': ['cm', 'cm', 'cm'], '__init__': <function T4Tally.__init__>, 'save_tracks': <function T4Tally.save_tracks>, 'plot': <function T4Tally.plot>, 'plot2D': <function T4Tally.plot2D>, '__dict__': <attribute '__dict__' of 'T4Tally' objects>, '__weakref__': <attribute '__weakref__' of 'T4Tally' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.tally'

__weakref__: list of weak references to the object (if defined)

kdsource.utils module

Module for utility functions

kdsource.utils.pt2pdg(pt)

Convert particle type to pdgcode.

Valid particle types are:

‘n’: neutron (pdgcode = 2112)

‘p’: photon (pdgcode = 22)

‘e’: electron (pdgcode = 11)

For other particle types, 0 is returned.

Parameters

pt: str: Particle type.

kdsource.utils.pdg2pt(pdgcode)

Convert pdgcode to particle type.

Available pdgcode to be converted are:

2112: neutron (pt = ‘n’)

22: photon (pt = ‘p’)

11: electron (pt = ‘e’)

For other pdgcodes, ‘0’ is returned.

Parameters

pdgcode: int: Particle PDG code.

kdsource.utils.H10(pt='n', ref='ICRP')

Object for loading and interpolating dosimetric factors.

Units are MeV for energy and pSv*cm^2 for dosimetric factors.

Parameters

pt: str: Particle type. Available are ‘n’ (neutron) and ‘p’ (photon).
ref: str: Reference for dosimetric factors. Available are: - ‘ICRP’: International Commission on Radiological Protection. - ‘ARN’: Nuclear Regulatory Authority of Argentina.

class kdsource.utils.Box(vec0, vec1)

Bases: object

__init__(vec0, vec1)

Functor class for box-style masks.

Parameters

vec0: array-like: Lower limit for each variable.
vec0: array-like: Upper limit for each variable.

__call__(vecs)

Box-style masking function.

Parameters

vecs: array-like: Array of points to evaluate mask.

Returns

mask: array-like

Array of bools with same len as vecs. Each element has the following value:

True: if the corresponding vecs element has all its
variables above vec0 values and below vec1 values.
False: otherwise.

__dict__ = mappingproxy({'__module__': 'kdsource.utils', '__init__': <function Box.__init__>, '__call__': <function Box.__call__>, '__dict__': <attribute '__dict__' of 'Box' objects>, '__weakref__': <attribute '__weakref__' of 'Box' objects>, '__doc__': None, '__annotations__': {}})

__module__ = 'kdsource.utils'

__weakref__: list of weak references to the object (if defined)