text:
\section{Introduction}

Machine learning (ML) has seen a surge in popularity in almost all fields of astronomy that involve massive amounts of complex data \citep{ball2010data, way2012advances, ivezic2014statistics, ntampaka2019role}.
Most ML methods target regression and classification, whose primary goal is to return a point estimate of an unknown response variable $\mathbf{y}$ given observed features $\mathbf{x}$, often falling short in quantifying nontrivial uncertainty in $\mathbf{y}$.
For instance, returning a point estimate for a supernova's type $\mathbf{y}$ given a supernova's light curve $\mathbf{x}$, or for a galaxy mass $\mathbf{y}$ given its light spectrum $\mathbf{x}$, fails to capture degeneracies in the mapping from $\mathbf{x}$ to $\mathbf{y}$.
Neglecting to propagate these uncertainties through down-stream  analyses may lead to imprecise or even inaccurate inferences of physical parameters.

Consider the following two examples of problems where uncertainty quantification can be impactful:
\begin{itemize}
    \item \textbf{Photometric redshift estimation.}
    In photometric redshift (photo-$z$) estimation, one attempts to constrain
    the cosmological {\em redshift} ($z$) of a galaxy  after observing the shifted spectrum using a handful of broadband filters, and sometimes additional variables such as morphology and environment.
    In a prediction setting, the response $\mathbf{y}$ could be the galaxy's true (i.e. spectroscopically observed) redshift but could also include other galaxy properties ($\alpha$) such as the galaxy's mass or age;
    the features $\mathbf{x}$ would represent the collection of directly observable inputs such as photometric magnitudes and colors used to predict $\mathbf{y} = z$ or more generally
    a multivariate response $\mathbf{y} = (z, \alpha)$.
    The use of photo-$z$\ posterior estimates --- that is, estimates of the probability density functions (PDFs) for individual galaxies --- is crucial for cosmological inference from photometric galaxy surveys, as $\gtrsim 99\%$ of currently cataloged galaxies are observed solely via photometry, a percentage that will only grow in the coming decade as the Large Synoptic Survey Telescope (\data{LSST}) begins gathering data \citep{Ivezic19}.
    In addition, widely different redshifts can be consistent with the observed colors of a galaxy;
    photo-$z$\ posterior estimates can capture such degeneracy or multimodality in the distribution whereas point estimates cannot.
    Though the benefits of using posteriors over $\alpha$ have yet to be fully exploited \citep{2018A&A...614A.129V}, it is thoroughly established that one can improve down-stream cosmological analysis by properly propagating photo-$z$\ estimate uncertainties via probability density functions (PDFs) rather than just using simple point predictions of $\mathbf{y}$ \citep{mandelbaum08precisionphotoz, wittman2009lies, sheldon2012photometric, carrasco2013tpz, graff2014skynet, desc_photoz}.
    \item \textbf{Forward-modeled observables in cosmology.}
    Some cosmological probes, such as the type Ia supernova (SN Ia) distance-redshift relationship, have observables that are straightforward to simulate in spite of an intractable likelihood.
    Likelihood-Free Inference (LFI) methods allow for parameter inference in settings where the likelihood function, which relates the observable  data $\mathbf{x}_{\rm obs}$ to the parameters of interest $\boldsymbol{\theta}$, is too complex to work with directly, but one is able to {\em simulate} $\mathbf{x}$ from a stochastic forward model at fixed parameter settings $\boldsymbol{\theta}$.
    The most common approach to LFI or simulation-based parameter inference is Approximate Bayesian Computation (ABC), whose many variations (see \citealt{beaumont2002approximate} and \citealt{sisson2018handbook} for a review) repeatedly draw from the simulator and compare the output with observed data $\ensuremath{\mathbf{x}_{\text{obs}}}$ in real time to arrive at a set of plausible parameter values consistent with $\ensuremath{\mathbf{x}_{\text{obs}}}$.
    With computationally intensive simulations, however, a classical ABC approach may not be practical.
    An alternative approach to ABC rejection sampling is to use faster training-based methods to assess the uncertainty about $\boldsymbol{\theta}$ for any $\mathbf{x}$ first, and then consider the specific case $\mathbf{x} = \ensuremath{\mathbf{x}_{\text{obs}}}$.
\end{itemize}

From a statistical perspective, the right tool for quantifying the uncertainty about $\mathbf{y}$ once $\mathbf{x}$ is observed is the {\em conditional density} $\pr{\mathbf{y} \gvn \mathbf{x}}$.
In a prediction context such as for photo-$z$\ problems, where heteroskedastic errors or multimodal responses may occur, conditional density estimation (CDE) of the density $\pr{z \gvn \mathbf{x}}$ for the redshift $z$ of individual galaxies given photometric data $\mathbf{x}$ provides a more nuanced accounting of uncertainty than a point estimate or prediction interval alone.
CDE can also be used in LFI where, in our notation, the parameters of interest $\boldsymbol{\theta}$ take the role of the ``response'' $y$.
In such settings, one can apply training-based approaches to forward-simulated data to estimate the posterior probability distribution $\pr{\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}}}$ of cosmological parameters $\boldsymbol{\theta}$ given observed data $\ensuremath{\mathbf{x}_{\text{obs}}}$, and from these posteriors then derive, e.g. posterior credible intervals of $\theta$.
Section \ref{sec:mult_response_mult_data} shows an example of cosmological parameter inference using mock weak lensing data.
Here we follow \citet{izbicki2014high, izbicki2019abc} to combine ABC and CDE by directly applying CDE techniques (Section \ref{sec:cde_methods}) and loss functions (Section \ref{sec:cde_loss}) to simulated data $\{(\boldsymbol{\theta}_i, \mathbf{x}_i)\}_{i=1}^{n}$.
Similar works include performing LFI via CDE using Gaussian copulas \citep{li2017gaussiancopula, chen2019gaussiancopula} and random forests \citep{marin2016abcrf}.
Other examples include neural density estimation in LFI via mixture density networks and masked autoregressive flows \citep{papamakarios2016fast, Lueckmann2017posterior, lueckmann2019likelihood, papamakarios2019sequential, greenberg2019automatic, alsing2018massive, alsing2019fast}.

Data in astronomy present a challenge to estimating conditional densities, due to both the complexity of the underlying physical processes and the complicated observational noise.
Precision cosmology, for example, requires combining data from different scientific probes, {each affected by unique sources of systematic uncertainty, to produce samples from} complicated joint likelihood functions {with nontrivial covariances} in a high-dimensional parameter space \citep{krause2017dark, Joudaki17, aghanim2018planck, vanUitert18, Abbott18}.
In such situations, CDE methods that target a variety of settings and non-standard data (images, correlation functions, mixed data types) become especially valuable.
However, for any given data type, there is no one-size-fits-all CDE method.
For example, deep neural networks often perform well in settings with large amounts of representative training data but in applications with smaller training samples one may need a different tool.
There is also additional value in models that are interpretable and easy to fit to the data at hand.

{\bf The goal of this paper is to provide statistical tools and open-source software for nonparametric CDE and method assessment appropriate for challenging data in a variety of inference scenarios.}

