Nothing
massProps
In massProps: Calculate Mass Properties and Uncertainties of Tree Structures
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
\newcommand{\sumwi}{\sum_{i=1}^{n}w_i}
\newcommand{\sumwivi}[1]{\sum_{i=1}^{n}{w_i}{{#1}_i}}
\newcommand{\Ixx}{I_{XX}}
\newcommand{\Iyy}{I_{YY}}
\newcommand{\Izz}{I_{ZZ}}
\newcommand{\Ixy}{I_{XY}}
\newcommand{\Ixz}{I_{XZ}}
\newcommand{\Iyz}{I_{YZ}}
\newcommand{\kx}{k_X}
\newcommand{\ky}{k_Y}
\newcommand{\kz}{k_Z}
\newcommand{\moia}[3]{
\sum_{i=1}^{n} \left[ {I_{#1}}_i
+ w_i \left( {#2}_i^2 + {#3}_i^2 \right)
- w_i \left( \bar{#2}^2 + \bar{#3}^2 \right)
\right]
}
\newcommand{\moib}[3]{
\sum_{i=1}^{n} \left{ {I_{#1}}_i
+ w_i \left[ \left( {#2}_i - \bar{#2} \right)^2 + \left( {#3}_i - \bar{#3} \right)^2 \right]
\right}
}
\newcommand{\di}[1]{({#1}_i - \bar{#1})}
\newcommand{\poia}[3]{\sum_{i=1}^{n} \left[ {I_{#1}}i + w_i {{#2}_i}{{#3}_i} -w_i (\bar{#2}\bar{#3})\right]}
\newcommand{\poib}[3]{\sum{i=1}^{n} \left[ {I_{#1}}i + w_i \di{#2}\di{#3}\right]}
\newcommand{\sigmacm}[1]{
\sqrt{ \sum{i=1}^n \left{ (w_i {{\sigma_\bar{#1}}}i)^2 + [{\sigma_w}_i \di{#1}]^2 \right} }
\bigg/ \sum{i=1}^{n}w_i \
}
\newcommand{\sigmamoi}[3]{\sqrt{ \sum_{i=1}^n \left{
\sigma_{{I_{#1}}i}^2 + \left[ 2 w_i \di{#2} \sigma{{#2}i} \right]^2 + \left[ 2 w_i \di{#3} \sigma{{#3}i} \right]^2 + \left[ \left(\di{#2}^2 + \di{#3}^2 \right)\sigma{w_i}\right]^2 \right} }
}
\newcommand{\sigmapoi}[3]{\sqrt{ \sum_{i=1}^n \left{
\sigma_{{I_{#1}}i}^2 + \left[ \di{#2} w_i \sigma{{#3}i}\right]^2 + \left[ \di{#2}\di{#3}\sigma{w_i} \right]^2 + \left[ \di{#3} w_i \sigma_{{#2}_i} \right]^2 \right} }
}
Overview
The massProps package extends rollupTree with functions to recursively calculate mass properties (and optionally, their uncertainties) for arbitrary decomposition trees. Formulas implemented are described in a technical paper published by the Society of Allied Weight Engineers [@zimmerman:05:sawe].
Synopsis
Data Structures
massProps operates on two fundamental data structures: a mass properties table and a tree. The mass properties table has an entry for every item in a tree structure of items; the edges of the tree convey the parent-child relations among items. The two data structures are linked by the id column of the data frame, which must be a character vector of unique item identifiers, and the vertex names of the tree. The sets of identifiers must be identical.
Mass Property Table
Required Columns for Mass Properties
The Mass Property Table must contain the following columns. Other columns may exist and will remain unmodified.
-
id unique identifier for each item (row)
-
mass mass of the item (numeric)
-
Cx $x$-component of center of mass (numeric)
-
Cy $y$-component of center of mass (numeric)
-
Cx $z$-component of center of mass (numeric)
-
Ixx moment of inertia about the $x$ axis (numeric)
-
Iyy moment of inertia about the $y$ axis (numeric)
-
Izz moment of inertia about the $z$ axis (numeric)
-
Ixy product of inertia relative to the $x$ and $y$ axes (numeric)
-
Ixz product of inertia relative to the $x$ and $z$ axes (numeric)
-
Iyz product of inertia relative to the $y$ and $z$ axes (numeric)
-
POIconv either '+' or '-', indicating the sign convention for products of inertia. In the negative convention, for example, $I_{XY} \equiv -\int{xy \rho \, dV}$. In the positive convention, $I_{XY} \equiv \int{xy \rho \, dV}$.
-
Ipoint logical indicator that this item is considered a point mass. The same algebraic result can be achieved by setting all moments and products of inertia to zero, but rollup_mass_props() by default ensures that all leaf items in the tree have mass properties that correspond to physically-realizable objects. A zero inertia tensor will fail this check. Rather than relax the check (which is essential for trustworthy results), a TRUE value for Ipoint indicates that the inertia tensor should be excluded from computations.
Required Columns for Mass Properties Uncertainties
The following columns are required for uncertainty calculations:
-
sigma_mass mass uncertainty (numeric)
-
sigma_Cx $x$-component of center of mass uncertainty (numeric)
-
sigma_Cy $y$-component of center of mass uncertainty (numeric)
-
sigma_Cx $z$-component of center of mass uncertainty (numeric)
-
sigma_Ixx moment of inertia about the $x$ axis uncertainty (numeric)
-
sigma_Iyy moment of inertia about the $y$ axis uncertainty (numeric)
-
sigma_Izz moment of inertia about the $z$ axis uncertainty (numeric)
-
sigma_Ixy product of inertia relative to the $x$ and $y$ axes uncertainty (numeric)
-
sigma_Ixz product of inertia relative to the $x$ and $z$ axes uncertainty (numeric)
-
sigma_Iyz product of inertia relative to the $y$ and $z$ axes uncertainty (numeric)
It is the caller's responsibility to ensure that all values are expressed in appropriate and compatible units.
Tree
The tree is an igraph::graph with vertices named by identifiers in the mass properties table. It can be of arbitrary depth and shape as long as it satisfies certain well-formedness properties:
-
it is connected and acyclic (as an undirected graph), i.e., it is a tree
-
it is directed, with edge direction going from child to parent
-
it contains neither loops (self-edges) nor multiple edges
-
it contains a single root vertex (i.e., one whose out degree is zero)
Invocation
library(rollupTree)
library(massProps)
suppressPackageStartupMessages({library(igraph)})
Suppose we have the following mass properties table:
test_table
Suppose we also have this tree:
test_tree
plot(test_tree,layout=layout_as_tree(test_tree, 2, mode="in"), vertex.shape = 'none', edge.arrow.mode = 0)
Then we can compute mass properties for non-leaf elements by calling rollup_mass_props():
rollup_mass_props(test_tree, test_table)
Note that, although the table shows the parent of each element for clarity of exposition, the child-parent relations are coneveyed only by the tree passed as the first argument.
The input may also contain uncertainties data. This example is from the Society of Allied Weight Engineers:
na_mass_props_and_unc <- function(d, t, v) {
xyz <- c("x", "y", "z")
list(
mass = NA,
center_mass = c(x = NA, y = NA, z = NA),
inertia = matrix(nrow = 3, ncol = 3, dimnames = list(xyz, xyz)),
poi_conv = "+",
point = FALSE,
sigma_mass = NA,
sigma_center_mass = c(x = NA, y = NA, z = NA),
sigma_inertia = matrix(nrow = 3, ncol = 3, dimnames = list(xyz, xyz))
)
}
na_mass_props_and_unc_update <- function(d, t, s) {
update_mass_props_and_unc(d, t, s, override = na_mass_props_and_unc)
}
sawe_input <- rollup(sawe_tree, sawe_table, update = na_mass_props_and_unc_update, validate_ds = validate_mass_props_and_unc_table)
sawe_input
rollup_mass_props_and_unc(sawe_tree, sawe_input)
Objectives and Strategy
The objective of this package is to provide a trustworthy, well-documented, reference implementation for computation of mass properties (and their uncertainties) of aggregate objects from those of their parts. Aggregation can be recursive (e.g., indentured Bill of Materials), so it must accommodate trees of arbitrary depth and shape.
Strategies for achieving the objective include
-
basing the calculations on published industry references,
-
re-casting those lengthy reference equations into concise vector or matrix forms to reduce the error surface for source code and exploit the capabilities of R, which treats vectors and matrices as first-class objects,
-
delegating orchestration to the rollupTree package, which, among other things, verifies that the input tree is well-formed and ensures proper ordering of computations,
-
ensuring that all asserted leaf mass properties and uncertainties correspond to physically-realizable objects,
-
coding in pure functional style, (i.e., avoiding mutable variables, implying iteration with Map() and Reduce()), and
-
covering the entire code base with unit tests.
The author has intentionally made no effort to micro-optimize for performance. In particular, the author is aware that representing the inertia and its uncertainty as 3 ⨉ 3 matrices is "inefficient" to the degree that it independently calculates values that are redundant by symmetry. "Inefficient", however, does not mean "slow". See [Performance Evaluation] below.
Theory
In this section, we state the reference equations [@zimmerman:05:sawe] and show, where applicable, how those equations can be rewritten in more concise form. The form of the equations actually implemented is displayed within a box, e.g. $\boxed{F = ma}$.
The reference uses the word weight and the symbol $w$ in equations. We interpret weight as mass. The reference refers to center of mass by its $x$, $y$, and $z$ components. Symbols for moments ($\Ixx$) and products ($\Ixy$) of inertia are conventional. Variables with $i$ subscripts designate properties of parts; those without designate properties of aggregates. The letter $\sigma$ denotes uncertainty. $\sigma_w$, for example, is the mass uncertainty.
Mass Properties
Mass
The mass equation is suitable as is.
$$
\boxed{
w = \sum_{i=1}^{n}w_i
}
$$
The corresponding R code is
amp$mass <- Reduce(`+`, Map(f = function(mp) mp$mass, mpl))
In this and the following code snippets, the variable mpl is a list of input mass property sets for parts, the variable mp is a formal parameter of an anonymous function applied to each member of mpl, and amp is the resulting aggregate mass property set. The line above is an R functional programming idiom for "set the mass value of the aggregate to the sum of the mass values of the parts".
Center of Mass
$$
\begin{align}
\bar{x} & = \sum_{i=1}^{n}{w_i}{{x}i} \bigg/ \sum{i=1}^{n}w_i\
\bar{y} & = \sum_{i=1}^{n}{w_i}{{y}i} \bigg/ \sum{i=1}^{n}w_i\
\bar{z} & = \sum_{i=1}^{n}{w_i}{{z}i} \bigg/ \sum{i=1}^{n}w_i\
\end{align}
$$
We can express center of mass as a 3-vector:
$$
\boxed{
\begin{align}
\boldsymbol{c}_i & = (x_i \quad y_i \quad z_i)^T \
\boldsymbol{\bar{c}}& = (\bar{x} \quad \bar{y} \quad \bar{z})^T
\end{align}
}
$$
Then
$$
\boxed{
\boldsymbol{\bar{c}} = \frac{1}{w} \sum_{i=1}^{n}{w_i}{{\boldsymbol{c}}_i}
}
$$
The corresponding R code is
amp$center_mass <- Reduce(`+`, Map(f = function(mp) mp$mass * mp$center_mass, mpl)) / amp$mass
Inertia Tensor
Moments of Inertia
$$
\begin{align}
\Ixx & = \moia{XX}{y}{z} & = \moib{XX}{y}{z} \
\Iyy & = \moia{YY}{x}{z} & = \moib{YY}{x}{z} \
\Izz & = \moia{ZZ}{x}{y} & = \moib{ZZ}{x}{y} \
\end{align}
$$
Products of Inertia
$$
\begin{align}
\Ixy & = \poia{XY}{x}{y} & = \poib{XY}{x}{y} \
\Ixz & = \poia{XZ}{x}{z} & = \poib{XZ}{x}{z} \
\Iyz & = \poia{YZ}{y}{z} & = \poib{YZ}{y}{z} \
\end{align}
$$
Matrix Formulation
Let $\boldsymbol{I}$ be the inertia tensor of the aggregate and $\boldsymbol{I}_i$ be that of part $i$. The equations for products of inertia above clearly follow the positive integral convention, so
$$
\boldsymbol{I} =
\left[ \begin{array}{rrr}
\Ixx & -\Ixy & -\Ixz \
-\Ixy & \Iyy & -\Iyz \
-\Ixz & -\Iyz & \Izz \
\end{array} \right]
$$
and similarly for $\boldsymbol{I}_i$.
Noting the repeated appearance of terms of the form $\di{x}\di{y}$, we form the outer product
$$
\boxed{
\begin{align}
\boldsymbol{d}_i & = (\di{x} \quad \di{y} \quad \di{z})^T \
\boldsymbol{Q}_i & = \boldsymbol{d}_i {\boldsymbol{d}_i}^T
\end{align}
}
$$ Then
$$
\begin{align}
\boldsymbol{Q}_i & =
\begin{bmatrix}
\di{x}^2 & \di{x}\di{y} & \di{x}\di{z} \
\di{y}\di{x} & \di{y}^2 & \di{y}\di{z} \
\di{z}\di{x} & \di{z}\di{y} & \di{z}^2 \
\end{bmatrix}
\end{align}
$$
Let $\boldsymbol{s}_i$ be the matrix of inertia tensor summands from the reference equations. That is,
$$
\boldsymbol{I} = \sum_{i=1}^{n} \boldsymbol{s}_i
$$
where $$
\begin{align}
\boldsymbol{s}_i & = \boldsymbol{I}_i \
& -
w_i \begin{bmatrix}
-\di{y}^2 - \di{z}^2 & \di{x}\di{y} & \di{x}\di{z} \
\di{x}\di{y} & -\di{x}^2 - \di{z}^2 & \di{y}\di{z} \
\di{x}\di{z} & \di{y}\di{z} & -\di{x}^2 - \di{y}^2 \
\end{bmatrix} \
& =
\boldsymbol{I}_i \
& -
w_i \begin{bmatrix}
\di{x}^2 & \di{x}\di{y} & \di{x}\di{z} \
\di{x}\di{y} & \di{y}^2 & \di{y}\di{z} \
\di{x}\di{z} & \di{y}\di{z} & \di{z}^2 \
\end{bmatrix} \
& -
w_i \begin{bmatrix}
-\di{x}^2 - \di{y}^2 - \di{z}^2 & 0 & 0 \
0 & -\di{x}^2 - \di{y}^2 - \di{z}^2 & 0 \
0 & 0 &-\di{x}^2 - \di{y}^2 - \di{z}^2 \
\end{bmatrix} \
& =
\boldsymbol{I}_i
-
w_i \left( \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i) \boldsymbol{1}_3 \right)
\end{align}
$$
where $\mathrm{tr}(\boldsymbol{Q}_i)$ is the trace of $\boldsymbol{Q}_i$, i.e., the sum of its diagonal elements, and $\boldsymbol{1}_3$ is the 3⨉3 identity matrix. Therefore
$$
\boxed{
\boldsymbol{I} = \sum_{i=1}^{n} \left( \boldsymbol{I}_i
-
w_i \boldsymbol{M}_i \right)
}
$$ where
$$
\boxed{
\boldsymbol{M}_i = \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i) \boldsymbol{1}_3
}
$$
The corresponding R code is
amp$inertia <- Reduce(`+`, Map(
f = function(mp) {
d <- amp$center_mass - mp$center_mass
Q <- outer(d, d)
M <- Q - sum(diag(Q)) * diag(3)
if (mp$point) -mp$mass * M else mp$inertia - mp$mass * M
},
mpl
))
Mass Property Uncertainties
Mass Uncertainty
The mass uncertainty equation is suitable as is.
$$
\boxed{
\sigma_w = \sqrt{ \sum_{i=1}^n {{\sigma_w}_i}^2 }
}
$$
The corresponding R code is
amp$sigma_mass = sqrt(Reduce(`+`, Map(f = function(v) v$sigma_mass^2, mpl)))
Center of Mass Uncertainty
$$
\begin{align}
\sigma_\bar{x} & = \sigmacm{x} \
\sigma_\bar{y} & = \sigmacm{y} \
\sigma_\bar{z} & = \sigmacm{z} \
\end{align}
$$
As before, we create a 3-vector for center of mass uncertainties. Let
$$
\boxed{
\begin{align}
\boldsymbol{\sigma_c} & = (\sigma_\bar{x} \quad \sigma_\bar{y} \quad \sigma_\bar{z})^T \
{\boldsymbol{\sigma_c}}i & = ({\sigma\bar{x}}i \quad {\sigma\bar{y}}i \quad {\sigma\bar{z}}_i)^T
\end{align}
}
$$
If we construe (as R does) squaring and taking square roots of vectors element-wise, then
$$
\boxed{
\boldsymbol{\sigma_c} = \frac{1}{w}
\sqrt{ \sum_{i=1}^n \left{ (w_i {{\boldsymbol{\sigma_c}}}_i)^2 + [{\sigma_w}_i ({\boldsymbol{c}}_i - \bar{\boldsymbol{c}})]^2 \right}}
}
$$
The corresponding R code is
amp$sigma_center_mass = sqrt(Reduce(`+`, Map(
f = function(v) {
(v$mass * v$sigma_center_mass)^2 +
(v$sigma_mass * (v$center_mass - amp$center_mass))^2
},
mpl
))) / amp$mass
Inertia Tensor Uncertainty
Moments of Inertia Uncertainties
$$
\begin{align}
\sigma_{\Ixx} & = \sigmamoi{XX}{y}{z} \
\sigma_{\Iyy} & = \sigmamoi{YY}{x}{z} \
\sigma_{\Izz} & = \sigmamoi{ZZ}{x}{y} \
\end{align}
$$
Products of Inertia Uncertainties
$$
\begin{align}
\sigma_{\Ixy} & = \sigmapoi{XY}{x}{y} \
\sigma_{\Ixz} & = \sigmapoi{XZ}{x}{z} \
\sigma_{\Iyz} & = \sigmapoi{YZ}{y}{z} \
\end{align}
$$
Matrix Formulation
Let
$$
\boxed{
\begin{align}
\boldsymbol{d}i & = (\di{x} \quad \di{y} \quad \di{z})^T \
{\boldsymbol{\sigma_c}}_i & = ({\sigma\bar{x}}i \quad {\sigma\bar{y}}i \quad {\sigma\bar{z}}_i)^T \
\boldsymbol{P}_i & = \boldsymbol{d}_i {\boldsymbol{\sigma_c}}_i^T \
\boldsymbol{Q}_i & = \boldsymbol{d}_i {\boldsymbol{d}_i}^T
\end{align}
}
$$
Then
$$
\begin{align}
\boldsymbol{P}i & =
\begin{bmatrix}
\di{x}\sigma{x_i} &\di{x}\sigma_{y_i} &\di{x}\sigma_{z_i} \
\di{y}\sigma_{x_i} & \di{y}\sigma_{y_i} & \di{y}\sigma_{z_i} \
\di{z}\sigma_{x_i} & \di{z}\sigma_{y_i} & \di{z}\sigma_{z_i} \
\end{bmatrix}
\
\
\boldsymbol{Q}_i & =
\begin{bmatrix}
\di{x}^2 &\di{x}\di{y} &\di{x}\di{z} \
\di{y}\di{x} & \di{y}^2 & \di{y}\di{z} \
\di{z}\di{x} & \di{z}\di{y} & \di{z}^2 \
\end{bmatrix}
\end{align}
$$
Let $\boldsymbol{s}_i^2$ be the matrix of inertia tensor uncertainty summands in the standard formulas for a given subcomponent $i$ above. That is,
$$
\boldsymbol{\sigma_I}^2 = \sum_{i=1}^{n} \boldsymbol{s}_i^2
$$
Let ${p_X}i$, ${p_Y}_i$, and ${p_Z}_i$ be the respective diagonal elements of $P_i$. Let $\boldsymbol{1}_3$ be the 3 ⨉ 3 identity matrix. If we interpret squaring a matrix as the Hadamard (element-wise) product with itself, then $$
\begin{align}
\boldsymbol{s}_i^2 & = {\boldsymbol{\sigma_I}}_i^2 \
& +
\begin{bmatrix}
2 w_i \di{y} \sigma{y_i} & w_i\di{x} \sigma_{y_i} & w_i\di{x} \sigma_{z_i} \
w_i\di{x} \sigma_{y_i} & 2 w_i\di{x} \sigma_{x_i} & w_i \di{y} \sigma_{z_i} \
w_i\di{x} \sigma_{z_i} & w_i \di{y} \sigma_{z_i} & 2 w_i\di{x} \sigma_{x_i}
\end{bmatrix}^2 \
& +
\begin{bmatrix}
2 w_i \di{z} \sigma_{z_i} & w_i \di{y} \sigma_{x_i} & w_i \di{z} \sigma_{x_i} \
w_i \di{y} \sigma_{x_i} & 2 w_i \di{z} \sigma_{z_i} & w_i \di{z} \sigma_{y_i} \
w_i \di{z} \sigma_{x_i} & w_i \di{z} \sigma_{y_i} & 2 w_i \di{y} \sigma_{y_i}
\end{bmatrix}^2 \
& +
\begin{bmatrix}
(\di{y}^2 + \di{z}^2)\sigma_{w_i} &\di{x}\di{y}\sigma_{w_i} &\di{x}\di{z}\sigma_{w_i} \
\di{y}\di{x}\sigma_{w_i} & (\di{x}^2 + \di{z}^2)\sigma_{w_i} & \di{y}\di{z}\sigma_{w_i} \
\di{z}\di{x}\sigma_{w_i} & \di{z}\di{y}\sigma_{w_i} & (\di{x}^2 + \di{y}^2)\sigma_{w_i} \
\end{bmatrix}^2 \ \
& = {\boldsymbol{\sigma_I}}i^2 \
& + w_i^2 \left(
\boldsymbol{P}_i -
\begin{bmatrix}
\di{x}\sigma{x_i} - 2 \di{y}\sigma_{y_i} & 0 & 0 \
0 & \di{y}\sigma_{y_i} - 2\di{x}\sigma_{x_i} & 0 \
0 & 0 & \di{z}\sigma_{y_i} - 2\di{x}\sigma_{x_i} \
\end{bmatrix}
\right) ^2 \
& + w_i^2 \left(
\boldsymbol{P}i^T -
\begin{bmatrix}
\di{x}\sigma{x_i} - 2 \di{z}\sigma_{y_i} & 0 & 0 \
0 & \di{y}\sigma_{y_i} - 2 \di{z}\sigma_{z_i} & 0 \
0 & 0 & \di{z}\sigma_{y_i} - 2 \di{y}\sigma_{y_i} \
\end{bmatrix}
\right) ^2 \
& + \sigma_{w_i}^2 \left(
\boldsymbol{Q}i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3
\right)^2 \ \
& = {\boldsymbol{\sigma_I}}_i^2 \
& + w_i^2 \left(
\boldsymbol{P}_i -
\begin{bmatrix}
{p_X}_i - 2 {p_Y}_i & 0 & 0 \
0 & {p_Y}_i - 2 {p_X}_i & 0 \
0 & 0 & {p_Z}_i - 2 {p_X}_i \
\end{bmatrix}
\right) ^2 \
& + w_i^2 \left(
\boldsymbol{P}_i^T -
\begin{bmatrix}
{p_X}_i - 2{p_Z}_i & 0 & 0 \
0 & {p_Y}_i - 2 {p_Z}_i & 0 \
0 & 0 & {p_Z}_i - 2 {p_Y}_i \
\end{bmatrix}
\right) ^2 \
& + \sigma{w_i}^2 \left(
\boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3
\right)^2
\end{align}
$$
Finally,
$$
\boxed{
\boldsymbol{\sigma_I} = \sqrt{
\sum_{i=1}^{n} {\left{
{\boldsymbol{\sigma_I}}_i^2
+{\boldsymbol{M}_4}_i
\right}}
}
}
$$
where
$$
\boxed{
\begin{align}
{\boldsymbol{M}1}_i & = \boldsymbol{P}_i -
\begin{bmatrix}
{p_X}_i - 2{p_Y}_i & 0 & 0 \
0 & {p_Y}_i - 2{p_X}_i & 0 \
0 & 0 & 2{p_Z}_i - 2{p_X}_i \
\end{bmatrix} \
{\boldsymbol{M}_2}_i & = \boldsymbol{P}_i^T -
\begin{bmatrix}
{p_X}_i - 2{p_Z}_i & 0 & 0 \
0 & {p_Y}_i - 2{p_Z}_i & 0 \
0 & 0 & 2{p_Z}_i - 2{p_Y}_i \
\end{bmatrix} \
{\boldsymbol{M}_3}_i & = \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3 \
{\boldsymbol{M}_4}_i & = w_i^2 \left(
{{\boldsymbol{M}_1}_i}^2 + {{\boldsymbol{M}_2}_i}^2
\right)
+ \left( \sigma{w_i} {\boldsymbol{M}_3}_i
\right)^2
\end{align}
}
$$
The corresponding R code is
amp$sigma_inertia = sqrt(Reduce(`+`, Map(
f = function(v) {
d <- v$center_mass - amp$center_mass
P <- outer(d, v$sigma_center_mass)
p <- diag(P)
Q <- outer(d, d)
M1 <- P - diag(p - 2 * p[c("y", "x", "x")])
M2 <- t(P) - diag(p - 2 * p[c("z", "z", "y")])
M3 <- Q - sum(diag(Q)) * diag(3)
M4 <- v$mass^2 * (M1^2 + M2^2) + (v$sigma_mass * M3)^2
if (v$point) M4 else v$sigma_inertia^2 + M4
},
mpl
)))
Radii of Gyration and Their Uncertainties
By definition:
$$
\begin{align}
\kx & = \sqrt{\Ixx / w} \
\ky & = \sqrt{\Iyy / w} \
\kz & = \sqrt{\Izz / w} \
\end{align}
$$ Let
$$
\boxed{
\begin{align}
\boldsymbol{k} & = (\kx \quad \ky \quad \kz)^T \
\boldsymbol{I} & = (\Ixx \quad \Iyy \quad \Izz)^T
\end{align}
}
$$ Then
$$
\boxed{
\boldsymbol{k} = \sqrt{ \boldsymbol{I} / w }
}
$$
The corresponding R code is
rg <- get_mass_props(d, i)
rg$radii_gyration <- sqrt(diag(rg$inertia) / rg$mass)
The SAWE reference gives equations for uncertainties of radii of gyration in recursive form, but as these radii are simply functions of moments of inertia and mass, we should be able to express their uncertainties in terms of uncertainties of moments of inertia and mass by applying standard uncertainty propagation theory [@wikipedia:uncprop].
Let
$$
\boxed{
\begin{align}
\sigma_\boldsymbol{k} & = (\sigma_{\kx} \quad \sigma_{\ky} \quad \sigma_{\kz})^T \
\sigma_\boldsymbol{I} & = (\sigma_{\Ixx} \quad \sigma_{\Iyy} \quad \sigma_{\Izz})^T
\end{align}
}
$$
Then
$$
\begin{align}
{\sigma_\boldsymbol{k}}^2 & \approx \left( \frac{\partial \boldsymbol{k}}{\partial \boldsymbol{I}} \right)^2 {\sigma_\boldsymbol{I}}^2 + \left( \frac{\partial \boldsymbol{k}}{\partial w} \right)^2 {\sigma_w}^2 + 2 \frac{\partial \boldsymbol{k}}{\partial \boldsymbol{I}} \frac{\partial \boldsymbol{k}}{\partial w} \sigma_{\boldsymbol{I}w} \
& \approx \left( \frac{1}{2 \sqrt{w \boldsymbol{I}}} \right)^2 {\sigma_\boldsymbol{I}}^2 + \left(\frac{- \sqrt{\boldsymbol{I}}}{2 w^{3/2}}\right)^2 {\sigma_w}^2 + 2 \left( \frac{1}{2 \sqrt{w \boldsymbol{I}}} \right) \left(\frac{- \sqrt{\boldsymbol{I}}}{2 w^{3/2}}\right) \sigma_{\boldsymbol{I}w} \
& \approx \frac{1}{4 w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{4 w^3} {\sigma_w}^2 - \frac{1}{2 w^2} \sigma_{\boldsymbol{I}w} \
\end{align}
$$
where $\sigma_{\boldsymbol{I}w}$ is the covariance between $\boldsymbol{I}$ and $w$.
Let $W$ and $X$ be random variables such that
$$
\begin{align}
W & = \sum_{i=1}^{n} (w_i + {\epsilon_w}i) \
X & = \sum{i=1}^{n} \left{ {\Ixx}i + {\epsilon\Ixx}i
+ (w_i + {\epsilon_w}_i) \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right]
\right}
\
\end{align}
$$
where $E[{\epsilon_w}_i] = 0$ and $E[{\epsilon\Ixx}i] = 0$ for all $i$,
$E[{\epsilon_w}_i {\epsilon_w}_j] = 0$ and $E[{\epsilon\Ixx}i {\epsilon\Ixx}j] = 0$ for all $i \neq j$, and
$E[{\epsilon_w}_i {\epsilon\Ixx}_j] = 0$ for all $i$ and $j$. It is clear from the linearity of $W$ and $X$ in $w_i$ and $\Ixx_i$ that
$$
E[W] = E \left[\sum_{i=1}^{n} (w_i + {\epsilon_w}i) \right] = w + \sum{i=1}^{n} {\epsilon_w}i = w
$$
and
$$
\begin{align}
E[X] & = E \left[ \sum{i=1}^{n} \left{ {\Ixx}i + {\epsilon{\Ixx}}i
+ (w_i + {\epsilon_w}_i) \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right] \
& = \Ixx + E \left[ \sum{i=1}^{n} \left{ {\epsilon_{\Ixx}}_i
+ {\epsilon_w}_i \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right] \
& = \Ixx
\end{align}
$$
Therefore,
$$
\begin{align}
\sigma_{\Ixx w} & = E \left[ \left( W - E[W] \right) \left( X - E[X] \right)\right] \
& = E \left[ \left( \sum_{i=1}^{n} {\epsilon_w}i \right)
\left( \sum{i=1}^{n} \left{ {\epsilon_{\Ixx}}i
+ {\epsilon_w}_i \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right) \right] \
& = E \left[ \sum{i=1}^{n} {\epsilon_w}^2_i \left( \left( {y}i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \right] \
& = \sum{i=1}^{n} E \left[ {\epsilon_w}^2_i \right] \left( \left( {y}i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \
& = \sum{i=1}^{n} {\sigma_w}^2_i \left( \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \
\end{align}
$$
Exploiting symmetry, let
$$
\boxed{
\begin{align}
\boldsymbol{d}_i & = (({x}_i - \bar{x}) \quad ({y}_i - \bar{y}) \quad ({z}_i - \bar{z}))^T \
\boldsymbol{s}_i & = (\boldsymbol{d}_i^T\boldsymbol{d}_i \quad \boldsymbol{d}_i^T\boldsymbol{d}_i \quad \boldsymbol{d}_i^T\boldsymbol{d}_i)^T \
\end{align}
}
$$
$$
\begin{align}
{\sigma_\boldsymbol{k}}^2 & \approx \frac{1}{4 w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{4 w^3} {\sigma_w}^2 -
\frac{1}{2 w^2} \sum_{i=1}^{n} {\sigma_w}^2_i (\boldsymbol{s}_i - \boldsymbol{d}_i^2 ) \
\end{align}
$$
Therefore, we define
$$
\boxed{
\begin{align}
{\sigma_\boldsymbol{k}} & = \frac{1}{2} \sqrt{ \frac{1}{w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{w^3} {\sigma_w}^2
- \frac{2}{w^2} \sum_{i=1}^{n} {\sigma_w}^2_i (\boldsymbol{s}_i - \boldsymbol{d}_i^2 )
}
\end{align}
}
$$
The corresponding R code is
amp <- get_mass_props_and_unc(ds, target)
I <- diag(amp$inertia)
sigma_I <- diag(amp$sigma_inertia)
amp$sigma_radii_gyration <- sqrt(
sigma_I^2 / (amp$mass * I) + (I * amp$sigma_mass^2) / amp$mass^3 -
2 / amp$mass^2 * Reduce(
`+`,
Map(
f = function(s) {
mp <- get_mass_props_and_unc(ds, s)
d2 <- (mp$center_mass - amp$center_mass)^2
mp$sigma_mass^2 * (sum(d2) - d2)
},
sources
),
init <- c(0, 0, 0)
)
) / 2
```
# Testing and Validation
## Comparison With Independently-Calculated Results
In this section we will calculate the results for the SAWE example step by step and compare them with the package results. The inputs are:
```r
sawe_input <- sawe_table[which(sawe_table$id != "Combined"), ]
sawe_input
Our computed result is
t <- rollup_radii_of_gyration_unc(sawe_tree,
add_radii_of_gyration(
rollup_mass_props_and_unc(sawe_tree, sawe_table)
)
)
sawe_result <- t[t$id == "Combined", ]
sawe_result
Mass
agree <- function(l) if (l) "agrees" else stop("disagreement")
mass <- sum(sawe_input$mass)
The independently-calculated mass is
mass
This r agree(isTRUE(all.equal(mass, sawe_result$mass))) with the computed result.
Center of Mass
C <- apply(sawe_input$mass / mass * sawe_input[, c("Cx", "Cy", "Cz")], 2, sum)
The independently-calculated center of mass is
C
This r agree(isTRUE(all.equal(unname(C), c(sawe_result$Cx, sawe_result$Cy, sawe_result$Cz)))) with the computed result.
Moments of Inertia
moi <- function(I, v1, v2, m, c1, c2) {
sum(I + m * ((v1^2 + v2^2) - (c1^2 + c2^2)))
}
MOI <- c(
Ixx = moi(sawe_input$Ixx, sawe_input$Cy, sawe_input$Cz, sawe_input$mass, C["Cy"], C["Cz"]),
Iyy = moi(sawe_input$Iyy, sawe_input$Cx, sawe_input$Cz, sawe_input$mass, C["Cx"], C["Cz"]),
Izz = moi(sawe_input$Izz, sawe_input$Cx, sawe_input$Cy, sawe_input$mass, C["Cx"], C["Cy"])
)
The independently-calculated moments of inertia are
MOI
This r agree(isTRUE(all.equal(MOI, c(Ixx = sawe_result$Ixx, Iyy = sawe_result$Iyy, Izz = sawe_result$Izz)))) with the computed result.
Products of Inertia
poi <- function(I, v1, v2, m, c1, c2) {
sum(I + m * (v1 * v2 - c1 * c2))
}
POI <- c(
Ixy = poi(sawe_input$Ixy, sawe_input$Cx, sawe_input$Cy, sawe_input$mass, C["Cx"], C["Cy"]),
Ixz = poi(sawe_input$Ixz, sawe_input$Cx, sawe_input$Cz, sawe_input$mass, C["Cx"], C["Cz"]),
Iyz = poi(sawe_input$Iyz, sawe_input$Cy, sawe_input$Cz, sawe_input$mass, C["Cy"], C["Cz"])
)
The independently-calculated products of inertia are
POI
This r agree(isTRUE(all.equal(POI, c(Ixy = sawe_result$Ixy, Ixz = sawe_result$Ixz, Iyz = sawe_result$Iyz)))) with the computed result.
Radii of Gyration
rog <- function(I, m) sqrt(I / m)
ROG <- c(
kx = rog(sawe_result$Ixx, sawe_result$mass),
ky = rog(sawe_result$Iyy, sawe_result$mass),
kz = rog(sawe_result$Izz, sawe_result$mass)
)
The independently-calculated radii of gyration are
ROG
This r agree(isTRUE(all.equal(ROG, c(kx = sawe_result$kx, ky = sawe_result$ky, kz = sawe_result$kz)))) with the computed result.
Mass Uncertainty
sigma_mass <- sqrt(sum(sawe_input$sigma_mass^2))
The independently-calculated mass uncertainty is
sigma_mass
This r agree(isTRUE(all.equal(sigma_mass, sawe_result$sigma_mass))) with the computed result.
Center of Mass Uncertainty
sigma_cm <- function(m, sigma_v, sigma_m, v, c, mass) {
sqrt(sum((m * sigma_v)^2 + (sigma_m * (v - c))^2)) / mass
}
sigma_C <- c(
sigma_Cx = sigma_cm(sawe_input$mass, sawe_input$sigma_Cx, sawe_input$sigma_mass, sawe_input$Cx, C["Cx"], mass),
sigma_Cy = sigma_cm(sawe_input$mass, sawe_input$sigma_Cy, sawe_input$sigma_mass, sawe_input$Cy, C["Cy"], mass),
sigma_Cz = sigma_cm(sawe_input$mass, sawe_input$sigma_Cz, sawe_input$sigma_mass, sawe_input$Cz, C["Cz"], mass)
)
The independently-calculated center of mass uncertainties are
sigma_C
This r agree(isTRUE(all.equal(sigma_C, c(sigma_Cx = sawe_result$sigma_Cx, sigma_Cy = sawe_result$sigma_Cy, sigma_Cz = sawe_result$sigma_Cz)))) with the computed result.
Moments of Inertia Uncertainties
sigma_moi <- function(sigma_I, mass, sigma_mass, v1, v2, c1, c2, sigma_v1, sigma_v2) {
sqrt(sum(
sigma_I^2 +
(2 * mass * (v1 - c1) * sigma_v1)^2 +
(2 * mass * (v2 - c2) * sigma_v2)^2 +
(((v1 - c1)^2 + (v2 - c2)^2) * sigma_mass)^2
))
}
sigma_MOI <- c(
sigma_Ixx = sigma_moi(sawe_input$sigma_Ixx, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cy,
sawe_input$Cz, C["Cy"], C["Cz"], sawe_input$sigma_Cy, sawe_input$sigma_Cz),
sigma_Iyy = sigma_moi(sawe_input$sigma_Iyy, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx,
sawe_input$Cz, C["Cx"], C["Cz"], sawe_input$sigma_Cx, sawe_input$sigma_Cz),
sigma_Izz = sigma_moi(sawe_input$sigma_Izz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx,
sawe_input$Cy, C["Cx"], C["Cy"], sawe_input$sigma_Cx, sawe_input$sigma_Cy)
)
The independently-calculated moments of inertia uncertainties are
sigma_MOI
This r agree(isTRUE(all.equal(sigma_MOI, c(sigma_Ixx = sawe_result$sigma_Ixx, sigma_Iyy = sawe_result$sigma_Iyy, sigma_Izz = sawe_result$sigma_Izz)))) with the computed result.
Products of Inertia Uncertainties
sigma_poi <- function(sigma_I, mass, sigma_mass, v1, v2, c1, c2, sigma_v1, sigma_v2) {
sqrt(sum(
sigma_I^2 +
((v1 - c1) * mass * sigma_v2)^2 +
((v1 - c1) * (v2 - c2) * sigma_mass)^2 +
((v2 - c2) * mass * sigma_v1)^2
))
}
sigma_POI <- c(
sigma_Ixy = sigma_poi(sawe_input$sigma_Ixy, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx,
sawe_input$Cy, C["Cx"], C["Cy"], sawe_input$sigma_Cx, sawe_input$sigma_Cy),
sigma_Ixz = sigma_poi(sawe_input$sigma_Ixz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx,
sawe_input$Cz, C["Cx"], C["Cz"], sawe_input$sigma_Cx, sawe_input$sigma_Cz),
sigma_Iyz = sigma_poi(sawe_input$sigma_Iyz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cy,
sawe_input$Cz, C["Cy"], C["Cz"], sawe_input$sigma_Cy, sawe_input$sigma_Cz)
)
The independently-calculated products of inertia uncertainties are
sigma_POI
This r agree(isTRUE(all.equal(sigma_POI, c(sigma_Ixy = sawe_result$sigma_Ixy, sigma_Ixz = sawe_result$sigma_Ixz, sigma_Iyz = sawe_result$sigma_Iyz)))) with the computed result.
Radii of Gyration Uncertainties
sigma_rog <- function(mt, It, sigma_mt, sigma_It, v1, c1, v2, c2, sigma_m) {
sqrt(sigma_It^2 / (4 * mt * It) + (It * sigma_mt^2) / (4 * mt^3) - sum(sigma_m^2 * ((v1 - c1)^2 + (v2 - c2)^2)) / (2 * mt^2))
}
sigma_ROG <- c(
sigma_kx = sigma_rog(sawe_result$mass, sawe_result$Ixx, sawe_result$sigma_mass, sawe_result$sigma_Ixx,
sawe_input$Cy, sawe_result$Cy, sawe_input$Cz, sawe_result$Cz, sawe_input$sigma_mass),
sigma_ky = sigma_rog(sawe_result$mass, sawe_result$Iyy, sawe_result$sigma_mass, sawe_result$sigma_Iyy,
sawe_input$Cx, sawe_result$Cx, sawe_input$Cz, sawe_result$Cz, sawe_input$sigma_mass),
sigma_kz = sigma_rog(sawe_result$mass, sawe_result$Izz, sawe_result$sigma_mass, sawe_result$sigma_Izz,
sawe_input$Cx, sawe_result$Cx, sawe_input$Cy, sawe_result$Cy, sawe_input$sigma_mass)
)
The independently-calculated radii of gyration uncertainties are
sigma_ROG
This r agree(isTRUE(all.equal(sigma_ROG, c(sigma_kx = sawe_result$sigma_kx, sigma_ky = sawe_result$sigma_ky, sigma_kz = sawe_result$sigma_kz)))) with the computed result.
result <- rollup_radii_of_gyration_unc(sawe_tree, add_radii_of_gyration(rollup_mass_props_and_unc(sawe_tree, sawe_table)))
parts <- result[1:2, ]
aggr <- result[3, ]
# does not agree
sigma_rog_sawe<- function(k, kt, v1, c1, v2, c2, sigma_m, m, mt, sigma_v1, sigma_v2, sigma_k) {
sqrt(sum(
((k^2 - kt^2 + (v1 - c1)^2 + (v2 - c2)^2) * sigma_m)^2 +
(2 * m * (v1 - c1) * sigma_v1)^2 +
(2 * m * (v2 - c2) * sigma_v2)^2 +
(2 * m * k * sigma_k)^2
)) / (2 * kt * mt)
}
sigma_ROG_sawe <- c(
sigma_kx = sigma_rog_sawe(parts$kx, aggr$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cy, parts$sigma_Cz, parts$sigma_kx),
sigma_ky = sigma_rog_sawe(parts$ky, aggr$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cz, parts$sigma_ky),
sigma_kz = sigma_rog_sawe(parts$kz, aggr$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cy, parts$sigma_kz)
)
# agrees
alt_moi_sawe <- function(m, mt, k, v1, c1, v2, c2) {
sum(m * (k^2 + v1^2 + v2^2)) - mt * (c1^2 + c2^2)
}
alt_MOI_sawe <- c(
Ixx = alt_moi_sawe(parts$mass, aggr$mass, parts$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz),
Iyy = alt_moi_sawe(parts$mass, aggr$mass, parts$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz),
Izz = alt_moi_sawe(parts$mass, aggr$mass, parts$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy)
)
# does not agree
alt_sigma_moi_sawe <- function(k, v1, c1, v2, c2, sigma_m, m, mt, sigma_v1, sigma_v2, sigma_k) {
sqrt(sum(
((k^2 + (v1 - c1)^2 + (v2 - c2)^2) * sigma_m)^2 +
(2 * m * (v1 - c1) * sigma_v1)^2 +
(2 * m * (v2 - c2) * sigma_v2)^2 +
(2 * m * k * sigma_k)^2
))
}
alt_sigma_MOI_sawe <- c(
sigma_Ixx = alt_sigma_moi_sawe(parts$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cy, parts$sigma_Cz, parts$sigma_kx),
sigma_Iyy = alt_sigma_moi_sawe(parts$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cz, parts$sigma_ky),
sigma_Izz = alt_sigma_moi_sawe(parts$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy,
parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cy, parts$sigma_kz)
)
Comparison With Published Results
The SAWE reference provides computed results for their example (excluding radii of gyration and their uncertainties). These results match those within a tolerance of 0.2%. The small differences are likely because their actual input values were not identical to the truncated values published in the article.
Performance Evaluation
mp_tree_depth <- max(dfs(mp_tree, 2, mode = "in", order=FALSE, dist=TRUE)$dist) + 1
nv <- length(igraph::V(mp_tree))
ne <- length(igraph::E(mp_tree))
nl <- length(which(igraph::degree(mp_tree, mode="in") == 0))
ni <- (nv - 1) * 20
no <- (nv - nl) * 20
mp_table and mp_tree are a synthesized data set representing a tree of depth r mp_tree_depth with r nv vertices and r ne edges. r nl vertices are leaves, the remaining r nv - nl are non-leaves. Rolling up mass properties and uncertainties for this data set combines r paste(ni) input values to produce r no output values. Mass properties alone halves those values.
Benchmarks were taken on a platform with these CPU characteristics:
Python Version: 3.12.7.final.0 (64 bit)
Cpuinfo Version: 9.0.0
Vendor ID Raw:
Hardware Raw:
Brand Raw: Apple M3
Hz Advertised Friendly:
Hz Actual Friendly: 2.4000 GHz
Hz Advertised:
Hz Actual: (2400000000, 0)
Arch: ARM_8
Bits: 64
Count: 8
Arch String Raw: arm64
L1 Data Cache Size:
L1 Instruction Cache Size:
L2 Cache Size:
L2 Cache Line Size:
L2 Cache Associativity:
L3 Cache Size:
Stepping:
Model:
Family: 6
Processor Type:
Flags: acpi, aes, apic, clfsh, cmov, cx16, cx8, de, ds, dscpl, dtse64, est, fpu, fxsr, htt, mca, mce, mmx, mon, msr, mtrr, pae, pat, pbe, pclmulqdq, pdcm, pge, pse, pse36, seglim64, sep, ss, sse, sse2, sse3, sse4.1, sse4.2, ssse3, tm, tm2, tpr, tsc, vme, vmx
Benchmark results for rollup of mass properties and uncertainties were taken with and without input validation:
benchmark('mp + unc no validation' = rollup_mass_props_and_unc_fast(mp_tree, mp_table),
'mp + unc validation' = rollup_mass_props_and_unc(mp_tree, mp_table),
'mp no validation' = rollup_mass_props_fast(mp_tree, mp_table),
'mp validation' = rollup_mass_props(mp_tree, mp_table),
replications = 1,
columns = c("test", "replications", "elapsed", "user.self", "sys.self")
)
Times reported are in seconds.
test replications elapsed user.self sys.self
4 mp validation 1 0.840 0.824 0.016
3 mp no validation 1 0.604 0.594 0.011
2 mp + unc validation 1 1.479 1.454 0.026
1 mp + unc no validation 1 0.993 0.974 0.019
References
Try the massProps package in your browser
Any scripts or data that you put into this service are public.
massProps documentation built on April 4, 2025, 1:45 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
\newcommand{\sumwi}{\sum_{i=1}^{n}w_i}
\newcommand{\sumwivi}[1]{\sum_{i=1}^{n}{w_i}{{#1}_i}}
\newcommand{\Ixx}{I_{XX}} \newcommand{\Iyy}{I_{YY}} \newcommand{\Izz}{I_{ZZ}} \newcommand{\Ixy}{I_{XY}} \newcommand{\Ixz}{I_{XZ}} \newcommand{\Iyz}{I_{YZ}}
\newcommand{\kx}{k_X} \newcommand{\ky}{k_Y} \newcommand{\kz}{k_Z}
\newcommand{\moia}[3]{ \sum_{i=1}^{n} \left[ {I_{#1}}_i + w_i \left( {#2}_i^2 + {#3}_i^2 \right) - w_i \left( \bar{#2}^2 + \bar{#3}^2 \right) \right] }
\newcommand{\moib}[3]{ \sum_{i=1}^{n} \left{ {I_{#1}}_i + w_i \left[ \left( {#2}_i - \bar{#2} \right)^2 + \left( {#3}_i - \bar{#3} \right)^2 \right] \right} }
\newcommand{\di}[1]{({#1}_i - \bar{#1})}
\newcommand{\poia}[3]{\sum_{i=1}^{n} \left[ {I_{#1}}i + w_i {{#2}_i}{{#3}_i} -w_i (\bar{#2}\bar{#3})\right]} \newcommand{\poib}[3]{\sum{i=1}^{n} \left[ {I_{#1}}i + w_i \di{#2}\di{#3}\right]} \newcommand{\sigmacm}[1]{ \sqrt{ \sum{i=1}^n \left{ (w_i {{\sigma_\bar{#1}}}i)^2 + [{\sigma_w}_i \di{#1}]^2 \right} } \bigg/ \sum{i=1}^{n}w_i \ } \newcommand{\sigmamoi}[3]{\sqrt{ \sum_{i=1}^n \left{ \sigma_{{I_{#1}}i}^2 + \left[ 2 w_i \di{#2} \sigma{{#2}i} \right]^2 + \left[ 2 w_i \di{#3} \sigma{{#3}i} \right]^2 + \left[ \left(\di{#2}^2 + \di{#3}^2 \right)\sigma{w_i}\right]^2 \right} } } \newcommand{\sigmapoi}[3]{\sqrt{ \sum_{i=1}^n \left{ \sigma_{{I_{#1}}i}^2 + \left[ \di{#2} w_i \sigma{{#3}i}\right]^2 + \left[ \di{#2}\di{#3}\sigma{w_i} \right]^2 + \left[ \di{#3} w_i \sigma_{{#2}_i} \right]^2 \right} } }
Overview
The massProps package extends rollupTree with functions to recursively calculate mass properties (and optionally, their uncertainties) for arbitrary decomposition trees. Formulas implemented are described in a technical paper published by the Society of Allied Weight Engineers [@zimmerman:05:sawe].
Synopsis
Data Structures
massProps operates on two fundamental data structures: a mass properties table and a tree. The mass properties table has an entry for every item in a tree structure of items; the edges of the tree convey the parent-child relations among items. The two data structures are linked by the id column of the data frame, which must be a character vector of unique item identifiers, and the vertex names of the tree. The sets of identifiers must be identical.
Mass Property Table
Required Columns for Mass Properties
The Mass Property Table must contain the following columns. Other columns may exist and will remain unmodified.
-
idunique identifier for each item (row) -
massmass of the item (numeric) -
Cx$x$-component of center of mass (numeric) -
Cy$y$-component of center of mass (numeric) -
Cx$z$-component of center of mass (numeric) -
Ixxmoment of inertia about the $x$ axis (numeric) -
Iyymoment of inertia about the $y$ axis (numeric) -
Izzmoment of inertia about the $z$ axis (numeric) -
Ixyproduct of inertia relative to the $x$ and $y$ axes (numeric) -
Ixzproduct of inertia relative to the $x$ and $z$ axes (numeric) -
Iyzproduct of inertia relative to the $y$ and $z$ axes (numeric) -
POIconveither '+' or '-', indicating the sign convention for products of inertia. In the negative convention, for example, $I_{XY} \equiv -\int{xy \rho \, dV}$. In the positive convention, $I_{XY} \equiv \int{xy \rho \, dV}$. -
Ipointlogical indicator that this item is considered a point mass. The same algebraic result can be achieved by setting all moments and products of inertia to zero, butrollup_mass_props()by default ensures that all leaf items in the tree have mass properties that correspond to physically-realizable objects. A zero inertia tensor will fail this check. Rather than relax the check (which is essential for trustworthy results), aTRUEvalue forIpointindicates that the inertia tensor should be excluded from computations.
Required Columns for Mass Properties Uncertainties
The following columns are required for uncertainty calculations:
-
sigma_massmass uncertainty (numeric) -
sigma_Cx$x$-component of center of mass uncertainty (numeric) -
sigma_Cy$y$-component of center of mass uncertainty (numeric) -
sigma_Cx$z$-component of center of mass uncertainty (numeric) -
sigma_Ixxmoment of inertia about the $x$ axis uncertainty (numeric) -
sigma_Iyymoment of inertia about the $y$ axis uncertainty (numeric) -
sigma_Izzmoment of inertia about the $z$ axis uncertainty (numeric) -
sigma_Ixyproduct of inertia relative to the $x$ and $y$ axes uncertainty (numeric) -
sigma_Ixzproduct of inertia relative to the $x$ and $z$ axes uncertainty (numeric) -
sigma_Iyzproduct of inertia relative to the $y$ and $z$ axes uncertainty (numeric)
It is the caller's responsibility to ensure that all values are expressed in appropriate and compatible units.
Tree
The tree is an igraph::graph with vertices named by identifiers in the mass properties table. It can be of arbitrary depth and shape as long as it satisfies certain well-formedness properties:
-
it is connected and acyclic (as an undirected graph), i.e., it is a tree
-
it is directed, with edge direction going from child to parent
-
it contains neither loops (self-edges) nor multiple edges
-
it contains a single root vertex (i.e., one whose out degree is zero)
Invocation
library(rollupTree) library(massProps) suppressPackageStartupMessages({library(igraph)})
Suppose we have the following mass properties table:
test_table
Suppose we also have this tree:
test_tree
plot(test_tree,layout=layout_as_tree(test_tree, 2, mode="in"), vertex.shape = 'none', edge.arrow.mode = 0)
Then we can compute mass properties for non-leaf elements by calling rollup_mass_props():
rollup_mass_props(test_tree, test_table)
Note that, although the table shows the parent of each element for clarity of exposition, the child-parent relations are coneveyed only by the tree passed as the first argument.
The input may also contain uncertainties data. This example is from the Society of Allied Weight Engineers:
na_mass_props_and_unc <- function(d, t, v) { xyz <- c("x", "y", "z") list( mass = NA, center_mass = c(x = NA, y = NA, z = NA), inertia = matrix(nrow = 3, ncol = 3, dimnames = list(xyz, xyz)), poi_conv = "+", point = FALSE, sigma_mass = NA, sigma_center_mass = c(x = NA, y = NA, z = NA), sigma_inertia = matrix(nrow = 3, ncol = 3, dimnames = list(xyz, xyz)) ) } na_mass_props_and_unc_update <- function(d, t, s) { update_mass_props_and_unc(d, t, s, override = na_mass_props_and_unc) } sawe_input <- rollup(sawe_tree, sawe_table, update = na_mass_props_and_unc_update, validate_ds = validate_mass_props_and_unc_table)
sawe_input
rollup_mass_props_and_unc(sawe_tree, sawe_input)
Objectives and Strategy
The objective of this package is to provide a trustworthy, well-documented, reference implementation for computation of mass properties (and their uncertainties) of aggregate objects from those of their parts. Aggregation can be recursive (e.g., indentured Bill of Materials), so it must accommodate trees of arbitrary depth and shape.
Strategies for achieving the objective include
-
basing the calculations on published industry references,
-
re-casting those lengthy reference equations into concise vector or matrix forms to reduce the error surface for source code and exploit the capabilities of
R, which treats vectors and matrices as first-class objects, -
delegating orchestration to the
rollupTreepackage, which, among other things, verifies that the input tree is well-formed and ensures proper ordering of computations, -
ensuring that all asserted leaf mass properties and uncertainties correspond to physically-realizable objects,
-
coding in pure functional style, (i.e., avoiding mutable variables, implying iteration with
Map()andReduce()), and -
covering the entire code base with unit tests.
The author has intentionally made no effort to micro-optimize for performance. In particular, the author is aware that representing the inertia and its uncertainty as 3 ⨉ 3 matrices is "inefficient" to the degree that it independently calculates values that are redundant by symmetry. "Inefficient", however, does not mean "slow". See [Performance Evaluation] below.
Theory
In this section, we state the reference equations [@zimmerman:05:sawe] and show, where applicable, how those equations can be rewritten in more concise form. The form of the equations actually implemented is displayed within a box, e.g. $\boxed{F = ma}$.
The reference uses the word weight and the symbol $w$ in equations. We interpret weight as mass. The reference refers to center of mass by its $x$, $y$, and $z$ components. Symbols for moments ($\Ixx$) and products ($\Ixy$) of inertia are conventional. Variables with $i$ subscripts designate properties of parts; those without designate properties of aggregates. The letter $\sigma$ denotes uncertainty. $\sigma_w$, for example, is the mass uncertainty.
Mass Properties
Mass
The mass equation is suitable as is.
$$ \boxed{ w = \sum_{i=1}^{n}w_i } $$
The corresponding R code is
amp$mass <- Reduce(`+`, Map(f = function(mp) mp$mass, mpl))
In this and the following code snippets, the variable mpl is a list of input mass property sets for parts, the variable mp is a formal parameter of an anonymous function applied to each member of mpl, and amp is the resulting aggregate mass property set. The line above is an R functional programming idiom for "set the mass value of the aggregate to the sum of the mass values of the parts".
Center of Mass
$$ \begin{align} \bar{x} & = \sum_{i=1}^{n}{w_i}{{x}i} \bigg/ \sum{i=1}^{n}w_i\ \bar{y} & = \sum_{i=1}^{n}{w_i}{{y}i} \bigg/ \sum{i=1}^{n}w_i\ \bar{z} & = \sum_{i=1}^{n}{w_i}{{z}i} \bigg/ \sum{i=1}^{n}w_i\ \end{align} $$
We can express center of mass as a 3-vector:
$$ \boxed{ \begin{align} \boldsymbol{c}_i & = (x_i \quad y_i \quad z_i)^T \ \boldsymbol{\bar{c}}& = (\bar{x} \quad \bar{y} \quad \bar{z})^T \end{align} } $$
Then
$$ \boxed{ \boldsymbol{\bar{c}} = \frac{1}{w} \sum_{i=1}^{n}{w_i}{{\boldsymbol{c}}_i} } $$
The corresponding R code is
amp$center_mass <- Reduce(`+`, Map(f = function(mp) mp$mass * mp$center_mass, mpl)) / amp$mass
Inertia Tensor
Moments of Inertia
$$ \begin{align} \Ixx & = \moia{XX}{y}{z} & = \moib{XX}{y}{z} \ \Iyy & = \moia{YY}{x}{z} & = \moib{YY}{x}{z} \ \Izz & = \moia{ZZ}{x}{y} & = \moib{ZZ}{x}{y} \ \end{align} $$
Products of Inertia
$$ \begin{align} \Ixy & = \poia{XY}{x}{y} & = \poib{XY}{x}{y} \ \Ixz & = \poia{XZ}{x}{z} & = \poib{XZ}{x}{z} \ \Iyz & = \poia{YZ}{y}{z} & = \poib{YZ}{y}{z} \ \end{align} $$
Matrix Formulation
Let $\boldsymbol{I}$ be the inertia tensor of the aggregate and $\boldsymbol{I}_i$ be that of part $i$. The equations for products of inertia above clearly follow the positive integral convention, so
$$ \boldsymbol{I} = \left[ \begin{array}{rrr} \Ixx & -\Ixy & -\Ixz \ -\Ixy & \Iyy & -\Iyz \ -\Ixz & -\Iyz & \Izz \ \end{array} \right] $$
and similarly for $\boldsymbol{I}_i$.
Noting the repeated appearance of terms of the form $\di{x}\di{y}$, we form the outer product
$$ \boxed{ \begin{align} \boldsymbol{d}_i & = (\di{x} \quad \di{y} \quad \di{z})^T \ \boldsymbol{Q}_i & = \boldsymbol{d}_i {\boldsymbol{d}_i}^T \end{align} } $$ Then
$$ \begin{align} \boldsymbol{Q}_i & = \begin{bmatrix} \di{x}^2 & \di{x}\di{y} & \di{x}\di{z} \ \di{y}\di{x} & \di{y}^2 & \di{y}\di{z} \ \di{z}\di{x} & \di{z}\di{y} & \di{z}^2 \ \end{bmatrix} \end{align} $$
Let $\boldsymbol{s}_i$ be the matrix of inertia tensor summands from the reference equations. That is,
$$ \boldsymbol{I} = \sum_{i=1}^{n} \boldsymbol{s}_i $$
where $$ \begin{align} \boldsymbol{s}_i & = \boldsymbol{I}_i \ & - w_i \begin{bmatrix} -\di{y}^2 - \di{z}^2 & \di{x}\di{y} & \di{x}\di{z} \ \di{x}\di{y} & -\di{x}^2 - \di{z}^2 & \di{y}\di{z} \ \di{x}\di{z} & \di{y}\di{z} & -\di{x}^2 - \di{y}^2 \ \end{bmatrix} \ & = \boldsymbol{I}_i \ & - w_i \begin{bmatrix} \di{x}^2 & \di{x}\di{y} & \di{x}\di{z} \ \di{x}\di{y} & \di{y}^2 & \di{y}\di{z} \ \di{x}\di{z} & \di{y}\di{z} & \di{z}^2 \ \end{bmatrix} \ & - w_i \begin{bmatrix} -\di{x}^2 - \di{y}^2 - \di{z}^2 & 0 & 0 \ 0 & -\di{x}^2 - \di{y}^2 - \di{z}^2 & 0 \ 0 & 0 &-\di{x}^2 - \di{y}^2 - \di{z}^2 \ \end{bmatrix} \ & = \boldsymbol{I}_i - w_i \left( \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i) \boldsymbol{1}_3 \right) \end{align} $$
where $\mathrm{tr}(\boldsymbol{Q}_i)$ is the trace of $\boldsymbol{Q}_i$, i.e., the sum of its diagonal elements, and $\boldsymbol{1}_3$ is the 3⨉3 identity matrix. Therefore
$$ \boxed{ \boldsymbol{I} = \sum_{i=1}^{n} \left( \boldsymbol{I}_i - w_i \boldsymbol{M}_i \right) } $$ where
$$ \boxed{ \boldsymbol{M}_i = \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i) \boldsymbol{1}_3 } $$
The corresponding R code is
amp$inertia <- Reduce(`+`, Map( f = function(mp) { d <- amp$center_mass - mp$center_mass Q <- outer(d, d) M <- Q - sum(diag(Q)) * diag(3) if (mp$point) -mp$mass * M else mp$inertia - mp$mass * M }, mpl ))
Mass Property Uncertainties
Mass Uncertainty
The mass uncertainty equation is suitable as is.
$$ \boxed{ \sigma_w = \sqrt{ \sum_{i=1}^n {{\sigma_w}_i}^2 } } $$
The corresponding R code is
amp$sigma_mass = sqrt(Reduce(`+`, Map(f = function(v) v$sigma_mass^2, mpl)))
Center of Mass Uncertainty
$$ \begin{align} \sigma_\bar{x} & = \sigmacm{x} \ \sigma_\bar{y} & = \sigmacm{y} \ \sigma_\bar{z} & = \sigmacm{z} \ \end{align} $$
As before, we create a 3-vector for center of mass uncertainties. Let
$$ \boxed{ \begin{align} \boldsymbol{\sigma_c} & = (\sigma_\bar{x} \quad \sigma_\bar{y} \quad \sigma_\bar{z})^T \ {\boldsymbol{\sigma_c}}i & = ({\sigma\bar{x}}i \quad {\sigma\bar{y}}i \quad {\sigma\bar{z}}_i)^T \end{align} } $$
If we construe (as R does) squaring and taking square roots of vectors element-wise, then
$$ \boxed{ \boldsymbol{\sigma_c} = \frac{1}{w} \sqrt{ \sum_{i=1}^n \left{ (w_i {{\boldsymbol{\sigma_c}}}_i)^2 + [{\sigma_w}_i ({\boldsymbol{c}}_i - \bar{\boldsymbol{c}})]^2 \right}} } $$
The corresponding R code is
amp$sigma_center_mass = sqrt(Reduce(`+`, Map( f = function(v) { (v$mass * v$sigma_center_mass)^2 + (v$sigma_mass * (v$center_mass - amp$center_mass))^2 }, mpl ))) / amp$mass
Inertia Tensor Uncertainty
Moments of Inertia Uncertainties
$$ \begin{align} \sigma_{\Ixx} & = \sigmamoi{XX}{y}{z} \ \sigma_{\Iyy} & = \sigmamoi{YY}{x}{z} \ \sigma_{\Izz} & = \sigmamoi{ZZ}{x}{y} \ \end{align} $$
Products of Inertia Uncertainties
$$ \begin{align} \sigma_{\Ixy} & = \sigmapoi{XY}{x}{y} \ \sigma_{\Ixz} & = \sigmapoi{XZ}{x}{z} \ \sigma_{\Iyz} & = \sigmapoi{YZ}{y}{z} \ \end{align} $$
Matrix Formulation
Let
$$ \boxed{ \begin{align} \boldsymbol{d}i & = (\di{x} \quad \di{y} \quad \di{z})^T \ {\boldsymbol{\sigma_c}}_i & = ({\sigma\bar{x}}i \quad {\sigma\bar{y}}i \quad {\sigma\bar{z}}_i)^T \ \boldsymbol{P}_i & = \boldsymbol{d}_i {\boldsymbol{\sigma_c}}_i^T \ \boldsymbol{Q}_i & = \boldsymbol{d}_i {\boldsymbol{d}_i}^T \end{align} } $$
Then
$$ \begin{align} \boldsymbol{P}i & = \begin{bmatrix} \di{x}\sigma{x_i} &\di{x}\sigma_{y_i} &\di{x}\sigma_{z_i} \ \di{y}\sigma_{x_i} & \di{y}\sigma_{y_i} & \di{y}\sigma_{z_i} \ \di{z}\sigma_{x_i} & \di{z}\sigma_{y_i} & \di{z}\sigma_{z_i} \ \end{bmatrix} \ \ \boldsymbol{Q}_i & = \begin{bmatrix} \di{x}^2 &\di{x}\di{y} &\di{x}\di{z} \ \di{y}\di{x} & \di{y}^2 & \di{y}\di{z} \ \di{z}\di{x} & \di{z}\di{y} & \di{z}^2 \ \end{bmatrix} \end{align} $$
Let $\boldsymbol{s}_i^2$ be the matrix of inertia tensor uncertainty summands in the standard formulas for a given subcomponent $i$ above. That is,
$$ \boldsymbol{\sigma_I}^2 = \sum_{i=1}^{n} \boldsymbol{s}_i^2 $$
Let ${p_X}i$, ${p_Y}_i$, and ${p_Z}_i$ be the respective diagonal elements of $P_i$. Let $\boldsymbol{1}_3$ be the 3 ⨉ 3 identity matrix. If we interpret squaring a matrix as the Hadamard (element-wise) product with itself, then $$ \begin{align} \boldsymbol{s}_i^2 & = {\boldsymbol{\sigma_I}}_i^2 \ & + \begin{bmatrix} 2 w_i \di{y} \sigma{y_i} & w_i\di{x} \sigma_{y_i} & w_i\di{x} \sigma_{z_i} \ w_i\di{x} \sigma_{y_i} & 2 w_i\di{x} \sigma_{x_i} & w_i \di{y} \sigma_{z_i} \ w_i\di{x} \sigma_{z_i} & w_i \di{y} \sigma_{z_i} & 2 w_i\di{x} \sigma_{x_i} \end{bmatrix}^2 \ & + \begin{bmatrix} 2 w_i \di{z} \sigma_{z_i} & w_i \di{y} \sigma_{x_i} & w_i \di{z} \sigma_{x_i} \ w_i \di{y} \sigma_{x_i} & 2 w_i \di{z} \sigma_{z_i} & w_i \di{z} \sigma_{y_i} \ w_i \di{z} \sigma_{x_i} & w_i \di{z} \sigma_{y_i} & 2 w_i \di{y} \sigma_{y_i} \end{bmatrix}^2 \ & + \begin{bmatrix} (\di{y}^2 + \di{z}^2)\sigma_{w_i} &\di{x}\di{y}\sigma_{w_i} &\di{x}\di{z}\sigma_{w_i} \ \di{y}\di{x}\sigma_{w_i} & (\di{x}^2 + \di{z}^2)\sigma_{w_i} & \di{y}\di{z}\sigma_{w_i} \ \di{z}\di{x}\sigma_{w_i} & \di{z}\di{y}\sigma_{w_i} & (\di{x}^2 + \di{y}^2)\sigma_{w_i} \ \end{bmatrix}^2 \ \ & = {\boldsymbol{\sigma_I}}i^2 \ & + w_i^2 \left( \boldsymbol{P}_i - \begin{bmatrix} \di{x}\sigma{x_i} - 2 \di{y}\sigma_{y_i} & 0 & 0 \ 0 & \di{y}\sigma_{y_i} - 2\di{x}\sigma_{x_i} & 0 \ 0 & 0 & \di{z}\sigma_{y_i} - 2\di{x}\sigma_{x_i} \ \end{bmatrix} \right) ^2 \ & + w_i^2 \left( \boldsymbol{P}i^T - \begin{bmatrix} \di{x}\sigma{x_i} - 2 \di{z}\sigma_{y_i} & 0 & 0 \ 0 & \di{y}\sigma_{y_i} - 2 \di{z}\sigma_{z_i} & 0 \ 0 & 0 & \di{z}\sigma_{y_i} - 2 \di{y}\sigma_{y_i} \ \end{bmatrix} \right) ^2 \ & + \sigma_{w_i}^2 \left( \boldsymbol{Q}i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3 \right)^2 \ \ & = {\boldsymbol{\sigma_I}}_i^2 \ & + w_i^2 \left( \boldsymbol{P}_i - \begin{bmatrix} {p_X}_i - 2 {p_Y}_i & 0 & 0 \ 0 & {p_Y}_i - 2 {p_X}_i & 0 \ 0 & 0 & {p_Z}_i - 2 {p_X}_i \ \end{bmatrix} \right) ^2 \ & + w_i^2 \left( \boldsymbol{P}_i^T - \begin{bmatrix} {p_X}_i - 2{p_Z}_i & 0 & 0 \ 0 & {p_Y}_i - 2 {p_Z}_i & 0 \ 0 & 0 & {p_Z}_i - 2 {p_Y}_i \ \end{bmatrix} \right) ^2 \ & + \sigma{w_i}^2 \left( \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3 \right)^2 \end{align} $$
Finally,
$$ \boxed{ \boldsymbol{\sigma_I} = \sqrt{ \sum_{i=1}^{n} {\left{ {\boldsymbol{\sigma_I}}_i^2 +{\boldsymbol{M}_4}_i \right}} } } $$
where
$$ \boxed{ \begin{align} {\boldsymbol{M}1}_i & = \boldsymbol{P}_i - \begin{bmatrix} {p_X}_i - 2{p_Y}_i & 0 & 0 \ 0 & {p_Y}_i - 2{p_X}_i & 0 \ 0 & 0 & 2{p_Z}_i - 2{p_X}_i \ \end{bmatrix} \ {\boldsymbol{M}_2}_i & = \boldsymbol{P}_i^T - \begin{bmatrix} {p_X}_i - 2{p_Z}_i & 0 & 0 \ 0 & {p_Y}_i - 2{p_Z}_i & 0 \ 0 & 0 & 2{p_Z}_i - 2{p_Y}_i \ \end{bmatrix} \ {\boldsymbol{M}_3}_i & = \boldsymbol{Q}_i - \mathrm{tr}(\boldsymbol{Q}_i)\boldsymbol{1}_3 \ {\boldsymbol{M}_4}_i & = w_i^2 \left( {{\boldsymbol{M}_1}_i}^2 + {{\boldsymbol{M}_2}_i}^2 \right) + \left( \sigma{w_i} {\boldsymbol{M}_3}_i \right)^2 \end{align} } $$
The corresponding R code is
amp$sigma_inertia = sqrt(Reduce(`+`, Map( f = function(v) { d <- v$center_mass - amp$center_mass P <- outer(d, v$sigma_center_mass) p <- diag(P) Q <- outer(d, d) M1 <- P - diag(p - 2 * p[c("y", "x", "x")]) M2 <- t(P) - diag(p - 2 * p[c("z", "z", "y")]) M3 <- Q - sum(diag(Q)) * diag(3) M4 <- v$mass^2 * (M1^2 + M2^2) + (v$sigma_mass * M3)^2 if (v$point) M4 else v$sigma_inertia^2 + M4 }, mpl )))
Radii of Gyration and Their Uncertainties
By definition:
$$ \begin{align} \kx & = \sqrt{\Ixx / w} \ \ky & = \sqrt{\Iyy / w} \ \kz & = \sqrt{\Izz / w} \ \end{align} $$ Let
$$ \boxed{ \begin{align} \boldsymbol{k} & = (\kx \quad \ky \quad \kz)^T \ \boldsymbol{I} & = (\Ixx \quad \Iyy \quad \Izz)^T \end{align} } $$ Then
$$ \boxed{ \boldsymbol{k} = \sqrt{ \boldsymbol{I} / w } } $$
The corresponding R code is
rg <- get_mass_props(d, i) rg$radii_gyration <- sqrt(diag(rg$inertia) / rg$mass)
The SAWE reference gives equations for uncertainties of radii of gyration in recursive form, but as these radii are simply functions of moments of inertia and mass, we should be able to express their uncertainties in terms of uncertainties of moments of inertia and mass by applying standard uncertainty propagation theory [@wikipedia:uncprop].
Let
$$ \boxed{ \begin{align} \sigma_\boldsymbol{k} & = (\sigma_{\kx} \quad \sigma_{\ky} \quad \sigma_{\kz})^T \ \sigma_\boldsymbol{I} & = (\sigma_{\Ixx} \quad \sigma_{\Iyy} \quad \sigma_{\Izz})^T \end{align} } $$ Then
$$ \begin{align} {\sigma_\boldsymbol{k}}^2 & \approx \left( \frac{\partial \boldsymbol{k}}{\partial \boldsymbol{I}} \right)^2 {\sigma_\boldsymbol{I}}^2 + \left( \frac{\partial \boldsymbol{k}}{\partial w} \right)^2 {\sigma_w}^2 + 2 \frac{\partial \boldsymbol{k}}{\partial \boldsymbol{I}} \frac{\partial \boldsymbol{k}}{\partial w} \sigma_{\boldsymbol{I}w} \ & \approx \left( \frac{1}{2 \sqrt{w \boldsymbol{I}}} \right)^2 {\sigma_\boldsymbol{I}}^2 + \left(\frac{- \sqrt{\boldsymbol{I}}}{2 w^{3/2}}\right)^2 {\sigma_w}^2 + 2 \left( \frac{1}{2 \sqrt{w \boldsymbol{I}}} \right) \left(\frac{- \sqrt{\boldsymbol{I}}}{2 w^{3/2}}\right) \sigma_{\boldsymbol{I}w} \ & \approx \frac{1}{4 w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{4 w^3} {\sigma_w}^2 - \frac{1}{2 w^2} \sigma_{\boldsymbol{I}w} \ \end{align} $$ where $\sigma_{\boldsymbol{I}w}$ is the covariance between $\boldsymbol{I}$ and $w$.
Let $W$ and $X$ be random variables such that
$$ \begin{align} W & = \sum_{i=1}^{n} (w_i + {\epsilon_w}i) \ X & = \sum{i=1}^{n} \left{ {\Ixx}i + {\epsilon\Ixx}i + (w_i + {\epsilon_w}_i) \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \ \end{align} $$ where $E[{\epsilon_w}_i] = 0$ and $E[{\epsilon\Ixx}i] = 0$ for all $i$, $E[{\epsilon_w}_i {\epsilon_w}_j] = 0$ and $E[{\epsilon\Ixx}i {\epsilon\Ixx}j] = 0$ for all $i \neq j$, and $E[{\epsilon_w}_i {\epsilon\Ixx}_j] = 0$ for all $i$ and $j$. It is clear from the linearity of $W$ and $X$ in $w_i$ and $\Ixx_i$ that
$$ E[W] = E \left[\sum_{i=1}^{n} (w_i + {\epsilon_w}i) \right] = w + \sum{i=1}^{n} {\epsilon_w}i = w $$ and $$ \begin{align} E[X] & = E \left[ \sum{i=1}^{n} \left{ {\Ixx}i + {\epsilon{\Ixx}}i + (w_i + {\epsilon_w}_i) \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right] \ & = \Ixx + E \left[ \sum{i=1}^{n} \left{ {\epsilon_{\Ixx}}_i + {\epsilon_w}_i \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right] \ & = \Ixx \end{align} $$ Therefore,
$$ \begin{align} \sigma_{\Ixx w} & = E \left[ \left( W - E[W] \right) \left( X - E[X] \right)\right] \ & = E \left[ \left( \sum_{i=1}^{n} {\epsilon_w}i \right) \left( \sum{i=1}^{n} \left{ {\epsilon_{\Ixx}}i + {\epsilon_w}_i \left[ \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right] \right} \right) \right] \ & = E \left[ \sum{i=1}^{n} {\epsilon_w}^2_i \left( \left( {y}i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \right] \ & = \sum{i=1}^{n} E \left[ {\epsilon_w}^2_i \right] \left( \left( {y}i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \ & = \sum{i=1}^{n} {\sigma_w}^2_i \left( \left( {y}_i - \bar{y} \right)^2 + \left( {z}_i - \bar{z} \right)^2 \right) \ \end{align} $$
Exploiting symmetry, let
$$ \boxed{ \begin{align} \boldsymbol{d}_i & = (({x}_i - \bar{x}) \quad ({y}_i - \bar{y}) \quad ({z}_i - \bar{z}))^T \ \boldsymbol{s}_i & = (\boldsymbol{d}_i^T\boldsymbol{d}_i \quad \boldsymbol{d}_i^T\boldsymbol{d}_i \quad \boldsymbol{d}_i^T\boldsymbol{d}_i)^T \ \end{align} } $$
$$ \begin{align} {\sigma_\boldsymbol{k}}^2 & \approx \frac{1}{4 w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{4 w^3} {\sigma_w}^2 - \frac{1}{2 w^2} \sum_{i=1}^{n} {\sigma_w}^2_i (\boldsymbol{s}_i - \boldsymbol{d}_i^2 ) \ \end{align} $$ Therefore, we define
$$ \boxed{ \begin{align} {\sigma_\boldsymbol{k}} & = \frac{1}{2} \sqrt{ \frac{1}{w \boldsymbol{I}} {\sigma_\boldsymbol{I}}^2 + \frac{\boldsymbol{I}}{w^3} {\sigma_w}^2 - \frac{2}{w^2} \sum_{i=1}^{n} {\sigma_w}^2_i (\boldsymbol{s}_i - \boldsymbol{d}_i^2 ) } \end{align} } $$
The corresponding R code is
amp <- get_mass_props_and_unc(ds, target) I <- diag(amp$inertia) sigma_I <- diag(amp$sigma_inertia) amp$sigma_radii_gyration <- sqrt( sigma_I^2 / (amp$mass * I) + (I * amp$sigma_mass^2) / amp$mass^3 - 2 / amp$mass^2 * Reduce( `+`, Map( f = function(s) { mp <- get_mass_props_and_unc(ds, s) d2 <- (mp$center_mass - amp$center_mass)^2 mp$sigma_mass^2 * (sum(d2) - d2) }, sources ), init <- c(0, 0, 0) ) ) / 2 ``` # Testing and Validation ## Comparison With Independently-Calculated Results In this section we will calculate the results for the SAWE example step by step and compare them with the package results. The inputs are: ```r sawe_input <- sawe_table[which(sawe_table$id != "Combined"), ] sawe_input
Our computed result is
t <- rollup_radii_of_gyration_unc(sawe_tree, add_radii_of_gyration( rollup_mass_props_and_unc(sawe_tree, sawe_table) ) ) sawe_result <- t[t$id == "Combined", ] sawe_result
Mass
agree <- function(l) if (l) "agrees" else stop("disagreement")
mass <- sum(sawe_input$mass)
The independently-calculated mass is
mass
This r agree(isTRUE(all.equal(mass, sawe_result$mass))) with the computed result.
Center of Mass
C <- apply(sawe_input$mass / mass * sawe_input[, c("Cx", "Cy", "Cz")], 2, sum)
The independently-calculated center of mass is
C
This r agree(isTRUE(all.equal(unname(C), c(sawe_result$Cx, sawe_result$Cy, sawe_result$Cz)))) with the computed result.
Moments of Inertia
moi <- function(I, v1, v2, m, c1, c2) { sum(I + m * ((v1^2 + v2^2) - (c1^2 + c2^2))) } MOI <- c( Ixx = moi(sawe_input$Ixx, sawe_input$Cy, sawe_input$Cz, sawe_input$mass, C["Cy"], C["Cz"]), Iyy = moi(sawe_input$Iyy, sawe_input$Cx, sawe_input$Cz, sawe_input$mass, C["Cx"], C["Cz"]), Izz = moi(sawe_input$Izz, sawe_input$Cx, sawe_input$Cy, sawe_input$mass, C["Cx"], C["Cy"]) )
The independently-calculated moments of inertia are
MOI
This r agree(isTRUE(all.equal(MOI, c(Ixx = sawe_result$Ixx, Iyy = sawe_result$Iyy, Izz = sawe_result$Izz)))) with the computed result.
Products of Inertia
poi <- function(I, v1, v2, m, c1, c2) { sum(I + m * (v1 * v2 - c1 * c2)) } POI <- c( Ixy = poi(sawe_input$Ixy, sawe_input$Cx, sawe_input$Cy, sawe_input$mass, C["Cx"], C["Cy"]), Ixz = poi(sawe_input$Ixz, sawe_input$Cx, sawe_input$Cz, sawe_input$mass, C["Cx"], C["Cz"]), Iyz = poi(sawe_input$Iyz, sawe_input$Cy, sawe_input$Cz, sawe_input$mass, C["Cy"], C["Cz"]) )
The independently-calculated products of inertia are
POI
This r agree(isTRUE(all.equal(POI, c(Ixy = sawe_result$Ixy, Ixz = sawe_result$Ixz, Iyz = sawe_result$Iyz)))) with the computed result.
Radii of Gyration
rog <- function(I, m) sqrt(I / m) ROG <- c( kx = rog(sawe_result$Ixx, sawe_result$mass), ky = rog(sawe_result$Iyy, sawe_result$mass), kz = rog(sawe_result$Izz, sawe_result$mass) )
The independently-calculated radii of gyration are
ROG
This r agree(isTRUE(all.equal(ROG, c(kx = sawe_result$kx, ky = sawe_result$ky, kz = sawe_result$kz)))) with the computed result.
Mass Uncertainty
sigma_mass <- sqrt(sum(sawe_input$sigma_mass^2))
The independently-calculated mass uncertainty is
sigma_mass
This r agree(isTRUE(all.equal(sigma_mass, sawe_result$sigma_mass))) with the computed result.
Center of Mass Uncertainty
sigma_cm <- function(m, sigma_v, sigma_m, v, c, mass) { sqrt(sum((m * sigma_v)^2 + (sigma_m * (v - c))^2)) / mass } sigma_C <- c( sigma_Cx = sigma_cm(sawe_input$mass, sawe_input$sigma_Cx, sawe_input$sigma_mass, sawe_input$Cx, C["Cx"], mass), sigma_Cy = sigma_cm(sawe_input$mass, sawe_input$sigma_Cy, sawe_input$sigma_mass, sawe_input$Cy, C["Cy"], mass), sigma_Cz = sigma_cm(sawe_input$mass, sawe_input$sigma_Cz, sawe_input$sigma_mass, sawe_input$Cz, C["Cz"], mass) )
The independently-calculated center of mass uncertainties are
sigma_C
This r agree(isTRUE(all.equal(sigma_C, c(sigma_Cx = sawe_result$sigma_Cx, sigma_Cy = sawe_result$sigma_Cy, sigma_Cz = sawe_result$sigma_Cz)))) with the computed result.
Moments of Inertia Uncertainties
sigma_moi <- function(sigma_I, mass, sigma_mass, v1, v2, c1, c2, sigma_v1, sigma_v2) { sqrt(sum( sigma_I^2 + (2 * mass * (v1 - c1) * sigma_v1)^2 + (2 * mass * (v2 - c2) * sigma_v2)^2 + (((v1 - c1)^2 + (v2 - c2)^2) * sigma_mass)^2 )) } sigma_MOI <- c( sigma_Ixx = sigma_moi(sawe_input$sigma_Ixx, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cy, sawe_input$Cz, C["Cy"], C["Cz"], sawe_input$sigma_Cy, sawe_input$sigma_Cz), sigma_Iyy = sigma_moi(sawe_input$sigma_Iyy, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx, sawe_input$Cz, C["Cx"], C["Cz"], sawe_input$sigma_Cx, sawe_input$sigma_Cz), sigma_Izz = sigma_moi(sawe_input$sigma_Izz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx, sawe_input$Cy, C["Cx"], C["Cy"], sawe_input$sigma_Cx, sawe_input$sigma_Cy) )
The independently-calculated moments of inertia uncertainties are
sigma_MOI
This r agree(isTRUE(all.equal(sigma_MOI, c(sigma_Ixx = sawe_result$sigma_Ixx, sigma_Iyy = sawe_result$sigma_Iyy, sigma_Izz = sawe_result$sigma_Izz)))) with the computed result.
Products of Inertia Uncertainties
sigma_poi <- function(sigma_I, mass, sigma_mass, v1, v2, c1, c2, sigma_v1, sigma_v2) { sqrt(sum( sigma_I^2 + ((v1 - c1) * mass * sigma_v2)^2 + ((v1 - c1) * (v2 - c2) * sigma_mass)^2 + ((v2 - c2) * mass * sigma_v1)^2 )) } sigma_POI <- c( sigma_Ixy = sigma_poi(sawe_input$sigma_Ixy, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx, sawe_input$Cy, C["Cx"], C["Cy"], sawe_input$sigma_Cx, sawe_input$sigma_Cy), sigma_Ixz = sigma_poi(sawe_input$sigma_Ixz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cx, sawe_input$Cz, C["Cx"], C["Cz"], sawe_input$sigma_Cx, sawe_input$sigma_Cz), sigma_Iyz = sigma_poi(sawe_input$sigma_Iyz, sawe_input$mass, sawe_input$sigma_mass, sawe_input$Cy, sawe_input$Cz, C["Cy"], C["Cz"], sawe_input$sigma_Cy, sawe_input$sigma_Cz) )
The independently-calculated products of inertia uncertainties are
sigma_POI
This r agree(isTRUE(all.equal(sigma_POI, c(sigma_Ixy = sawe_result$sigma_Ixy, sigma_Ixz = sawe_result$sigma_Ixz, sigma_Iyz = sawe_result$sigma_Iyz)))) with the computed result.
Radii of Gyration Uncertainties
sigma_rog <- function(mt, It, sigma_mt, sigma_It, v1, c1, v2, c2, sigma_m) { sqrt(sigma_It^2 / (4 * mt * It) + (It * sigma_mt^2) / (4 * mt^3) - sum(sigma_m^2 * ((v1 - c1)^2 + (v2 - c2)^2)) / (2 * mt^2)) } sigma_ROG <- c( sigma_kx = sigma_rog(sawe_result$mass, sawe_result$Ixx, sawe_result$sigma_mass, sawe_result$sigma_Ixx, sawe_input$Cy, sawe_result$Cy, sawe_input$Cz, sawe_result$Cz, sawe_input$sigma_mass), sigma_ky = sigma_rog(sawe_result$mass, sawe_result$Iyy, sawe_result$sigma_mass, sawe_result$sigma_Iyy, sawe_input$Cx, sawe_result$Cx, sawe_input$Cz, sawe_result$Cz, sawe_input$sigma_mass), sigma_kz = sigma_rog(sawe_result$mass, sawe_result$Izz, sawe_result$sigma_mass, sawe_result$sigma_Izz, sawe_input$Cx, sawe_result$Cx, sawe_input$Cy, sawe_result$Cy, sawe_input$sigma_mass) )
The independently-calculated radii of gyration uncertainties are
sigma_ROG
This r agree(isTRUE(all.equal(sigma_ROG, c(sigma_kx = sawe_result$sigma_kx, sigma_ky = sawe_result$sigma_ky, sigma_kz = sawe_result$sigma_kz)))) with the computed result.
result <- rollup_radii_of_gyration_unc(sawe_tree, add_radii_of_gyration(rollup_mass_props_and_unc(sawe_tree, sawe_table))) parts <- result[1:2, ] aggr <- result[3, ] # does not agree sigma_rog_sawe<- function(k, kt, v1, c1, v2, c2, sigma_m, m, mt, sigma_v1, sigma_v2, sigma_k) { sqrt(sum( ((k^2 - kt^2 + (v1 - c1)^2 + (v2 - c2)^2) * sigma_m)^2 + (2 * m * (v1 - c1) * sigma_v1)^2 + (2 * m * (v2 - c2) * sigma_v2)^2 + (2 * m * k * sigma_k)^2 )) / (2 * kt * mt) } sigma_ROG_sawe <- c( sigma_kx = sigma_rog_sawe(parts$kx, aggr$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cy, parts$sigma_Cz, parts$sigma_kx), sigma_ky = sigma_rog_sawe(parts$ky, aggr$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cz, parts$sigma_ky), sigma_kz = sigma_rog_sawe(parts$kz, aggr$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cy, parts$sigma_kz) ) # agrees alt_moi_sawe <- function(m, mt, k, v1, c1, v2, c2) { sum(m * (k^2 + v1^2 + v2^2)) - mt * (c1^2 + c2^2) } alt_MOI_sawe <- c( Ixx = alt_moi_sawe(parts$mass, aggr$mass, parts$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz), Iyy = alt_moi_sawe(parts$mass, aggr$mass, parts$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz), Izz = alt_moi_sawe(parts$mass, aggr$mass, parts$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy) ) # does not agree alt_sigma_moi_sawe <- function(k, v1, c1, v2, c2, sigma_m, m, mt, sigma_v1, sigma_v2, sigma_k) { sqrt(sum( ((k^2 + (v1 - c1)^2 + (v2 - c2)^2) * sigma_m)^2 + (2 * m * (v1 - c1) * sigma_v1)^2 + (2 * m * (v2 - c2) * sigma_v2)^2 + (2 * m * k * sigma_k)^2 )) } alt_sigma_MOI_sawe <- c( sigma_Ixx = alt_sigma_moi_sawe(parts$kx, parts$Cy, aggr$Cy, parts$Cz, aggr$Cz, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cy, parts$sigma_Cz, parts$sigma_kx), sigma_Iyy = alt_sigma_moi_sawe(parts$ky, parts$Cx, aggr$Cx, parts$Cz, aggr$Cz, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cz, parts$sigma_ky), sigma_Izz = alt_sigma_moi_sawe(parts$kz, parts$Cx, aggr$Cx, parts$Cy, aggr$Cy, parts$sigma_mass, parts$mass, aggr$mass, parts$sigma_Cx, parts$sigma_Cy, parts$sigma_kz) )
Comparison With Published Results
The SAWE reference provides computed results for their example (excluding radii of gyration and their uncertainties). These results match those within a tolerance of 0.2%. The small differences are likely because their actual input values were not identical to the truncated values published in the article.
Performance Evaluation
mp_tree_depth <- max(dfs(mp_tree, 2, mode = "in", order=FALSE, dist=TRUE)$dist) + 1 nv <- length(igraph::V(mp_tree)) ne <- length(igraph::E(mp_tree)) nl <- length(which(igraph::degree(mp_tree, mode="in") == 0)) ni <- (nv - 1) * 20 no <- (nv - nl) * 20
mp_table and mp_tree are a synthesized data set representing a tree of depth r mp_tree_depth with r nv vertices and r ne edges. r nl vertices are leaves, the remaining r nv - nl are non-leaves. Rolling up mass properties and uncertainties for this data set combines r paste(ni) input values to produce r no output values. Mass properties alone halves those values.
Benchmarks were taken on a platform with these CPU characteristics:
Python Version: 3.12.7.final.0 (64 bit) Cpuinfo Version: 9.0.0 Vendor ID Raw: Hardware Raw: Brand Raw: Apple M3 Hz Advertised Friendly: Hz Actual Friendly: 2.4000 GHz Hz Advertised: Hz Actual: (2400000000, 0) Arch: ARM_8 Bits: 64 Count: 8 Arch String Raw: arm64 L1 Data Cache Size: L1 Instruction Cache Size: L2 Cache Size: L2 Cache Line Size: L2 Cache Associativity: L3 Cache Size: Stepping: Model: Family: 6 Processor Type: Flags: acpi, aes, apic, clfsh, cmov, cx16, cx8, de, ds, dscpl, dtse64, est, fpu, fxsr, htt, mca, mce, mmx, mon, msr, mtrr, pae, pat, pbe, pclmulqdq, pdcm, pge, pse, pse36, seglim64, sep, ss, sse, sse2, sse3, sse4.1, sse4.2, ssse3, tm, tm2, tpr, tsc, vme, vmx
Benchmark results for rollup of mass properties and uncertainties were taken with and without input validation:
benchmark('mp + unc no validation' = rollup_mass_props_and_unc_fast(mp_tree, mp_table), 'mp + unc validation' = rollup_mass_props_and_unc(mp_tree, mp_table), 'mp no validation' = rollup_mass_props_fast(mp_tree, mp_table), 'mp validation' = rollup_mass_props(mp_tree, mp_table), replications = 1, columns = c("test", "replications", "elapsed", "user.self", "sys.self") )
Times reported are in seconds.
test replications elapsed user.self sys.self 4 mp validation 1 0.840 0.824 0.016 3 mp no validation 1 0.604 0.594 0.011 2 mp + unc validation 1 1.479 1.454 0.026 1 mp + unc no validation 1 0.993 0.974 0.019
References
Try the massProps package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.