To our knowledge, existing CDE software either targets discrete $\mathbf{y}$ (e.g., probabilistic classifiers or ordinal classification \citep{eibe2001OrdinalClass}) or uses kernel density estimation (KDE)
across all data points to provide an estimate for a continuous $\mathbf{y}$ \citep{hdrcde}.
What distinguishes our methodology work from others is that we are able to adapt virtually any training-based prediction method to the problem of estimating full probability distributions.
By leveraging existing ML methods, we are hence able to provide uncertainty prediction methods and software for more general and complex data settings than `the computational tools currently available
in the literature.
In this paper, we showcase and provide \code{Python}\ and \code{R}\ code for four flexible CDE methods:
\code{NNKCDE}, \code{RFCDE}/\code{fRFCDE}, \code{FlexCode}\ and \code{DeepCDE}.\footnote{All of these methods, except for \code{DeepCDE}, have occurred in previously published or archived papers.
Hence, we only briefly review the methods in this paper and instead focus on software usage, algorithmic aspects, and the settings under which each method can be applied.}
Each CDE approach has particular usage for different settings of response dimensionality, feature types, and computational requirements.
Table \ref{tab:method-strenghts} provides a high-level overview in the top panel, and lists some properties of each method in the bottom panel.

A highlight of our software is that it makes uncertainty quantification straightforward for users of standard open-source machine learning \code{Python}\ packages.
As \code{NNKCDE}, \code{RFCDE}/\code{fRFCDE}, \code{FlexCode}\ share the \code{sklearn} \citep{Pedregosa2011sklearn} API (with \texttt{fit} and \texttt{predict} methods), our methods are usable within the \code{sklearn} ecosystem for, e.g., cross-validation and model selection.
In addition, \code{FlexCode}\ is essentially a plug-in method where the user can utilize
any \code{sklearn}-compatible regression function.
\code{DeepCDE}\ has implementations for both \code{Tensorflow} \citep{tensorflow2015-whitepaper} and \code{Pytorch} \citep{paszke2019pytorch}, two of the most widely used deep learning frameworks.

\begin{table}
\caption{\textit{Top:} Naming convention, high-level summary and hyper-parameters of CDE methods, along with references for further details and code examples.
\textit{Bottom:} Comparison of CDE methods in terms of training capacity and compatibility with multivariate response and different types of features, with capacities estimated based on input with around 100 features and a standard i5/i7/quad-core processor with $16$GB of RAM.
Note that less complex methods (such as \code{NNKCDE}) tend to be easier to use, easier to interpret, and often perform better in settings with smaller training sets, whereas more complex methods (such as \code{DeepCDE}) perform better in settings with larger (representative) training sets.}

\begin{tabular}{|c|c|l|l|c|}
\multicolumn{1}{c}{\textbf{Method}} & \multicolumn{1}{c}{\textbf{Name}}  & \multicolumn{1}{c}{\textbf{Summary}} & \multicolumn{1}{c}{\textbf{Hyper-parameters}} & \multicolumn{1}{c}{\textbf{Details}} \\ \hline
\multirow{4}{*}{\code{NNKCDE}} & & Computes a KDE estimate of & & \\
& \code{N}earest \code{N}eighbor & $\,$ multivariate $\mathbf{y}$ using the
nearest &  $\bullet$ Number of neighbors $k$ & Section \ref{sec:NNKCDE} \\
& $\,$ \code{K}ernel \code{CDE} & $\,$ neighbors of the evaluation point  $\mathbf{x}$ & $\bullet$ Kernel bandwidth $h$ & (Code: \ref{app:nnkcde})  \\
& & $\,$ in feature space. &  & \\ \hline
\multirow{5}{*}{\code{RFCDE}} &  & Random forest  that partitions & &  \\
&  &  $\,$ the feature space using a CDE loss.
 & $\bullet$ Random forest hyperparams. &  Section \ref{sec:RFCDE}  \\
& \code{R}andom \code{F}orest \code{CDE} & $\,$
Constructs a weighted KDE estimate & $\bullet$ Kernel bandwidth $h$ & (Code: \ref{app:rfcde}) \\
&  & $\,$ of multivariate $\mathbf{y}$ with weights & & \\
&  &  $\,$ defined by leaves in the forest. & & \\
\hline
\multirow{4}{*}{\code{fRFCDE}} &  & \code{RFCDE}\ version suitable for functional
 & $\bullet$ Random forest hyperparams. &  \\
& \code{f}unctional \code{R}andom &  $\,$ features $\mathbf{x}$. Partitions the feature & $\bullet$ Kernel bandwidth $h$ & Section \ref{sec:frfcde}\\
& $\,$ \code{F}orest \code{CDE} & $\,$ space directly rather than & $\bullet$ Partition parameter $\lambda$ & (Code: \ref{app:frfcde}) \\
& & $\,$ representing $\mathbf{x}$ as a vector. & & \\ \hline
\multirow{3}{*}{\code{FlexCode}} & \code{Flex}ible \code{Co}nditional & Uses basis expansion of univariate $y$ & $\bullet$ Number of expansion coeffs. & Section \ref{sec:FlexCode}  \\
& \code{D}ensity \code{E}stimation &  $\,$ to turn CDE into a series of &  $\bullet$ Selected regression method & (Code: \ref{app:flexcode})\\
&  & $\,$ univariate regression problems. & $\quad$ hyperparams. & \\  \hline
\multirow{5}{*}{\code{DeepCDE}} & & Uses basis expansion of univariate
 & & \\
 & \code{Deep} Neural & $\,$ $y$ similar to \code{FlexCode}, but learns
 & $\bullet$ Number of expansion coeffs. & Section \ref{sec:DeepCDE} \\
& Networks \code{CDE} &  $\,$ the expansion coefficients & $\bullet$ Selected deep neural network & (Code: \ref{sec:code_examples}) \\
& & $\,$ simultaneously using a deep &  $\quad$ architecture hyperparams. & \\
& & $\,$ neural network. & & \\ \hline
\end{tabular}

\vspace{0.25 cm}

\begin{minipage}[b]{0.12\textwidth}
\centering
    \begin{tikzpicture}
\coordinate [label=., align=right]  (A) at (0,0)  ;
\coordinate [label=., align=right]  (B) at (0,-2)  ;
\draw [->,very thick,
]($ (A)!1.5mm!(0,0) $) -- node [left, pos=0.5, align=center] {\textit{Method} \\ \textit{Complexity}}(B) ;
\end{tikzpicture
\end{minipage}
\begin{minipage}[b]{0.6\textwidth}
    \begin{tabular}{ccccc}
Method   &  Capacity (\# Training Pts) & Multivariate Response & Functional Features & Image Features \\ \hline
\code{NNKCDE}   & Up to $\sim 10^5$ & \checkmark             &                       &                   \\
(f)\code{RFCDE} & Up to $\sim 10^6$& \checkmark      & \checkmark            &                  \\
\code{FlexCode} & Up to $\sim 10^6$&    & \checkmark            &                  \\
\code{DeepCDE}  & Up to $\sim 10^8$ &    & \checkmark                      & \checkmark
\end{tabular}
\end{minipage}
\label{tab:method-strenghts}
\end{table}

In addition to the CDE methods above, we provide the package \code{cdetools}, which can be used for tuning and assessing the performance of CDE models on held-out validation data.
CDE method assessment is challenging per se because we never observe the true conditional probability density, merely samples (observations) from  it.
Furthermore, whereas loss functions such as the root-mean-squared error (RMSE) are typically used in regression problems, they are not appropriate for the task of uncertainty quantification of estimated probability densities.
The \code{cdetools}\ package provides two types of functions for method assessment.
First, it provides functions for computing a so-called CDE loss function (defined by Equation~\ref{eq:cde-loss} in Section \ref{sec:cde_loss}) for tuning and assessing the quality of individual PDFs.
Second, it provides diagnostic functions that probe the population-level performance of the PDFs.
More specifically, we have included functions for computing the Probability Integral Transform (PIT) and Highest Posterior Density (HPD);
these metrics check how well the final density estimates on average fit the data in  the tail and highest-density regions, respectively (see Section \ref{sec:hpd}, and Figure 2 for a visual sketch).

\textbf{Organization of the paper.}
In Section \ref{sec:cde_methods}, we introduce tools for conditional density estimation (\code{NNKCDE}, \code{RFCDE}/\code{fRFCDE}, \code{FlexCode}, \code{DeepCDE}).
In Section \ref{sec:assessment}, we describe tools for model selection and diagnostics.
Then, in  Section \ref{sec:examples}, we illustrate our CDE and method assessment tools for three different applications:
photo-$z$ estimation, likelihood-free cosmological inference and spec-$z$ estimation.
\code{Python}\ and \code{R}\ code usage examples can be found in Appendix \ref{sec:code_examples}.

\textbf{Notation.}
We denote the true (unknown) conditional density by $p(\mathbf{y}|\mathbf{x})$, whereas an estimate of the density is
denoted by $\hat p(\mathbf{y}|\mathbf{x})$.
We represent the CDE loss function that measures the discrepancy between the conditional density $p$ and its estimate $\hat p$ by $L(p,\hat p)$.
Typically this loss cannot be computed directly because it depends on unknown quantities;
an estimate of the  loss is denoted by $\hat L(p,\hat p)$.
As before, we continue to use bold-faced letters to denote vectors.

\section{Overview of Conditional Density Estimation Tools}
\label{sec:cde_methods}

We start by briefly describing the conditional density estimators in Table \ref{tab:method-strenghts}.
Unless otherwise stated, we choose the tuning or hyper-parameters by minimizing the CDE empirical loss in Equation~\ref{eq:estimated-cde-loss} using cross-validation.

\subsection{\code{NNKCDE}}
\label{sec:NNKCDE}

Nearest-Neighbors Kernel CDE (NNKCDE; \citealt{izbicki2017photo}, \citealt{freeman2017unified}) is a simple  and easily interpretable CDE method.
It computes a kernel density estimate of $\mathbf{y}$ using the $k$ nearest neighbors of the evaluation point $\mathbf{x}$.
The model has only two tuning parameters:
the number of nearest neighbors $k$ and the  bandwidth $h$ of the smoothing kernel in $\mathbf{y}$-space.
Both tuning parameters are chosen in a principled way by minimizing the CDE loss on validation data.

More specifically, the kernel density estimate of $\mathbf{y}$ given $\mathbf{x}$ is defined as
\begin{equation}
\widehat{p}(\mathbf{y} \gvn \mathbf{x}) = \frac{1}{k} \sum_{i=1}^k K_h\left[\rho(\mathbf{y},\mathbf{y}_{s_i(\mathbf{x})})\right],
\end{equation}
where $K_h$ is a normalized kernel (e.g., a Gaussian function) with bandwidth $h$, $\rho$ is a distance metric, and $s_i(\mathbf{x})$ is the index of the $i^{\rm th}$ nearest neighbor of $\mathbf{x}$.
It is essentially a smoother version of the histogram estimator proposed by \citet{cunha2009estimating} in that it approximates the density with a smooth continuous function rather than by binning.

We provide the \code{NNKCDE}\footnote{\url{https://github.com/tpospisi/nnkcde}} software in both \code{Python}\ and \code{R}\ \citep{NNKCDE_code}, with examples in Section~\ref{app:nnkcde}.
\code{NNKCDE}\ can accommodate any number of training examples.
However, selecting tuning parameters scales quadratically in $k$, the number of nearest neighbors, which prohibits using $k \gtrsim 10^{3}$.
Our implementation (\citealt[Appendix D]{izbicki2019abc}) is computationally more efficient than standard nearest-neighbor kernel smoothers.
For instance, we are able to efficiently evaluate the loss function in Equation~\ref{eq:estimated-cde-loss} on large validation samples by expressing the first integral in terms of convolutions of the kernel function;
these have a fast-to-compute closed solution for a Gaussian kernel.

\subsection{\code{RFCDE}}
\label{sec:RFCDE}

Random forests (RFs, \citealt{breiman2001random}) is one of the best off-the-shelf solutions for regression and classification problems.
It builds a large collection of decorrelated trees, where each tree is a data-based partition of the feature space.
The trees are then averaged to yield a prediction.
RFCDE, introduced by \cite{pospisil2018rfcde}, is an extension of random forests to conditional density estimation and multivariate responses.
Like \code{NNKCDE}, it computes a kernel density estimate of $\mathbf{y}$ but with nearest neighbor weightings defined by the location of the evaluation point $\mathbf{x}$ relative to the leaves in the random forest.
RFCDE inherits the advantages of random forests in that it can handle mixed-typed data.
It also does not require the user to specify distances or similarities between data points, and it has good performance while remaining relatively interpretable.

The main departure from other random forest algorithms is our criterion for feature space partitioning decisions.
In regression contexts, the splitting variable and split point are typically chosen so as to minimize the mean-squared-error loss.
In classification contexts, the splits are typically chosen so as to minimize a classification error.
Existing random forest density estimation methods such as quantile regression forests by \cite{meinshausen2006quantile} and the TPZ algorithm by \cite{carrasco2013tpz} use the same tree structure as regression and classification random forests, respectively.
\code{RFCDE}, however, builds trees that minimize the CDE loss (see Equation~\ref{eq:estimated-cde-loss}), allowing the forest to adapt to structures in the conditional density; hence overcoming some of the limitations of the usual regression approach for data with heteroskedasticity and multimodality.
In addition, \code{RFCDE}\ does not require discretizing the response as in TPZ, thereby providing more accurate results at a lower cost for continuous responses, especially in the case of multivariate continuous responses where binning is problematic.
See \cite{pospisil2018rfcde} for further examples and comparisons.

Another unique feature of \code{RFCDE}\ is that it can handle multivariate responses with joint densities by applying a weighted kernel smoother to $\mathbf{y}$.
This added feature enables analysis of complex conditional distributions that describe relationships between multiple responses and features, or equivalently between multiple parameters and observables in an LFI setting.
Like quantile regression forests, the \code{RFCDE}\ algorithm takes advantage of the fact that random forests can be viewed as a form of adaptive nearest-neighbor method with the aggregated tree structures determining a weighting scheme.
This weighting scheme can then be used to estimate the conditional density $\pr{\mathbf{y} \gvn \mathbf{x}}$, as well as the conditional mean and quantiles, as in quantile regression forests (but for CDE-optimized trees).
As mentioned above, \code{RFCDE}\ computes the latter density by a weighted kernel density estimate (KDE) in $\mathbf{y}$ using training points near the evaluation point $\mathbf{x}$.
These distances are effectively defined by how often a training point $\mathbf{x}_i$ belongs to the same leaf node as $\mathbf{x}$ in the forest (see  Equation 1 in \citealt{pospisil2018rfcde} for details).

Despite the increased complexity of our CDE trees, \code{RFCDE}\ still scales to large data sets because of an efficient computation of splits via orthogonal series.
Moreover, \code{RFCDE}\ extends the density estimates on new $\mathbf{x}$ to the multivariate case through the use of multivariate kernel density estimators \citep{epanechnikov1969non}.
In both the univariate and multivariate cases, bandwidth selection can be handled by either plug-in estimators or by tuning using a density estimation loss.

For ease of use in the statistics and astronomy communities, we provide \code{RFCDE}\footnote{\url{https://github.com/tpospisi/RFCDE}} in both \code{Python}\ and \code{R}, which call a common \code{C++} implementation of the core training functions that can easily be wrapped in other languages \citep{RFCDE_code}.

\textbf{Remark:}
Estimating a CDE loss is an inherently harder task than calculating the mean squared error (MSE) in regression.
As a consequence, \code{RFCDE}\ might not provide meaningful tree splits in settings with a large number of noisy features.
In such settings, one may benefit from combining a regular random forest tree structure (optimized for regression) with a weighted kernel density estimator (for the density calculation).
See Section \ref{sec:functional} for an application to functional data along with software implementation.\footnote{Available at \url{https://github.com/Mr8ND/cdetools_applications/spec_z/}}

\newpage

\subsubsection{\code{fRFCDE}}
\label{sec:frfcde}

In addition the \code{RFCDE}\ package includes \code{fRFCDE}, a variant of \code{RFCDE}\ \citep{pospisil2019}, that can accommodate functional features $\mathbf{x}$ by partitioning in the  continuous domain of such features.
The spectral energy distribution (SED) of a galaxy is its energy as a function of continuous wavelength $\lambda$
of light; hence it can be viewed as functional data.
Another example of functional data is the shear correlation function of weak lensing, which measures the mean product of the shear at two points as a function of a continuous distance $r$ between those points. Similarly, any function of continuous time is an example of functional data.
Treating functional features (like spectra, correlation functions or images) as unordered multivariate vectors on a grid suffers from a curse of dimensionality.
As the resolution of the grid becomes finer the dimensionality of the data increases but little additional information is added, due to high correlation between nearby grid points.
\code{fRFCDE}\ adapts to this setting by partitioning the  domain of each functional feature (or curve) into intervals, and passing the mean values of the function in each interval as input to \code{RFCDE}.
Feature selection is then effectively done over regions of the domain rather than over single points.
More specifically, the partitioning in \code{fRFCDE}\ is governed by the
parameter
$\mu$ of a Poisson process, with each functional feature entering as a high-dimensional vector $\mathbf{x} = (x_1,\ldots,x_d)$.
Starting with the first element of the vector, we group the first
\(\operatorname{Poisson}(\mu)\) elements together.
We then repeat the procedure sequentially until we have assigned all $d$ elements into a group;
this effectively partitions the function domain into disjoint intervals $\left\{(l_{i}, h_{i})\right\}$.
The function mean values
or smoothed brightness measurements $\widetilde{x}_i \equiv \int_{l_{i}}^{h_{i}} f(\lambda) d\lambda$ of each interval are finally treated as new inputs to a standard  (vectorial) \code{RFCDE}\ tree.
The splitting of the smoothed predictors $\widetilde{x}_i$ is done independently for each tree in the forest.
Other steps of \code{fRFCDE}, such as the computation of variable importance, also proceed as in (vectorial) \code{RFCDE}\ but with the averaged values of a region as inputs.
As a result, \code{fRFCDE}\ has the capability of identifying the functional inputs and the regions in the input domain that are important for estimating the response $\mathbf{y}$.
Figure \ref{fig:rfcde_viz} shows schematically the differences and similarities in construction between standard \code{RFCDE}\ and its \code{fRFCDE}\ variant.

The \code{fRFCDE}\ method is as scalable as standard random forests, accommodating training sets on the order of $10^6$.
As the examples in Section~\ref{sec:functional} show, we can obtain substantial gains with a functional approach, both in terms of statistical performance (that is, CDE loss) and computational time.
In addition, the change in code syntax is minimal, as one only needs to pass the
Poisson $\mu$ parameter as the \texttt{flambda} argument during the forest initialization.
Examples in \code{Python}\ and \code{R}\ are provided in Appendix~\ref{app:frfcde}.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{figs/spectrum_viz_plot.pdf}%
    \includegraphics[width=0.5\textwidth]{figs/trees_viz_tibs_style.pdf}\\
    \includegraphics[width=0.5\textwidth]{figs/spectrum_viz_plot_rfcde.pdf}%
    \includegraphics[width=0.5\textwidth]{figs/trees_viz_tibs_style_rfcde.pdf}
    \caption{
    A schematic diagram of \code{RFCDE}\ (top row) and \code{fRFCDE}\ (bottom row) applied to a galaxy spectrum from \cite{astroML2012}.
    \emph{Top row}:
    \code{RFCDE}\ treats the intensity $x_i$ at each recorded wavelength $\lambda_i$ of the spectrum as a feature or ``input'' to the random forests algorithm --- the blue vertical dashed lines indicate every 100$^{\rm th}$ recorded wavelength.
    \code{RFCDE}\ then builds an ensemble of CDE trees, where each tree partitions the feature space according to the CDE loss, as illustrated in the top right figure for features $x_1$ and $x_2$.
    \emph{Bottom row}:
    \code{fRFCDE}\ instead groups nearby measurements together where the group divisions are defined by a Poisson process with parameter $\mu$ (vertical green dashed lines, left figure).
    The new smoothed features $\widetilde{x_i}$ are computed by integrating the intensity over the grouped wavelengths.
    A forest of CDE trees is then built using the same construction as in \code{RFCDE}\ but with the smoothed features as inputs (bottom right figure).
    }\label{fig:rfcde_viz}
\end{figure}

\subsection{\code{FlexCode}}
\label{sec:FlexCode}

Introduced by \cite{izbicki2017converting}, \code{FlexCode}\footnote{ \url{https://github.com/tpospisi/FlexCode} (\code{Python}) and \url{https://github.com/rizbicki/FlexCoDE} (\code{R})}
is a CDE method that uses a basis expansion for the univariate response $y$ and poses CDE as a series of univariate regression problems.
The main advantage of this method is its flexibility as any regression method can be applied towards CDE, enabling us to tailor our choice of regression  method to the intrinsic structure and type of data at hand.

More precisely, let $\left\{\phi_{j} (y) \right\}_j$ be an orthonormal basis like a Fourier or wavelet basis for
functions of $y \in \mathbb{R}$.
The key idea of \code{FlexCode}\ is to express the unknown density $\pr{y \gvn \mathbf{x}}$ as a basis expansion
\begin{equation}
\label{eq:flexcode}
  p(y \gvn \mathbf{x}) = \sum_{j} \beta_{j}(\mathbf{x}) \phi_{j}(y) .
\end{equation}
By the orthogonality property of the basis, the (unknown) expansion coefficients $\{\beta_{j}(\mathbf{x})\}_j$ are then just orthogonal projections of $p(y|\mathbf{x})$ onto the basis vectors.
We can estimate these coefficients using a training set of $(\mathbf{x},y)$ data by {\em regressing} the transformed response variables $\phi_j(y)$ on predictors $\mathbf{x}$ for every basis function $j$ (see \cite{izbicki2017converting} Equation 2.2, for details).
The number of basis function $n_{basis}$ is chosen by minimizing a CDE loss function on validation data.
The estimated density, $\sum_{j=1}^{n_{basis}} \widehat{\beta}_j(\mathbf{x}) \phi_j(y)$, may contain small spurious bumps induced by the Fourier approximation and may not integrate to one.
We remove such artifacts as described in \citet{izbicki2016nonparametric} by applying a thresholding parameter $\delta$ chosen via cross-validation.
\code{FlexCode}\ turns a challenging density estimation problem into a simpler regression problem, where we can choose any regression method that fits the problem at hand.

To provide a concrete example, for high-dimensional $\mathbf{x}$ (such as galaxy spectra) we can use methods such as sparse or spectral series (``manifold'') regression methods \citep{tibshirani1996regression,ravikumar2009sparse,lee2016spectral};
see Section~\ref{sec:functional} for an example with \code{FlexCode}\code{-Series}.
For multi-probe studies with mixed data types, we can use random forests regression \citep{breiman2001random}.
On the other hand, large-scale photometric galaxy surveys such as \data{LSST}\
require methods that can work with data from millions, if not billions, of galaxies.
\citet{desc_photoz} present the results of an initial study of the \data{LSST}\ Dark Energy Science Collaboration (\textsc{LSST-DESC}).
Their initial data challenge (``Photo-$z$\ DC 1'') compares the CDEs of a dozen photo-$z$\ codes run on simulations of \data{LSST}\ galaxy photometry catalogs in the presence of complete, correct, and representative training data.
\code{FlexZBoost}, a version of \code{FlexCode}\ based on the scalable gradient boosting regression technique by \cite{chen2016xgboost}, was entered into the data challenge because of the method's ability to scale to massive data.
In the DC1 analysis, \code{FlexZBoost}\ was among the the strongest performing codes according to established performance metrics of such PDFs and was
one of only two codes to beat the experimental control under a more discriminating metric, the CDE loss.

For massive surveys such as \data{LSST}, \code{FlexCode}\ also has another advantage compared to other CDE methods, namely its compact, lossless storage format.
\citet{juric_data_2017} establishes that \data{LSST}\ has allocated $\sim$100 floating point numbers to quantify the redshift of each galaxy.
As is shown in \citet{desc_photoz}, the myriad methods for deriving photo-$z$\ PDFs yield radically different results, motivating a desire to store the results of more than one algorithm in the absence of an obvious best choice.
For the photo-$z$\ PDFs of most codes, one may need to seek a clever storage parameterization to meet \data{LSST}'s constraints \citep{carrasco_kind_sparse_2014, malz_approximating_2018}, but \code{FlexCode}\ is virtually immune to this restriction.
Since \code{FlexCode}\ relies on a basis expansion, one only needs to store $n_{basis}$ coefficients per target density for a lossless compression of
the estimated PDF with no need for binning.
Indeed, for DC1, we can with \code{FlexZBoost}\ reconstruct our estimate $\widehat{p}(z \gvn \mathbf{x})$ at any resolution from estimates of the first 35 coefficients in a Fourier basis expansion.
In other words, \code{FlexZBoost}\ enables the creation and storage of high-resolution photo-$z$\ catalogs for several billion galaxies at no added cost.

Our public implementation of \code{FlexCode}
--- available in both \code{Python}\ \citep{Flexcode_code} and \code{R}\ \citep{Flexcode_codeR},  respectively --- cross-validates over regression tuning parameters (such as the number of nearest neighbors $k$ in \code{FlexCode}\code{-kNN}) and computes the \code{FlexCode}\ coefficients in parallel for further time savings.
The computational complexity of \code{FlexCode}\ will be the same as $n_{basis}$ parallel individual regressions.
In particular, the scalability of \code{FlexCode}\ is determined by the underlying regression method.
A scalable method like XGBoost leads to scalable \code{FlexCode}\ fitting.
In the \code{Python}\ version of the code, the user can choose between the following regression methods: XGBoost for \code{FlexZBoost}\, but also nearest neighbors, LASSO \citep{tibshirani1996regression}, and random forests regression \citep{breiman2001random}.
In the \code{R}\ version, the following choices are available:
XGBoost, nearest neighbors, LASSO, random forests, Nadaraya-Watson kernel smoothing \citep{nadaraya1964estimating,watson1964smooth}, sparse additive models \citep{ravikumar2009sparse}, and spectral series (``manifold'') regression \citep{lee2016spectral}.
In both implementations, the user may also use a custom regression method;
we illustrate how this can be done with a vignette in both packages.
For the \code{Python}\ implementation, the user can add any custom regression method following the \code{sklearn} API, i.e., with \texttt{fit} and \texttt{predict} methods.
Examples in both languages are presented in Appendix~\ref{app:flexcode}.

\subsection{\code{DeepCDE}}
\label{sec:DeepCDE}

Recently, neural networks have reemerged as a powerful tool for prediction settings where large amounts of representative training data are available;
see \cite{lecun2015deep} and \cite{goodfellow2016deep} for a full review.
Neural networks for CDE, such as Mixture Density Networks (MDNs; \citealt{bishop1994mixture}) and variational methods \citep{Tang2013varationalCDE, Sohn2015generativeCDE}, usually assume a Gaussian or some other parametric form of the conditional density.
MDNs have lately also been used for photometric redshift estimation \citep{d2018photometric, pasquet2019photometric} and for direct estimation of likelihoods and posteriors in cosmological parameter inference (see \citealt{alsing2019fast} and references within).

\code{DeepCDE}\footnote{\url{https://github.com/tpospisi/DeepCDE}} \citep{DeepCDE_code} takes a different, {\em fully nonparametric} approach to CDE.
It combines the advantages of basis expansions with the flexibility of neural network architectures, allowing for data types like image features and time-series data.
\code{DeepCDE}\ is based on the orthogonal series representation in \code{FlexCode}, given in Equation \ref{eq:flexcode}, but rather than relying on regression methods to estimate the expansion coefficients in Equation~\ref{eq:flexcode}, \code{DeepCDE}\ computes the coefficients $\left\{\beta_{i}(\mathbf{x})\right\}_{i=1}^{B}$ {\em jointly} with a  neural network that minimizes the CDE loss in Equation~\ref{eq:cde-loss}.
Indeed, one can show that for an orthogonal basis, the problem of minimizing this CDE loss is (asymptotically) equivalent to finding the best basis coefficients in \code{FlexCode}\ under mean squared error loss for the individual regressions;
see Appendix~\ref{app:cdeloss_equiv} for a proof.
The value of this result is that \code{DeepCDE}\ with a CDE loss directly connects prediction with uncertainty quantification, implying that one can leverage the state-of-the-art deep architectures for an application at hand toward uncertainty quantification for the same prediction setting.

From a neural network architecture perspective, \code{DeepCDE}\ only adds a linear output layer of coefficients for a series expansion of the density according to
\begin{equation}
 \label{eq::expansion_DeepCDE}

   \hat{p}(y \gvn \mathbf{x}) = \sum_{j=1}^{B} \widehat{\beta}_{j}(\mathbf{x}) \phi_{j}(y) ,
\end{equation}
where $\{\phi_{j}(y)\}_{j=1}^B$ is an orthogonal basis for functions of $y \in \mathbb{R}$.
Like \code{FlexCode}, we normalize and remove spurious bumps from the final density estimates according to the procedure in Section 2.2 of \cite{izbicki2016nonparametric}.

The greatest benefit of \code{DeepCDE}\ is that it can be implemented with both convolutional and recurrent neural network architectures, extending to both image and sequential data.
For most deep architectures, adding a linear layer represents a small modification, and a negligible increase in the number of parameters.
For instance, with the AlexNet architecture \citep{Krizhevsky2012Alexnet},
a widely used, relatively shallow convolutional neural network, adding a final layer with 30 coefficients for a cosine basis only adds $\sim 120,000$ extra parameters.
This represents a $0.1\%$ increase over the $62$ million already existing parameters, and hence a negligible increase in training and prediction time.
Moreover, the CDE loss for \code{DeepCDE}\ is especially easy to evaluate; see Appendix \ref{app:cdeloss_equiv} for details.

We provide both \code{TensorFlow} and \code{Pytorch} implementations of \code{DeepCDE}\ \citep{DeepCDE_code}.
We also include examples that shows how one can easily build \code{DeepCDE}\ on top of AlexNet \citep{Krizhevsky2012Alexnet};
in this case, for the task of estimating the probability distribution $\pr{y \gvn \mathbf{x}}$ of the correct orientation $y$ of a color image $\mathbf{x}$.
Note that \cite{Fischer2015ImageOE} use AlexNet for the corresponding regression task of predicting the orientation $y$ for a color image $\mathbf{x}$ but without quantifying the uncertainty in the predictions.

\section{How to Assess Method Performance}
\label{sec:assessment}

After fitting CDEs, it is important to assess the quality of our models of uncertainty.
For instance, after computing photo-$z$\ PDF estimates $\widehat{p}(z \gvn \mathbf{x})$ for some galaxies, one may ask whether these estimates accurately quantify the true uncertainty distributions $\pr{z\gvn \mathbf{x}}$.
Similarly, in the LFI task, a key question is whether an estimate of the posterior distribution, $\widehat{p}(\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}})$ of the cosmological parameters is close enough to the true posterior $\pr{\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}}}$ given the observations $\ensuremath{\mathbf{x}_{\text{obs}}}$.

We present two method-assessment tools suitable to different situations, which are complementary and can be performed simultaneously, with public implementations in the \code{cdetools}\footnote{\url{https://github.com/tpospisi/cdetools}} package in both \code{Python}\ and \code{R}\ \citep{cdetools_code0}.
First, we describe a CDE loss function that directly provides relative comparisons between conditional density estimators, such as the methods presented in Section \ref{sec:cde_methods} or, equivalently, between a set of models (for the same method) with different tuning parameters.
Second, we describe visual diagnostic tools, such as Probability Integral Transforms (PIT) and Highest Probability Density (HPD) plots, that can provide insights on the overall goodness-of-fit of a given estimator to observed data.

\subsection{CDE loss}
\label{sec:cde_loss}

Here we briefly review the CDE loss from \cite{izbicki2016nonparametric} for assessing conditional density estimators and discuss it in the context of the cosmology LFI case.

The goal of a loss function is to provide relative comparisons between different estimators, so that it is easy to directly choose the best fitted model among a list of candidates.
Given an estimate $\widehat{p}$ of $p$, we define the CDE loss by
\begin{equation}
  \label{eq:cde-loss}
  L(\widehat{p}, p) = \int \int \left[\widehat{p}(\mathbf{y} | \mathbf{x}) - p(\mathbf{y} | \mathbf{x})\right]^{2} d\mathbf{y} dP(\mathbf{x}) ,
\end{equation}
where $P(\mathbf{x})$ is the marginal distribution of the features $\mathbf{x}$.
This loss is the CDE analog to the standard mean squared error (MSE) in standard regression.
The weighting by the marginal distribution of the features emphasizes that errors in the estimation of $\mathbf{y}$ for unlikely features $\mathbf{x}$ are less important.
The CDE loss cannot be directly evaluated because it depends on the unknown true density $p(z | \mathbf{x})$. However, one can estimate the loss (up to a constant determined by the true $p$) by
\begin{equation}
  \label{eq:estimated-cde-loss}
  \hat{L}(\widehat{p}, p) = \frac{1}{n} \sum_{i=1}^{n} \int \widehat{p}(\mathbf{y} | \mathbf{x}_{i}^{te})^{2} d\mathbf{y} - \frac{2}{n} \sum_{i=1}^{n} \widehat{p}(\mathbf{y}_{i}^{te} | \mathbf{x}_{i}^{te}) ,
\end{equation}
where $\left\{(\mathbf{x}_{i}^{te}, \mathbf{y}_{i}^{te})\right\}_{i=1}^{n}$ represents our validation or test data, i.e., a held-out set not used to construct $\widehat{p}$.
In our implementation, the function \texttt{cde\_loss} returns the estimated CDE loss as well as an estimate of the standard deviation or the {\em standard error} (SE) of the estimated loss.

{\bf CDE loss for LFI.} In LFI settings, we use a slightly different version of the CDE loss in Eq.~\ref{eq:cde-loss}. Because the goal (in ABC) is to approximate the posterior density $p(\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}})$ at fixed $\mathbf{x}=\ensuremath{\mathbf{x}_{\text{obs}}}$, a natural evaluation metric is the integrated squared error loss
\begin{equation}
  \label{eq:abc-loss}
  \int \left[\hat{p}(\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}}) - p(\boldsymbol{\theta} \gvn \ensuremath{\mathbf{x}_{\text{obs}}})\right]^{2} d\boldsymbol{\theta}
\end{equation}
of the conditional density \emph{at $\ensuremath{\mathbf{x}_{\text{obs}}}$ only}.
Estimating this loss can however be tricky as only a single instance of data with $x=\ensuremath{\mathbf{x}_{\text{obs}}}$ is available in practice.
Hence, \cite{izbicki2019abc} approximates Equation \ref{eq:abc-loss} by computing the empirical loss $\widehat{L}(\widehat{p}, p)$ in Eq.~\ref{eq:estimated-cde-loss} over a restricted subset of the validation data that only includes the  $\mathbf{x}_{i}^{te}$ points that fall in an $\epsilon$-neighborhood of $\ensuremath{\mathbf{x}_{\text{obs}}}$, where $\epsilon$ is the tolerance of the ABC rejection algorithm.
The detailed analysis of this approach can be found in \cite{izbicki2019abc}.

\subsection{PIT and HPD diagnostics}
\label{sec:hpd}

The CDE loss function is a relative measure of performance that cannot address absolute goodness-of-fit.
To quantify overall goodness-of-fit, we examine how well calibrated an ensemble of conditional density estimators are on average, over validation or test data $\left\{(\mathbf{x}_{i}^{te}, \mathbf{y}_{i}^{te})\right\}_{i=1}^{n}$. For ease of notation, we will in this section denote $\mathbf{x}_{i}^{te}$ and  $\mathbf{y}_{i}^{te}$ for a generic $i$ by $\mathbf{x}_{val}$ and $\mathbf{y}_{val}$.

Given a true probability density $p(\mathbf{y} \gvn \mathbf{x}) = \gamma_{0}$ of a variable $\mathbf{y}$ conditioned on data $\mathbf{x}$, an estimated probability density $\widehat{p}(\mathbf{y} \gvn \mathbf{x}) = \gamma$ cannot be {\em well-calibrated} unless $\gamma \approx \gamma_{0}$.
Built on the same logic, the probability integral transform (PIT; \citealt{Polsterer})
\begin{equation}
PIT(\ensuremath{\mathbf{x}_{\text{val}}}, \ensuremath{y_{\text{val}}}) = \int_{-\infty}^{\ensuremath{y_{\text{val}}}} \hat{p}(y \gvn \ensuremath{\mathbf{x}_{\text{val}}}) dy
\end{equation}
assesses the calibration quality of an ensemble of CDEs for scalar $y$ representing the cumulative distribution function (CDF) of $\hat{p}(y \gvn \ensuremath{\mathbf{x}_{\text{val}}})$ evaluated at $y = \ensuremath{y_{\text{val}}}$;  this PIT value corresponds to the shaded area in Figure \ref{fig:pit_hpd_visualization}, left.
A statistically self-consistent population of densities has a uniform distribution of PIT values, and deviations from uniformity indicate inaccuracies of the estimated PDFs.
Overly broad CDEs manifest as under-representation of the lowest and highest PIT values, whereas overly narrow CDEs manifest as over-representation of the lowest and highest values.

However, PIT values do not easily generalize to multiple responses.
For instance, for a bivariate response $\mathbf{y} = (z, \eta)$, the quantity $\int_{-\infty}^{\ensuremath{z_{\text{val}}}} \int_{-\infty}^{\ensuremath{\eta_{\text{val}}}} \pr{z, \eta \gvn \ensuremath{\mathbf{x}_{\text{val}}}} dz d\eta$ is not in general uniformly distributed \citep{genest2001multivariate}.
An alternative statistic that easily generalizes to multivariate $\mathbf{y}$ is the highest probability density value (HPD; \citealt[Appendix A]{izbicki2017photo}):
\begin{equation} \label{eq:HPD}
\xi(\ensuremath{\mathbf{x}_{\text{val}}}, \ensuremath{\mathbf{y}_{\text{val}}}) = \int_{\textbf{y}:\hat{p}(\textbf{y} \gvn \ensuremath{\mathbf{x}_{\text{val}}}) \ge \hat{p}(\ensuremath{\mathbf{y}_{\text{val}}} \gvn \ensuremath{\mathbf{x}_{\text{val}}})} \hat{p}(\textbf{y} \gvn \ensuremath{\mathbf{x}_{\text{val}}}) d\textbf{y} .
\end{equation}
The HPD value is based on the definition of the {\em highest density region} (HDR, \citealt{hyndman1996computing}) of a random variable $\mathbf{y}$;
that is, the  subset of the sample space of $\mathbf{y}$ where all points in the region have a probability above a certain value.
The HDR of $\textbf{y} \gvn \ensuremath{\mathbf{x}_{\text{val}}}$ can be seen as a region estimate of  $\textbf{y}$ when $\mathbf{x}=\ensuremath{\mathbf{x}_{\text{val}}}$ is observed.
In words, the set $\{\textbf{y}:\hat{p}(\textbf{y} \gvn \ensuremath{\mathbf{x}_{\text{val}}}) \ge \hat{p}(\ensuremath{\mathbf{y}_{\text{val}}} \gvn \ensuremath{\mathbf{x}_{\text{val}}})\}$ is the smallest HDR containing the point $\ensuremath{\mathbf{y}_{\text{val}}}$, and the HPD value
is simply the probability
of such a region.
Figure \ref{fig:pit_hpd_visualization}, right, shows a schematic diagram of the HPD value (green shaded area) and HDR region (highlighted segments on the y-axis) for the estimated density $\hat{p}(\textbf{y} \gvn \ensuremath{\mathbf{x}_{\text{val}}})$.
The HPD value $\xi(\ensuremath{\mathbf{x}_{\text{obs}}}, \ensuremath{\mathbf{y}_{\text{val}}})$ can also be viewed as a measure of how plausible $\ensuremath{\mathbf{y}_{\text{val}}}$ is according to $\hat p (\mathbf{y} | \ensuremath{\mathbf{x}_{\text{val}}})$ and is directly related to the Bayesian analog of p-values or {\em the e-value} \citep{pereira1999}.
One can show \citep{harrison2015validation} that even for multivariate $\mathbf{y}$, the
HPD values for validation data  follow a $U(0,1)$  distribution if the CDEs are well-calibrated on the population level. Thus, these values can also be used for assessing the fit of conditional densities in the same way as PIT values.
In our implementation the functions \texttt{pit\_coverage(cde\_test, y\_grid, y\_test)} and \texttt{hpd\_coverage(cde\_test, y\_grid, y\_test)}, respectively, calculate PIT and HPD values for CDE estimates \texttt{cde\_test} using the grid \texttt{y\_grid} with observed values \texttt{y\_test}.

The PIT and HPD are not without their limitations, however, as demonstrated in the control case of \citet{desc_photoz} and Figure~\ref{fig:setting1_diagnostics} of Section~\ref{sec:univ_response_mult_data} here.
Because the PIT and HPD values can be uniformly distributed even if $\pr{\mathbf{y} \gvn \mathbf{x}}$ is not well estimated, they must be used in conjunction with loss functions for method assessment.
A popular way of visualizing PIT and HPD diagnostics for the entire population is through   {\em  probability-probability plots} or P-P plots of the empirical distribution of the  (PIT or HPD) statistic  versus its distribution under the hypothesis that $\hat p(\mathbf{y} | \mathbf{x})= p(\mathbf{y} | \mathbf{x})$; henceforth, we will refer to the latter  Uniform(0,1) distribution as the ``theoretical'' distribution of PIT or HPD.  An ideal P-P plot has all points close to the identity line where the ``empirical'' and ``theoretical'' distributions are the same. Note that HPD P-P plots, in particular, are valuable calibration tools if our goal is to calibrate the estimated densities so that the computed  predictive regions have the right coverage.

Figure \ref{fig:teddy-pvals} (c)-(d) contains example P-P plots for the PIT and HPD values of photo-$z$\ estimators, where the empirical (observed) distribution is close to the theoretical (ideal) distribution.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{figs/pit_plot_viz_val.pdf}%
    \includegraphics[width=0.5\textwidth]{figs/hpd_plot_viz_val.pdf}
    \caption{
    Schematic diagram of the construction of the Probability Integral Transform (PIT, left) and  the Highest Probability Density (HPD, right) values for  the estimated density $\widehat{p}(y|\mathbf{x})$ at $\mathbf{x}=\ensuremath{\mathbf{x}_{\text{val}}}$, where  $\ensuremath{y_{\text{val}}}$ is the response at $x=\ensuremath{\mathbf{x}_{\text{val}}}$. In the  plot to the right, the highlighted segments on the $y$-axis form the so-called highest density region (HDR) of $y|\ensuremath{\mathbf{x}_{\text{val}}}$. The PIT and HPD values correspond to the area of the tail versus highest density region, respectively, of the estimate; here indicated by the blue versus green shaded areas.
    }\label{fig:pit_hpd_visualization}
\end{figure}

\newpage

\section{Applications in Astronomy}
\label{sec:examples}

We demonstrate\footnote{Code for these examples is publicly available at \url{https://github.com/Mr8ND/cdetools_applications}.} the breadth of our CDE methods in three different astronomical use cases.

\begin{enumerate}
  \item
  {\bf Photo-$z$\ estimation: univariate response with multivariate input.}
  This is the standard prediction setting in which all four methods apply.
  In Section \ref{sec:univ_response_mult_data}, we apply the methods to the \data{Teddy}\ photometric redshift data \citep{beck2017realistic}, and illustrate the need for loss functions to properly assess PDF estimates of redshift $z$ given photometric colors $\mathbf{x}$.
  \item
  {\bf Likelihood-free cosmological inference: multivariate response.}
  For multiple response components we want to model the often complicated dependencies between these components;
  this is in contrast to approaches which model each component separately, implicitly introducing an assumption of conditional independence.
  In Section \ref{sec:mult_response_mult_data}, we use an example of LFI for simulated weak lensing shear data to show how \code{NNKCDE}\ and \code{RFCDE}\ can capture more challenging bivariate distributions with curved structures;
  in this toy example $\mathbf{y}$ represents cosmological parameters $(\Omega_M, \sigma_8)$ in the $\Lambda$CDM-model, and $\mathbf{x}$ represents (coarsely binned) weak lensing shear correlation functions.
  \item
  {\bf Spec-$z$ estimation: functional input.}
  Standard prediction methods, such as random forests, do not typically fare well with functional features, simply treating them as unordered vectorial data and ignoring the functional structure.
  However, often there are substantial benefits to explicitly taking advantage of that structure as in \code{fRFCDE}.
  In Section \ref{sec:functional}, we compare the performance of a vectorial implementation of \code{RFCDE}\ with \code{fRFCDE}\ and \code{FlexCode}\code{-Series} for a spectroscopic sample from the Sloan Digital Sky Survey (\data{SDSS}; \citealt{alam2015eleventh}). The input $\mathbf{x}$ is here a high-resolution spectrum of a galaxy, and the response is the galaxy's redshift  $z$.
\end{enumerate}

Throughout this section we will report CDE loss mean and standard error for each method, using $95\%$ Gaussian confidence intervals for performance comparison.
Since the CDE loss is an empirical mean and the test size is reasonably large in all the examples, the validity of this approximation is guaranteed by the central limit theorem.

\subsection{Photo-$z$\ estimation: Univariate response with multivariate input}
\label{sec:univ_response_mult_data}

Here we estimate photo-$z$\ posterior PDFs $\pr{z \gvn \mathbf{x}}$ using representative training and test data (Samples A and B, respectively) from the \data{Teddy}\ catalog by \cite{beck2017realistic}.\footnote{Data available at \url{https://github.com/COINtoolbox/photoz_catalogues}} These data include 74309 training and 74557 test observations.
Each observation $\mathbf{x}$ has five features:
the magnitude of the $r$-band and the pairwise differences or ``colors'' $u-g$, $g-r$, $r-i$, and $i-z$.

Among the chief sources of uncertainty affecting photo-$z$s estimated by ML techniques are the incompleteness and nonrepresentativity of training sets, defined by the mismatch in the distributions of training and test data in $z$ and $\mathbf{x}$, which may be extreme to the point of not guaranteeing mutual coverage.
Realistically modeling incompleteness is highly challenging, requiring both simulations of SEDs and of the observational conditions of a given survey, which is outside the scope of this work.
Accounting for redshift incompletenss is not a lost cause and may be accomplished by extrapolating outside of the training range by abandoning standard instance-based ML algorithms (see e.g., \cite{leistedt2016_gp}).
Certain types of selection bias, known as covariate shift, can also be corrected by importance weights in the CDE loss \citep{izbicki2017photo, freeman2017unified};
see Appendix~\ref{sec:cov_shift_photoz} for details and code.
However, for simplicity, in this paper we consider only representative training sets with no disparity in color-space coverage, putting this demonstration on equal footing with all previous comparisons of photo-$z$ PDF methods.

We fit \code{NNKCDE}, \code{RFCDE}, \code{DeepCDE}, and {\code{FlexZBoost}}
to the data.
In addition we fit an \code{RFCDE}\ model, ``\code{RFCDE}\code{-Limited},'' restricted to the first three of the five features, as a toy model which fails to extract some information from the features, allowing us to showcase the difference between PIT or HPD diagnostics and the CDE loss function.
For \code{DeepCDE}\ we use a three-layer deep neural network with linear layers and reLu activations with $25$ neurons per layers, trained for $10,000$ epochs with Adam \citep{Kingma2014AdamAM}.
As our goal is to showcase its applicability we do not optimize the neural architecture (e.g., number of layers, number of neurons per layer, activation functions) nor the learning parameters (i.e., number of epochs, learning rate, momentum).
To illustrate our validation methods, we also include the marginal distribution $\widehat{p}(z) = \frac{1}{n}\sum_{i=1}^n \widehat{p}(z | \mathbf{x}_i)$ as an estimate of individual photo-$z$\ distributions $p(z | \mathbf{x})$.
This estimate will be the same regardless of $\mathbf{x}$.

\begin{figure}[htbp]
\centering
\includegraphics[width=1\linewidth]{figs/teddy/cdf-hpd-pvals-horizontal.pdf}

\vspace{0.75cm}

\includegraphics[width=1\linewidth]{figs/teddy/pit-hpd-pplot-distr.pdf}
\caption{
(a)-(b) Histograms of PIT and HPD values for \data{Teddy}\ photo-$z$\ data.
Both versions of \code{RFCDE}\ as well as the marginal distribution have a uniform distribution indicating  that the conditional density estimates are well-calibrated on average.
\code{FlexZBoost}\ exhibits an overrepresentation of 0 and 1 values which indicates overly narrow CDEs.
(c)-(d) Probability-Probability (P-P) plots of HPD and PIT values. Both sets of values are computed over   data in the test set  and their empirical distribution is plotted against the uniform $U(0,1)$ distribution; i.e., the ``theoretical distribution'' of the HPD and PIT values when $\hat{p}(z|\mathbf{x})=p(z|\mathbf{x})$. If the estimated CDEs are well calibrated, the  empirical and theoretical distributions should coincide and all points should be close to the identity line.
As in the top panel, the P-P plots indicate a good fit for all methods including the clearly misspecified ``Marginal'' model.
 \label{fig:setting1_diagnostics}}\label{fig:teddy-pvals}
\end{figure}

\begin{figure}[htbp]
\begin{minipage}[t]{0.5\textwidth}
\vspace{0pt}
\centering
\includegraphics[width=1\textwidth]{figs/teddy/teddy_performance_cdeplot_large.pdf}
\end{minipage}%
\begin{minipage}[t]{0.4\textwidth}
\vspace{10pt}
{\footnotesize
\begin{tabular}{llr}
Method & CDE Loss $\pm$ SE & Storage (single CDE)\\
\hline
Marginal & -3.192 $\pm$ 0.007 & 200 floats (1.6 KB)\\
RFCDE-Limited & -11.038 $\pm$ 0.047 & 200 floats (1.6 KB)\\
NNKCDE & -12.139 $\pm$ 0.043 & 200 floats (1.6 KB)\\
FlexZBoost & -12.272 $\pm$ 0.044 & 30 coeffs (0.24 KB)\\
DeepCDE & -12.821 $\pm$ 0.04 & 31 coeffs (0.25 KB)\\
RFCDE & -13.108 $\pm$ 0.052 & 200 floats (1.6 KB)
\end{tabular}
}
\end{minipage}
\caption{Method comparison via the CDE loss of Equation~\ref{eq:estimated-cde-loss} with estimated standard error (SE) and the storage space for each galaxy's photo-$z$\ CDE.
As \code{FlexZBoost}\ and \code{DeepCDE}\ are basis expansion methods, we need only to store the estimated coefficients for a lossless compression of the CDEs;
the other CDEs are discretized to 200 bins.}\label{tab:teddy-cde-loss}
\end{figure}

Figure~\ref{tab:teddy-cde-loss} presents the estimated CDE loss of Equation~\ref{eq:estimated-cde-loss} with estimated standard error (SE), as well as the storage space for each galaxy's photo-$z$\ CDE.
Figure \ref{fig:teddy-pvals} shows that the different models, including the clearly misspecified ``Marginal'' model, achieve comparable performance on goodness-of-fit diagnostics.
However, Figure~\ref{tab:teddy-cde-loss} shows the discriminatory power of the CDE loss function, which distinguishes the methods from one another.
This emphasizes the need for method comparison through loss functions in addition to goodness-of-fit diagnostics.
We note that the ranking in CDE loss also correlates roughly with the quality of the point estimates, as shown in Appendix \ref{sec:app_pointestimates}.
The \data{LSST} -DESC PZ DC1 paper \citep{desc_photoz} draws similar conclusions from a comprehensive photo-$z$ code comparison where the ``Marginal'' model (there referred to as \texttt{trainZ} photo-$z$ PDF estimator, the experimental control) outperformed all codes when using traditional metrics for assessing photo-$z$ PDF acccuracy.
Indeed, of the metrics considered in DC1, the CDE loss was the only metric that could appropriately penalize the pathological \texttt{trainZ}.

\subsection{Likelihood-free cosmological inference: Multivariate response}
\label{sec:mult_response_mult_data}

To showcase the ability to target joint distributions, we apply \code{NNKCDE}\ and \code{RFCDE}\ to the problem of estimating multivariate posteriors $p(\boldsymbol{\theta} | \mathbf{x}_{\rm obs})$ of the cosmological parameters in a likelihood-free setting via ABC-CDE~\citep{izbicki2019abc}.

ABC is an approach to parameter inference in settings where the likelihood is not tractable but we have a forward model that can simulate data $\mathbf{x}$ under fixed parameter settings $\boldsymbol{\theta}$.
The simplest form of ABC is the ABC rejection sampling algorithm, where a set of parameters is first drawn from a prior distribution.
The simulated data $\mathbf{x}$ is accepted with tolerance $\varepsilon \geq 0$ if $d(\mathbf{x},\mathbf{x}_{\rm obs}) \leq \varepsilon$ for some distance metric $d$ (e.g., the Euclidean distance) that measures the discrepancy between the simulated data $\mathbf{x}$ and observed data $\mathbf{x}_{\rm obs}$.
The outcome of the ABC rejection algorithm for small enough $\epsilon$ is a sample of parameter values approximately distributed according to the desired posterior distribution $p(\boldsymbol{\theta} | \mathbf{x}_{\rm obs})$.

The basic idea of ABC-CDE is to improve the ABC estimate --- and hence reduce the number of required simulations --- by using the ABC sample/output as input to a CDE method tuned with the CDE loss restricted to a neighborhood of $\ensuremath{\mathbf{x}_{\text{obs}}}$ defined by the tolerance $\epsilon$ (\cite{izbicki2019abc}, Equation 3).
Hence, in ABC-CDE, our CDE method can be seen as a post-adjustment method:
it returns an estimate $\widehat{p}(\boldsymbol{\theta} | \mathbf{x})$ which we evaluate at the point $\mathbf{x}=\mathbf{x}_{\rm obs}$ to obtain a more accurate approximation of $p(\boldsymbol{\theta} | \mathbf{x}_{\rm obs})$.
This could also be beneficial in an active learning setting where the posterior distribution is used to identify relevant regions of the parameter space (e.g., \citep{lueckmann2019likelihood, papamakarios2019sequential, alsing2019fast}).

In this example, we consider the problem of cosmological parameter inference via {\em cosmic shear}, caused by weak gravitational lensing inducing a distortion in images of distant galaxies.
The size and direction of the distortion is directly related to the size and shape of the matter distribution along the line of sight, which varies across the universe.
We use shear correlation functions to constrain the dark matter density \(\Omega_{M}\) and matter power spectrum normalization \(\sigma_{8}\) parameters of the \(\Lambda\)CDM cosmological model, which predicts the properties and evolution of the large scale structure of the matter distribution.
For further background see \cite{hoekstra2008weak}, \cite{munshi2008cosmology} and \citet{mandelbaum2017weak}.
Here we use the \code{GalSim}\footnote{\url{https://github.com/GalSim-developers/GalSim}} toolkit \citep{rowe2015galsim} to generate simplified galaxy shears distributed according to a Gaussian random field determined by $(\Omega_{M}, \sigma_{8})$.
The binned shear correlation functions serve as our input data or summary statistics $x$.
For the inference, we assume uniform priors \(\Omega_{M} \sim U(0.1, 0.8)\) and \(\sigma_{8} \sim U(0.5, 1.0)\) and fix \(h = 0.7\), \(\Omega_{b} = 0.045\), \(z = 0.7\).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figs/lensing_posterior_with_true_point.pdf}
\caption{\label{fig:orge983804}
 Cosmological parameter inference in an LFI setting with simulated weak lensing data.

 The top row shows the estimated bivariate distribution of $\Omega_M$ and $\sigma_8$ for ABC rejection sampling at different acceptance rates (20\%, 50\%, and 100\%).
 The middle row shows the estimated posterior densities after applying \code{NNKCDE}\ to the ABC sample, and the bottom row when applying \code{RFCDE}.\
 Both \code{NNKCDE}\ and \code{RFCDE}\ tuned with a CDE loss improve on ABC across acceptance rate levels, i.e., they provide approximate posteriors that are more concentrated around the true observed parameter values and the degeneracy curve on which the data are indistinguishable (shown here as a black diamond and dashed line, respectively).
 We even see some structure at an ABC acceptance rate of 1 (right column);
 that is, at an ABC threshold of $\epsilon \rightarrow \infty$ for which the entire sample is accepted by ABC and passed to our CDE code.
 HPD values were not included as they fail to distinguish between conditional and marginal distribution, as mentioned above.}
\end{figure}

The top row of Figure~\ref{fig:orge983804} shows the estimated bivariate posterior distribution of $\boldsymbol{\theta}=(\Omega_{M},\sigma_{8})$ from ABC rejection sampling only, at varying acceptance rates (20\%, 50\%, and 100\%).
An acceptance rate of 100\% just returns the (uniform) ABC prior distribution, whereas the ABC posteriors for smaller acceptance rates (that is, smaller values of $\varepsilon$) concentrate around the parameter degeneracy curve (shown as a dashed line) on which the data are indistinguishable.
The second and third rows show the estimated posteriors when we, respectively, improve the initial ABC estimate by applying \code{NNKCDE}\ and \code{RFCDE}\ tuned with our CDE surrogate loss.
A
result that is apparent from the figure
is that \code{NNKCDE}\ and \code{RFCDE}\ fitted with a CDE loss are able to capture the degeneracy curve at a larger acceptance rate, that is, for a smaller number of simulations, than when using ABC only.

\begin{table}[htbp]
  \label{tab:lfi-experiment-cdeloss}
  \caption{\small Performance of ABC, \code{NNKCDE}, \code{RFCDE}\ in LFI settings with simulated weak lensing data in terms of the surrogate CDE loss.}
  \centering
  \begin{tabular}{lrrr}
    \hline
    & \multicolumn{3}{c}{CDE Loss $\pm$ SE ($\times 10^{-5}$)} \\ \hline
    Method \textbackslash Acceptance Rate & 20\% & 50\% & 100\% \\ \hline
     ABC & -0.686 $\pm$ 0.009 & -0.392 $\pm$ 0.004 & -0.227 $\pm$ 0.001 \\
     \code{NNKCDE}\ & -1.652 $\pm$ 0.022 & -1.199 $\pm$ 0.016 & -0.844 $\pm$ 0.010\\
     \code{RFCDE}\ & -4.129 $\pm$ 0.063 & -3.698 $\pm$ 0.064 & -2.817 $\pm$ 0.055 \\
  \end{tabular}
\end{table}

We also have a direct measure of performance via the CDE loss and can adapt to different types of data by leveraging different CDE codes.
Table \ref{tab:lfi-experiment-cdeloss} shows how the methods compare in terms of the surrogate loss.
While decreasing the acceptance rate benefits all methods, it is clear that CDE-based approaches have better performance in all cases.

\subsection{Spec-$z$ estimation: Functional input}
\label{sec:functional}
In this example, we compare CDE methods in the context of spectroscopic redshift prediction using 2812 high-resolution spectra from the Sloan Digital Sky Survey (\textsc{SDSS}) Data Release 6 (preprocessed with the cuts of \citealt{Richards:EtAl:2009}), corresponding to features $\mathbf{x}$ of flux measurements at $d=3501$ wavelengths.
The high-resolution spectra $\mathbf{x}$ can be seen as functional inputs on a continuum of wavelength values.
Spectroscopic redshifts (or redshifts $z$ predicted from spectra $\mathbf{x}$) tend to be both accurate and precise, so the density $p(z \gvn \mathbf{x})$ is well-approximated by a delta function at the true redshift.
For the purposes of illustrating the use of the CDE codes, we define a noisified redshift $z_i= z_i^{\rm SDSS}+\epsilon_i$, where $\epsilon_i$ are independent and identically distributed variables drawn from a normal distribution $N(0,0 .02)$ and $z_i^{\rm SDSS}$ is the true redshift of galaxy $i$ provided by \textsc{SDSS}.
Thus the conditional density $p(z_i \gvn \mathbf{x}_i)$ of this example is a Gaussian distribution with mean $z_i^{\rm SDSS}$ and variance 0.02.

We compare \code{fRFCDE}, a ``Functional'' adaptation of \code{RFCDE}, with a standard ``Vector'' implementation of \code{RFCDE}, which treats functional data as a vector.
For completeness, we also compare against standard regression random forest combined with KDE, as well as \code{FlexCode}\code{-Spec}, an extension of \code{FlexCode}\ with a Spectral Series regression~\citep{Richards:EtAl:2009,Freeman, lee2016spectral} for determining the expansion coefficients.
We train on 2000 galaxies and test on the remaining galaxies.
The \code{RFCDE}\ and the \code{fRFCDE}\ trees are both trained with $n_{trees} = 1000$, $n_{basis}  = 31$, and bandwidths chosen by plug-in estimators.
We use \(\lambda=50\) for the \code{fRFCDE}\ rate parameter of the Poisson process that defines the variable groupings.

\begin{table}[htbp]
  \label{tab:frfcde-experiment}
  \caption{\small Performance of \code{RFCDE}\ for functional data.
  We achieve both a lower CDE loss and computational time for Functional \code{RFCDE} (as compared to Vector \code{RFCDE}) by leveraging the functional nature of the data.}
  \centering
  \begin{tabular}{lrr}
    Method & Train Time (in sec) & CDE Loss $\pm$ SE\\
    \hline
     Functional \code{RFCDE}\ & 24.89 & -3.38 $\pm$ 0.155\\
     Vector \code{RFCDE}\ & 41.60 & -2.52 $\pm$ 0.104\\
     Regression RF + KDE & 50.63 & -3.42 $\pm$ 0.120 \\
     \code{FlexCode}\code{-Spec} & 4017.93 & -3.53 $\pm$ 0.136 \\
  \end{tabular}
\end{table}

Table \ref{tab:frfcde-experiment} contains the CDE loss and train time for both the vector-based and functional-based \code{RFCDE}\ models on the \textsc{SDSS} data, as well as for \code{FlexCode}\code{-Spec}.
We obtain substantial gains when incorporating functional features both in terms of CDE loss and computational time.
The computational gains are attributed to requiring fewer searches for each split point as the default value of $m_{try} = \sqrt{d}$ is reduced.
As anticipated, vector-based \code{RFCDE}\ underperforms with functional data due to the tree splits on CDE loss struggling to pick up signals in the data and returning almost the same conditional density regardless of the input.
In contrast, splitting on mean squared error is an easier task and we include the results of regular regression random forest and KDE, which are comparable to the results of \code{fRFCDE}.
Whereas random forests rely on variable selection (either individual variables as in vectorial \code{RFCDE}\ or grouped together as in \code{fRFCDE}) for dimensional reduction, \code{FlexCode}\code{-Spec} is based on the Spectral Series mechanism for dimension reduction that finds (potentially) nonlinear, sparse structure in the data distribution~\citep{izbicki2016nonparametric}.
Both types of dimension reduction are reasonable for spec-$z$ estimation, as particular wavelength locations or regions in the galaxy SED could carry information of the galaxy's true redshift;
\code{RFCDE}\ and \code{fRFCDE}\ are effective in finding such locations by variable selection.
Spectral Series on the other hand are able to recover low-dimensional manifold structure in the entire data ensemble;
as \citealt{Richards:EtAl:2009} show, the main direction of variation in SDSS spectra is directly related to the spectroscopic redshift.

\section{Conclusions}
This paper presents statistical tools and software for uncertainty quantification in complex regression and parameter inference tasks.
Given a set of features $\mathbf{x}$ with associated response $\mathbf{y}$, our methods extend the usual point prediction task of classification and regression to nonparametric estimation of the entire (conditional) probability density $p(\mathbf{y} \gvn \mathbf{x})$.
The described CDE methods are meant to handle a range of different data settings as outlined in Table~\ref{tab:method-strenghts}. This paper includes examples
of code usage in the contexts of photo-$z$\ estimation, likelihood-free inference for cosmology analysis, and spec-z estimation.
In addition, it provides tools for CDE method assessment and for choosing tuning parameters of the CDE methods in a principled way.

Our software includes four packages for CDE --- \code{NNKCDE}, \code{RFCDE}, \code{FlexCode}\ and \code{DeepCDE} --- each using a different machine learning algorithm, as well as a package for  model assessment, \code{cdetools}.  All packages are implemented in both \code{Python}\ and \code{R}, with no dependence on proprietary software, making them compatible with any operating system supporting either language (e.g., Windows, MacOS, Ubuntu, and any other Unix-based OS).
The code is provided in publicly available Github repositories, documented using \code{Python}\ and \code{R}\ function documentation standards, and equipped with \code{Travis CI}\footnote{\url{https://docs.travis-ci.com/}} automatic bug-tracking.
Finally, our software makes uncertainty quantification straightforward for those used to standard open-source machine learning \code{Python}\ packages:
\code{NNKCDE}, \code{RFCDE}, {\code{FlexCode}} share the \code{sklearn} API (with \texttt{fit} and \texttt{predict} methods), which makes our methods compatible with \code{sklearn} wrapper functions (for e.g., cross-validation and model ensembling).
\code{DeepCDE}\ has implementations for both \code{Tensorflow} and \code{Pytorch}, two of the most widely used deep learning frameworks, and can easily be combined with most existing network architectures.

\section*{Acknowledgements}

We would like to thank the two anonymous reviewers for their insightful comments that helped improve the manuscript.
ND, TP and ABL were partially supported by the National Science Foundation under Grant Nos.\ DMS1521786.
RI was supported by Conselho Nacional de Desenvolvimento Cient\'ifico e Tecnol\'ogico (grant number 306943/2017-4) and Funda\c{c}\~ao de Amparo \`a Pesquisa do Estado de S\~ao Paulo (grants number 2017/03363-8
and 2019/11321-9).
AIM acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research.
During the completion of this work, AIM was advised by David W. Hogg and supported in part by National Science Foundation grant AST-1517237.

\newpage


text:
Q: how can I add git submodule into git repo as normal directory? Let's see what it is

*

*Create two git repo sub & test

mkdir test sub
cd test && git init && touch README && git add README && git commit -m "initialize the git repo" && cd ..
cd sub && git init && touch README && git add README && git commit -m "initialize the sub git repo" && cd ..

*

*Move the sub repo into test

mv sub test
cd test
git add sub
git commit -m "add sub directory"

I want to treat them as one git repo and push them remotely, but now the files under sub directory can not be included ?
How can I achieve this in simple way like treat sub as normal directory ?
Use case for this
I try to add my jenkins data folder (JENKINS_HOME) into docker images using Dockerfile for demo. (ADD JEKINS_HOME /opt/jenkins)
JENKINS_HOME
Dockerfile

My jenkin has scriptler plugin which contains the git repo for its purpose. Then it exists in my docker image git repo like below

$ find jenkins-docker
./.git
./.git/.. (skipped
./Dockerfile
./JENKINS_HOME
./JENKINS_HOME/scriptler
./JENKINS_HOME/scriptler/scripts
./JENKINS_HOME/scriptler/scripts/.git
./JENKINS_HOME/scriptler/scripts/.git/... (skipped)
./JENKINS_HOME/scriptler/scripts/Sample.groovy
./JENKINS_HOME/... (skipped)
./README

A: Seems you cannot do that. The name .git is hard-coded in the source code: https://github.com/git/git/blob/fe9122a35213827348c521a16ffd0cf2652c4ac5/dir.c#L1260
Probably one way is to make a script which renames .git to something else and back before and after adding it into repo like
In working directory under scripts
mv .git hidden-git

In Dockerfile
RUN mv $JENKINS_HOME/scriptler/scripts/hidden-git
$JENKINS_HOME/scriptler/scripts/.git

Alternatively, probably it's possible to pass GIT_DIR environment variable into the plugin, so it could use another name.

A: git add sub
git commit -m "add sub directory"

I want to treat them as one git repo and push them remotely, but now the files under sub directory can not be included ?

They are not included because test sees repo sub as nested git repo, and records only its gitlink, or SHA1, and not its url as it would have if sub had been added as a submodule.
You would need to push sub to a remote url first, and then add it as submodule for test to see sub files.
 cd test
 git submodule add -- /url/to/sub

Or you would need to use a subtree
cd test
git subtree --prefix sub /url/to/sub master --squash

redpajama_set_name: RedPajamaStackExchange


text:
Высокое — упразднённый населённый пункт в Костромской области, Нерехтский район, находившийся на территории Воскресенского сельского поселения.

География
В 1,5 км к востоку протекает река Чёрная. На западе от селения лежит лес. Так же лес располагается на северо-западе. На северо-востоке торфяное поле. На юге от села небольшой лес за которым следуют поля. За ними же стоит соседняя деревня Челпаново. Само Высокое построено на возвышенности, от этого и идёт его название.
Сегодня стоит экологическая проблема из-за организации песчаного карьера на западе Залесья.

Находится в 1,5 км к северо-востоку от посёлка Космынино.

Ссылки

Примечания

Исчезнувшие населённые пункты Костромской области
Урочища Костромской области
Охраняемые природные территории Костромской области

redpajama_set_name: RedPajamaWikipedia


West Plum Studio is an ethical design studio based in York. Inspired by nature, architecture and other cultures Sophie creates surface pattern designs with triangles - she's a little obsessed. She adds these designs to cards and on organic fabric and tells us about how she's growing her business.
More than 13 MILLION people in the UK are living below the poverty line and one million people are fed by food banks, yet we're the 7th richest country in the world. Listen to hear how Caroline tackles this through Foodinate.


Heat marshmallows and butter in 3-quart saucepan over low heat, stirring constantly, until marshmallows are melted and mixture is smooth.
Stir in cereal, raisins and nuts until evenly coated.
Insert wooden stick into middle of each ball.
Let stand for about 30 minutes or until firm.
Store balls wrapped in plastic wrap.


Reina Dracon is the daughter of the dragonslayer who turned into Acnologia.
Reina was distraught when her father turned into a dragon. She was a dragonslayer at this time, and used the same type of dragonslayer magic as her father. When her village was ravaged by Acnologia, she put herself into a sleep which lasted many eras. She woke up in the year 788 as a 17-year old.
When Reina woke up, she searched desperately for new of her father, and was terrified (and, in an odd way, slightly relieved) that her father was alive. Reina started training relentlessly, concinved that she would rid the world of her father.
Reina realized that although she was strong, would need others to help her defeat Acnologia. She evantually joined Sabertooth and became friends with Sting and Rogue.
A year later, Reina started dating Rogue.
Reina is dark, brooding, and quiet. Many people are intimidated by her. She often bottles her feelings up inside, and usually prefers to be alone.
Reina uses the same type of Dragonslayer magic as her father. She is often referred to as the Black Dragonslayer.


text:
Oregon Snow Kirby Wright Micro-Fiction
pages Oregon Snow
by Kirby Wright
Published in Issue No. 274 ~ March, 2020
photo_camera by Photo by Jason Blackeye on Unsplash
A truck with Oregon plates arrives in San Diego at dusk. When the driver slides open the cargo door, his Christmas trees are coated with lumps of snow and sheets of ice. There's so much white he imagines an artificial forest smothered in flocking spray.
"How'd this happen?" asks Home Depot man.
The driver tucks his hands in the kangaroo pocket of his sweatshirt. "Forgot to close that damn vent," he mutters.
The driver climbs on the bumper, grabs a chain hanging from the roof, and swings up. He shakes branches—white tumbles onto the steel bed. He digs in a blade and shovels snow and ice into the parking lot. A mogul forms. Kids gather. Most are boys but there are twin sisters. They grab the cold and shape without gloves. Snowballs sail through the half-light and splatter on asphalt. A few fathers join in. One father under-hands a snowball that bounces off his wife's head.
"Wanna play rough, eh?" she giggles and retaliates with ice.
The driver enjoys this unexpected winter. A boy leaps over the mogul and others follow. The snow and ice don't seem to melt. The driver's happy his vent was open as he churned through the passes north of Eugene. The evergreen smell reminds him of home. He feels like a boy again.
account_box More About Kirby Wright
Kirby Wright was born and raised in Honolulu, Hawaii. He is a graduate of Punahou School in Honolulu and the University of California at San Diego. He received his MFA in Creative Writing from San Francisco State University. Wright has been nominated for two Pushcart Prizes and is a past recipient of the Jodi Stutz Memorial Prize in Poetry, the Ann Fields Poetry Prize, the Academy of American Poets Award, the Robert Browning Award for Dramatic Monologue, and Arts Council Silicon Valley Fellowships in Poetry and The Novel. BEFORE THE CITY, his first poetry collection, took First Place at the 2003 San Diego Book Awards. Wright is also the author of the companion novels PUNAHOU BLUES and MOLOKA'I NUI AHINA, both set in Hawaii. He was a Visiting Fellow at the 2009 International Writers Conference in Hong Kong, where he represented the Pacific Rim region of Hawaii. He was also a Visiting Writer at the 2010 Martha's Vineyard Residency in Edgartown, Mass., and the 2011 Artist in Residence at Milkwood International, Czech Republic. His futuristic novel THE END, MY FRIEND was published in 2013. He published SQUARE DANCING AT THE ASYLUM, a collection of flash fiction, in 2014.

redpajama_set_name: RedPajamaCommonCrawl


text:
Translational Infectious Disease Research Program
Antimicrobial resistance is considered by the National Foundation for Infectious Diseases to be the number one long term infectious disease problem in the U.S. Bacteria such as Pseudomonas, Klebsiella, and Acinetobacter are presently resistant to all therapy except for the polymyxins, but polymyxin resistance is now being reported as well. New antibiotics are not being rapidly developed. One approach to help combat resistance is to determine synergy with the use of two or more presently available agents but doing so would not have practical value in caring for very sick patients, unless the synergy could be determined quickly.
Dr. George Pankey
gpankey@ochsner.org
Our Vision/Goals
The goal of the Infectious Disease Translational Research Laboratory is to develop rapid diagnostic methods for the determination of synergy or antagonism between two or more antimicrobial agents against multidrug-resistant bacteria and fungi. We developed an Etest® methodology that is very promising in this regard. Our research also involves evaluation of rapid diagnostic testing for direct-from-blood organisms and resistance genes.
Bacterial Synergy Testing
The laboratory is equipped with a Vitek® 2 instrument (which performs rapid automated identification and antimicrobial susceptibility testing of bacteria and fungi). An Eddy Jet Spiral plater with a Sphere Flash automated colony counter are used for time-kill studies to determine interactions between antimicrobial agents, along with the Etest and checkerboard methods. Recent synergy studies include: nisin plus polymyxin B against pan drug-resistant Acinetobacter baumannii; fosfomycin plus doxycycline against linezolid-resistant, vancomycin-resistant Enterococcus faecium; fosfomycin plus meropenem against metallo-ß-lactamase producing Pseudomonas aeruginosa; and lefamulin plus doxycycline against linezolid-resistant, vancomycin-resistant E. faecium. Our laboratory is also involved in the study of how quickly bacteria are killed by antibiotics and non-antibiotics with antibacterial activity.
Thomas VM, Brown RM, Ashcraft DS, Pankey GA. 2019. Synergistic effect between nisin and polymyxin B against pandrug-resistant and extensively drug-resistant Acinetobacter baumannii. Int J Antimicrob Agents. 53:663-668. https://doi.org/10.1016/j.ijantimicag.2019.03.009
Davis H, Brown R, Ashcraft D, Pankey G. 2020. In Vitro Synergy with Fosfomycin Plus Doxycycline Against Linezolid and Vancomycin-resistant Enterococcus faecium. J Global Antimicrob Resist.;22:78-83. https://doi.org/10.1016/j.jgar.2020.01.014
Candida Studies
Candida species are the fourth most common cause of health care-associated bloodstream infections in the US (2015), with Candida albicans predominating. However, at Ochsner since 2010, Candida glabrata has been the most common species. C. glabrata is more resistant to certain antifungals, creating a concern regarding appropriate therapy, because susceptibility studies are not routinely performed by the diagnostic microbiology laboratory. The Infectious Disease research laboratory is performing antifungal susceptibility testing on these bloodstream isolates. Specific aims of the study include determining if there is a trend in Candida species being isolated more frequently from the blood and if any antifungal resistance is developing. We will determine if there is any effect of the time of diagnosis and time of initial antifungal therapy on mortality of patients with Candida bloodstream infections. This study will provide useful information that may improve treatments for future patients. Synergy studies include fluconazole plus polymyxin B or doxycycline or tigecycline or caspofungin against C. glabrata; and fluconazole plus polymyxin B or trimethoprim-sulfamethoxazole against Candida auris.
Hooper RW, Ashcraft DS, Pankey GA. 2018. In vitro synergy with fluconazole plus doxycycline or tigecycline against clinical Candida glabrata isolates. Med Mycol. Mar:1-5. https://doi.org/10.1093/mmy/myy008
Adams EK, Ashcraft DS, Pankey GA. 2016. In vitro Synergistic Activity of Caspofungin Plus Polymyxin B Against Fluconazole-Resistant Candida glabrata. Am J Med Sciences. Mar:351(3):265-9. https://doi.org/10.1016/j.amjms.2015.12.014
Pankey G, Ashcraft D, Kahn H, Ismail A. 2014. Time-kill assay and Etest evaluation for synergy with polymyxin B and fluconazole against Candida glabrata. Antimicrob. Agents Chemother. Oct;58(10):5795-800. https://doi.org/10.1128/AAC.03035-14
Molecular Cardiovascular
The Infectious Disease research laboratory has a T2MR® instrument (T2Dx®, T2 Biosystems) that will rapidly detect Candida spp. and 90% of deadly ESKAPE pathogens (E. faecium, Staphylococcus aureus, Klebsiella pneumoniae, P. aeruginosa, and Escherichia coli) in the blood in 3-5 hrs. This will enable early and more appropriate treatment decisions, since up to 5 days are required for identification from blood culture. Our laboratory was one of 11 U.S. acute care hospitals who participated in the clinical trial for the T2Bacteria Panel, which was FDA-cleared in 2018. The T2Resistance Panel- Research use only (RUO), which can detect 13 resistance genes from both gram-positive and gram-negative pathogens directly from blood, is available for future research studies. Other products available with this T2MR technology include the T2Cauris Panel RUO, which detects the rapidly emerging multidrug-resistant Candida auris in patient skin and blood, in addition to hospital environmental samples in <5 hrs and the T2SARS-CoV-2 Panel, designed for direct detection of nucleic acid from SARS-CoV-2 in upper respiratory specimens in <2 hrs.
M. Hong Nguyen, MD; Cornelius J. Clancy, MD; A. William Pasculle, ScD; Peter G. Pappas, MD; George Alangaden, MD; George A. Pankey, MD; Bryan H. Schmitt, DO; Altaf Rasool, MD; Melvin P. Weinstein, MD; Raymond Widen, MD; Diana R. Hernandez, PhD; Donna M. Wolk, PhD; Thomas J. Walsh, MD; John R. Perfect, MD; Mollie N. Wilson, MS; and Eleftherios Mylonakis, MD. Performance of the T2Bacteria Panel for Diagnosing Bloodstream Infections. A Diagnostic Accuracy Study. 2019. Ann Intern Med. https://doi.org/10.7326/M18-2772
Voigt C, Silbert S, Widen RH, Marturano JE, Lowery TJ, Ashcraft D, Pankey G. 2020. The T2Bacteria Assay Is a Sensitive and Rapid Detector of Bacteremia That Can Be Initiated in the Emergency Department and Has Potential to Favorably Influence Subsequent Therapy. J Emerg Med May;58(5):785-796. https://doi.org/10.1016/j.jemermed.2019.11.028 Epub 2020 Jan PMID: 31982197
Other Studies
The Infectious Disease research laboratory continues to collect and bank bacterial and fungal isolates from blood cultures—as well as unusual or multidrug-resistant isolates from other specimen sources. Over 25,000 isolates are available for use in future studies. This research laboratory serves as a reference laboratory to the Ochsner diagnostic microbiology laboratory to expand testing capabilities on patients' cultures and as a teaching laboratory for undergraduate and medical students, residents, and fellows interested in Infectious Disease translational research.
The laboratory also is a regional site for the collection of bacterial and fungal isolates for the SENTRY Antimicrobial Surveillance Program for USA/Core Lab, representing the SW region of the U.S. (JMI Laboratories, North Liberty, IA).
SENTRY Program Participating Sites (1997–2016)
Open Forum Infectious Diseases, Volume 6, Issue Supplement_1, March 2019, Pages S95–S102, https://doi.org/10.1093/ofid/ofz058
For more information, please contact Dr. George Pankey (gpankey@ochsner.org).
Research Programs ▼
Core Laboratories ▼
Biomarker Core Laboratory
Flow Cytometry and In Vivo Imaging Core Laboratory
Microscopy Core Laboratory


text:
Home / Money & Career /
12 Vehicles Racing to Be 2014′s Truck/SUV of the Year
Justin Lloyd-Miller Google+ Twitter
It's that time of year again, when an official consortium of auto journalists and writers get together to determine what the best car and truck for 2014 will be. A total of 48 members make up the panel that will ultimately bestow the honor on a lucky vehicle.
"The awards — first given in 1994 — are unique in North America because — instead of being given by a single publication, web site, radio or television program — they are given by automotive journalists representing all of those outlets in two countries," the North American Car and Truck of the Year website says. The two countries are the United States and Canada.
"The awards recognize the most outstanding vehicles of the year based on factors including innovation, design, safety, handling, driver satisfaction and value for the dollar," the website reads. "To be eligible, vehicles must be 'all new' or 'substantially changed.'"
We've already covered the 12 cars up for the Car of the Year Award, but that's only half the story. Check out the 12 trucks and SUVs competing for the 2014 Truck of the Year Award.
1. Acura MDX
Acura (NYSE:HMC) seems to be still fine-tuning its sedan lineup, but when it comes to utilities, Honda's luxury arm, has it down. The new MDX, which starts at around $42,000, is a tech-laden, seven-passenger coach with a 290 horsepower beating heart. Each MDX comes standard with Acura's Jewel Eye headlights, a multi-view rear camera, and 28 miles per gallon on the highway.
2. BMW X5
Having been around for nearly 14 years now, the BMW X5 has established itself as a benchmark of luxury SUVs in North America. It comes standard with an inline-six and 300 horsepower at $52,800, but buyers can also spring for a diesel version that gets 31 miles per gallon highway or the more powerful, 445 horsepower X5 xDrive50i.
3. Buick Encore
Buick's (NYSE:GM) compact crossover offers about the most luxury for the dollar as you'll find, with a base price of $24,160. If you're looking to gussy it up a bit, though, Buick offers a lengthy menu of premium options that will make the Encore — which manages 33 miles per gallon highway and 25 in the city — a truly competitive luxury crossover.
4. Chevrolet Silverado
The new Silverado is Chevrolet's latest foray in the highly competitive light-duty pickup market, and it's proven to be a formidable foe. Though it doesn't sell in Ford F-150 numbers, the 5.3-liter EcoTec V8-fitted Silverado is selling faster than General Motors can produce them. Chevrolet managed to put a modern spin on the vehicle without sacrificing the rugged appearance that a truck should have.
5. GMC Sierra
Pitted against the Silverado for the Truck of the Year Award is its corporate sibling, the GMC Sierra. The Sierra starts at a base price of $26,075 and comes standard with a 4.3-liter EcoTec3 V6, which is good for 285 horsepower and 305 pound-feet of torque. The GMC is aimed more at professionals, though underneath the trimmings and trappings, the Sierra and Silverado are quite similar.
6. Hyundai Santa Fe LWB
The Hyundai Santa Fe has become a worthy contender for the award through the value that it brings consumers, offering leading horsepower figures (290) for its class, room for seven passengers, and a healthy list of standard options, all for less than $30,000 at base.
7. Jeep Cherokee
The new Jeep Cherokee certainly looks striking, and it represents Jeep's renewed commitment to the compact SUV segment. With a base MSRP of $22,995, it's also one of the most affordable new Jeep models, offering a more civilized take on an offroad SUV that won't deplete the bank account.
8. Kia Sorento
Hand in hand with the Hyundai comes the Kia Sorento, which has matured into a fine vehicle that — like the Hyundai — goes above and beyond what would be insinuated from its MSRP of $24,100. The Kia comes standard with heated mirrors and solar control glass, which reportedly keeps the interior cooler on hot, sunny days.
9. Land Rover Range Rover Sport
Long the standard of the luxury SUV, the Range Rover Sport is back in a new suit and all the comfort that one would expect from Land Rover. A 3-liter V6 is good for 340 horsepower, though a 5-liter supercharged V8 is also optional (510 horsepower). Alongside its raw power, the Range Rover also boasts a sumptuous interior, decked out with just about every luxury one can ask for in an automobile.
10. Nissan Rogue
Nissan's (NSANY.PK) previous Rogue proved to be an enduring model that remained popular right through this year (Nissan will continue to produce it alongside the 2014 model), but it's hard to argue that it wasn't due for a refresh. The latest update sees the bolder, more aggressive look that's becoming an industry standard, though it still starts at a reasonable $22,490 base price.
11. Subaru Forester
As far as the apex of capability, efficiency, and value go, it's hard to do much better than the Subaru Forester, which only seemed to get better with the 2014 model year. It now boasts improved mileage (32 highway, 24 in the city), retains the standard all-wheel drive that Subaru is renowned for, and starts at a reasonable $21,995.
12. Toyota Tundra
Toyota's (NYSE:TM) Tundra was badly needing a redesign to remain competitive, and the 2014 model appears to have been given the model the shot it needed. Powertrain options are still limited but the Tundra is as muscular as ever and looks the part. It starts at $25,920, and decked with a 5.7-liter V8, it'll tow just about anything you throw at it.
Don't Miss: 9 Upscale Small Cars to Splurge On.


text:
The venue was filled with over 35 corgi dogs racing for a treat at Martensville Curling Club last Saturday.
"We just love Corgis, we think they're a really unique breed, and they're actually becoming more and more trendy these days, so we just wanted people to come out and enjoy a laugh," said Kate Kading and Sask Corgi Racers Group Media Liaison, in which this breed is getting a lot of attention on the event during the Queen's Platinum Jubilee.
Kading said that this event was held to raise money and the amount earned will be donated for the Saskatoon SPCA.
However, she didn't expect that there will be over 600 people attending the said event, she added. "I was expecting around 50-100 people. It was actually overwhelming how many people showed up today to support this," Kading said.
Totally, they raised around $3,000 after the event. There were different types of corgis that compete with each other just like Pembroke welsh corgis, cowboy corgis, cardigan corgis and a lab-corgi mix to name some.
The winner of the event, Stevie, was asked how it felt to clinch the title. Stevie just said,"*Lick*". Stevie's owners did expect that they will win the competition, thus were very proud of the result and the performance of Stevie.
"Stevie dominated the competition at her first Sask dog race, and it was very exhilarating and clearly the strict training regimen we put her through paid off this year, it was a tough camp, but clearly she made it worth it," Kade said, the owner of Stevie.
Hamilton, another competitor was crowned 2nd place and left his mark.."*Lick*," Hamilton responded.
"I've been wanting Corgi get-togethers in Saskatoon or the area for quite a long time, I've had these guys for 4-5 years, and you always see the meet ups in other places, so the videos are great, and the people are always really good," Megan said, the owner of Hamilton.
They can possibly invite one of Queen Elizabeth's corgis to be Grand Marshall in another
corgi race that will happen next year.
Source: Global News