lsst::ap::cluster Namespace Reference

Namespaces
	config
	detail

Classes
class	SourceAndExposure
struct	ClusteringControl
	Parameters for the clustering algorithm and its internals. More...
class	SourceClusterColumnViewT
struct	KeyTuple
	A three-element tuple of mean, uncertainty, and count keys. More...
class	SourceClusterRecord
	Record class that contains measurement averages on clusters of single exposure sources. More...
class	SourceClusterTable
	Table class that contains measurement means on clusters of single exposure sources. More...
class	SourceClusterIdFactory
struct	SourceProcessingControl
	Parameters for source processing. More...

Typedefs
typedef SourceClusterColumnViewT < SourceClusterRecord >	SourceClusterColumnView
typedef lsst::afw::table::SortedCatalogT < SourceClusterRecord >	SourceClusterCatalog
typedef lsst::afw::table::SortedCatalogT < SourceClusterRecord const >	ConstSourceClusterCatalog

Functions
boost::shared_ptr< std::vector < SourceAndExposure > > const	computeBasicAttributes (SourceClusterRecord &cluster, SourceCatalog const &sources, ExposureInfoMap const &exposures, std::string const &exposurePrefix)
void	computeFluxMean (SourceClusterRecord &cluster, std::vector< SourceAndExposure > const &sources, std::string const &fluxDef, std::vector< Key< Flag > > const &skipFlags, double fluxScale)
void	computeShapeMean (SourceClusterRecord &cluster, std::vector< SourceAndExposure > const &sources, std::string const &shapeDef, std::vector< lsst::afw::table::Key< lsst::afw::table::Flag > > const &skipFlags)
std::pair< boost::shared_ptr < SourceTable >, SchemaMapper > const	makeOutputSourceTable (SourceTable const &prototype, SourceProcessingControl const &control)
boost::shared_ptr < SourceClusterTable > const	makeSourceClusterTable (SourceTable const &prototype, boost::shared_ptr< IdFactory > const &idFactory, SourceProcessingControl const &control)
void	processSources (SourceCatalog const &expSources, ExposureInfo const &expInfo, lsst::ap::utils::PT1SkyTile const *skyTile, SourceProcessingControl const &control, SchemaMapper const &mapper, SourceCatalog &sources, SourceCatalog &badSources, SourceCatalog &invalidSources)
std::vector< SourceCatalog > const	cluster (SourceCatalog const &sources, ClusteringControl const &control)
void	setClusterFields (SourceCatalog &sources, SourceClusterRecord const &record, SourceProcessingControl const &control)
KeyTuple< lsst::afw::table::Shape >	addShapeFields (lsst::afw::table::Schema &schema, std::string const &filter, std::string const &name, std::string const &doc)
	Convenience function to setup fields for shapes. More...
KeyTuple< lsst::afw::table::Flux >	addFluxFields (lsst::afw::table::Schema &schema, std::string const &filter, std::string const &name, std::string const &doc, std::string const &unit)
	Convenience function to setup fields for fluxes. More...
boost::shared_ptr< std::vector < SourceAndExposure > > const	computeBasicAttributes (SourceClusterRecord &cluster, lsst::afw::table::SourceCatalog const &sources, lsst::ap::match::ExposureInfoMap const &exposures, std::string const &exposurePrefix)
void	computeFluxMean (SourceClusterRecord &cluster, std::vector< SourceAndExposure > const &sources, std::string const &fluxDef, std::vector< lsst::afw::table::Key< lsst::afw::table::Flag > > const &skipFlags, double fluxScale)
boost::shared_ptr < SourceClusterTable > const	makeSourceClusterTable (lsst::afw::table::SourceTable const &prototype, boost::shared_ptr< lsst::afw::table::IdFactory > const &idFactory, SourceProcessingControl const &control)
void	processSources (lsst::afw::table::SourceCatalog const &expSources, lsst::ap::match::ExposureInfo const &expInfo, lsst::ap::utils::PT1SkyTile const *skyTile, SourceProcessingControl const &control, lsst::afw::table::SchemaMapper const &mapper, lsst::afw::table::SourceCatalog &sources, lsst::afw::table::SourceCatalog &badSources, lsst::afw::table::SourceCatalog &invalidSources)
std::vector < lsst::afw::table::SourceCatalog > const	cluster (lsst::afw::table::SourceCatalog const &sources, ClusteringControl const &control)
void	setClusterFields (lsst::afw::table::SourceCatalog &sources, SourceClusterRecord const &record, SourceProcessingControl const &control)

Typedef Documentation

typedef lsst::afw::table::SortedCatalogT<SourceClusterRecord const> lsst::ap::cluster::ConstSourceClusterCatalog

Definition at line 858 of file SourceCluster.h.

typedef lsst::afw::table::SortedCatalogT<SourceClusterRecord> lsst::ap::cluster::SourceClusterCatalog

Definition at line 857 of file SourceCluster.h.

typedef SourceClusterColumnViewT<SourceClusterRecord> lsst::ap::cluster::SourceClusterColumnView

Definition at line 826 of file SourceCluster.h.

Function Documentation

KeyTuple< lsst::afw::table::Flux > lsst::ap::cluster::addFluxFields	(	lsst::afw::table::Schema &	schema,
		std::string const &	filter,
		std::string const &	name,
		std::string const &	doc,
		std::string const &	unit
	)

Convenience function to setup fields for fluxes.

Definition at line 435 of file SourceCluster.cc.

{
    using lsst::afw::table::Flux;
    KeyTuple<Flux> kt;
    kt.mean = schema.addField<Flux::MeasTag>(
        filter + "." + name, doc, unit);
    kt.err = schema.addField<Flux::ErrTag>(
        filter + "." + name + ".err",
        "uncertainty for " + filter + "." + name,
        unit);
    kt.count = schema.addField<int>(
        filter + "." + name + ".count",
        "Number of samples used to compute the " + filter +
        "." + name + " mean");
    return kt;
}

KeyTuple< lsst::afw::table::Shape > lsst::ap::cluster::addShapeFields	(	lsst::afw::table::Schema &	schema,
		std::string const &	filter,
		std::string const &	name,
		std::string const &	doc
	)

Convenience function to setup fields for shapes.

Definition at line 414 of file SourceCluster.cc.

{
    using lsst::afw::table::Shape;
    KeyTuple<Shape> kt;
    kt.mean = schema.addField<Shape::MeasTag>(
        filter + "." + name, doc, "rad^2");
    kt.err = schema.addField<Shape::ErrTag>(
        filter + "." + name + ".err",
        "covariance matrix for " + filter + "." + name,
        "rad^4");
    kt.count = schema.addField<int>(
        filter + "." + name + ".count",
        "Number of samples used to compute the " + filter + 
        "." + name + " mean");
    return kt;
}

std::vector<lsst::afw::table::SourceCatalog> const lsst::ap::cluster::cluster	(	lsst::afw::table::SourceCatalog const &	sources,
		ClusteringControl const &	control
	)

Spatially cluster sources using the OPTICS algorithm. For details, see the following paper:

"OPTICS: Ordering Points To Identify the Clustering Structure". Mihael Ankerst, Markus M. Breunig, Hans-Peter Kriegel, Jorg Sander (1999). ACM SIGMOD international conference on Management of data. ACM Press. pp. 49-60.

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.129.6542 http://en.wikipedia.org/wiki/OPTICS_algorithm

Parameters

[in]	sources	Sources to cluster.
[in]	control	Clustering parameters.

std::vector<SourceCatalog> const lsst::ap::cluster::cluster	(	SourceCatalog const &	sources,
		ClusteringControl const &	control
	)

Definition at line 578 of file clustering.cc.

{
    typedef SourceCatalog::const_iterator Iter;

    if (sources.size() > MAX_SOURCES) {
        throw LSST_EXCEPT(InvalidParameterError, "too many sources to cluster");
    }
    control.validate();
    boost::scoped_array<OpticsPoint> entries(new OpticsPoint[sources.size()]);
    std::vector<SourceCatalog> clusters;
    int i = 0;
    for (Iter s = sources.begin(), e = sources.end(); s != e; ++s, ++i) {
        entries[i].coords = s->getCoord().getVector().asEigen();
        entries[i].data = s;
    }
    if (i > 0) {
        // Convert epsilon and leafExtentThreshold to radians, and account
        // for the fact that our metric is the squared euclidian distance,
        // not angular separation.
        double eps = std::sin(0.5 * control.getEpsilon().asRadians());
        eps = 4.0 * eps * eps;
        double let = control.getLeafExtentThreshold().asRadians();
        if (let > 0.0) {
            let = std::sin(0.5 * let);
            let = 4.0 * let * let;
        }
        Optics optics(entries.get(), i, control.minNeighbors, eps, let, control.pointsPerLeaf);
        optics.run(sources.getTable(), clusters, detail::SquaredEuclidianDistanceOverSphere());
    }
    return clusters;
}

boost::shared_ptr<std::vector<SourceAndExposure> > const lsst::ap::cluster::computeBasicAttributes	(	SourceClusterRecord &	cluster,
		lsst::afw::table::SourceCatalog const &	sources,
		lsst::ap::match::ExposureInfoMap const &	exposures,
		std::string const &	exposurePrefix
	)

Compute basic cluster attributes.

The Coord slot of cluster is set to the unweighted mean of the coordinates of the associated sources. This fiducial cluster position is used as the center of a N,E tangent plane projection that serves as a common coordinate system for averaging other measurements. Exposure information for each input source is looked up, and used to create a vector of tuples that bundle a source with it's originating exposure and an affine transform from the pixel space coordinate system centered on the source's centroid to the common coordinate system.

This vector of tuples is sorted such that sources originating from exposures in the same filter are all consecutive.

Finally, several other cluster attributes are computed/set:

• "obs.count" : number of sources in the cluster
• "obs.time.min", "obs.time.max", "obs.time.mean" : observation time range
• "coord.weighted", "coord.weighted.err", "coord.weighted.count" : inverse variance weighted mean coordinate, error, and sample count (computed only if input sources have a valid CoordErr slot mapping).
• "\<filter\>.obs.count" : filter specific source count
• "\<filter\>.time.min", "\<filter\>.time.max" : filter specific observation time range

Parameters

[out]	cluster	SourceClusterRecord to store results in.
[in]	sources	Sources in cluster.
[in]	exposures	Map from exposure ID to exposure information.
[in]	exposurePrefix	Exposure column name prefix.

Returns

A list of (source, exposure information, transform) tuples, grouped by filter.

boost::shared_ptr<std::vector<SourceAndExposure> > const lsst::ap::cluster::computeBasicAttributes	(	SourceClusterRecord &	cluster,
		SourceCatalog const &	sources,
		ExposureInfoMap const &	exposures,
		std::string const &	exposurePrefix
	)

Definition at line 334 of file attributes.cc.

{
    typedef SourceCatalog::const_iterator SourceIter;
    typedef std::vector<SourceAndExposure>::const_iterator SeIter;

    if (sources.empty()) {
        throw LSST_EXCEPT(InvalidParameterError, "No sources in cluster");
    }
    // unweighted mean coordinates
    meanCoord(cluster, sources);

    boost::shared_ptr<std::vector<SourceAndExposure> > se(
        new std::vector<SourceAndExposure>);
    NeTanProj const proj(cluster.getCoord());
    Schema const schema = sources.getSchema();
    Key<int64_t> const expIdKey = schema.find<int64_t>(exposurePrefix + ".id").key;
    Key<double> const expTimeMidKey = schema.find<double>(exposurePrefix + ".time.mid").key;
    // fill in vector of source,exposure pairs
    se->reserve(sources.size());
    for (SourceIter i = sources.begin(), e = sources.end(); i != e; ++i) {
        boost::shared_ptr<ExposureInfo> exp =
            const_cast<ExposureInfoMap &>(exposures).get(i->get(expIdKey));
        if (!exp) {
            throw LSST_EXCEPT(NotFoundError, "No ExposureInfo for source");
        }
        se->push_back(SourceAndExposure(
            i, exp, proj.pixelToNeTransform(i->getCentroid().asEigen(), *exp->getWcs())));
    }
    // source counts and observation time range
    StatTuple t;
    t.compute(se->begin(), se->end(), expTimeMidKey);
    cluster.setNumSources(static_cast<int>(sources.size()));
    cluster.setTimeMin(t.min);
    cluster.setTimeMean(t.mean);
    cluster.setTimeMax(t.max);
    // inverse variance weighted mean coordinates
    if (sources.getTable()->getCentroidKey().isValid() &&
        sources.getTable()->getCentroidErrKey().isValid()) {
        weightedMeanCoord(cluster, *se, proj);
    }
    // sort source,exposure vector by filter
    std::sort(se->begin(), se->end(), CmpSourceAndExposure());
    // per-filter source counts and observation time range
    SeIter i = se->begin();
    SeIter const e = se->end();
    while (i != e) {
        int const filterId = i->getExposureInfo()->getFilter().getId();
        SeIter j = i;
        for (++j; j != e && filterId == j->getExposureInfo()->getFilter().getId(); ++j) { }
        std::string const filterName = i->getExposureInfo()->getFilter().getName();
        cluster.setNumSources(filterName, static_cast<int>(j - i));
        t.compute(i, j, expTimeMidKey);
        cluster.setTimeMin(filterName, t.min);
        cluster.setTimeMax(filterName, t.max);
        i = j;
    }
    return se;
}

void lsst::ap::cluster::computeFluxMean	(	SourceClusterRecord &	cluster,
		std::vector< SourceAndExposure > const &	sources,
		std::string const &	fluxDef,
		std::vector< lsst::afw::table::Key< lsst::afw::table::Flag > > const &	skipFlags,
		double	fluxScale
	)

Compute per-filter calibrated flux means for a cluster of sources.

The input source vector is assumed to be grouped by filter. Flux measurments and errors are obtained from fields named '<fluxDef>' and '<fluxDef>.err'. Results are stored in fields named '<filter>.<fluxDef>' and '<filter>.<fluxDef>.err', with a sample count stored in '<filter>.<fluxDef>.count' field. For an input source to be included in the mean, the following must hold:

• the source must have finite flux
• the source must have positive flux error
• none of the flag fields in skipFlags must be set.

If, after applying these filtering critera, 2 or more flux samples are available, then the inverse variance weighted mean of the samples (after calibration) is computed and stored. If only one sample is available, the input flux and error are calibrated and stored directly. If no samples are available, the output flux, error, and sample count fields are set to NaN, NaN, and 0.

Parameters

[out]	cluster	Cluster to store results in.
[in]	sources	List of sources and their exposures.
[in]	fluxDef	Name of flux field to operate on.
[in]	skipFlags	Flags marking sources to skip.
[in]	fluxScale	Flux scaling factor.

void lsst::ap::cluster::computeFluxMean	(	SourceClusterRecord &	cluster,
		std::vector< SourceAndExposure > const &	sources,
		std::string const &	fluxDef,
		std::vector< Key< Flag > > const &	skipFlags,
		double	fluxScale
	)

Definition at line 439 of file attributes.cc.

{
    typedef std::vector<SourceAndExposure>::const_iterator Iter;
    typedef std::vector<Key<Flag > >::const_iterator FlagIter;

    if (sources.empty()) {
        throw LSST_EXCEPT(InvalidParameterError, "No sources in cluster");
    }
    Schema const sourceSchema = sources[0].getSource()->getSchema();
    Schema const clusterSchema = cluster.getSchema();
    Flux::MeasKey const sourceFluxKey = sourceSchema[fluxDef];
    Flux::ErrKey const sourceFluxErrKey = sourceSchema[fluxDef + ".err"];

    Iter i = sources.begin();
    Iter const e = sources.end();
    while (i != e) {
        // Determine range [i, j) containing sources from the same filter
        int const filterId = i->getExposureInfo()->getFilter().getId();
        Iter j = i;
        for (++j; j != e && filterId == j->getExposureInfo()->getFilter().getId(); ++j) { }
        // iterate over [i, j), storing all valid flux/error pairs
        // after calibrating them.
        std::vector<FluxAndVariance> samples;
        samples.reserve(j - i);
        std::string const filter = i->getExposureInfo()->getFilter().getName();
        for (; i < j; ++i) {
            SourceRecord const * r = i->getSource().get();
            double flux = r->get(sourceFluxKey);
            double fluxErr = r->get(sourceFluxErrKey);
            if (!lsst::utils::isfinite(flux) || lsst::utils::isnan(fluxErr) || fluxErr <= 0.0) {
                continue;
            }
            bool skip = false;
            for (FlagIter f = skipFlags.begin(), fe = skipFlags.end(); f != fe; ++f) {
                if (r->get(*f)) {
                    skip = true;
                    break;
                }
            }
            if (skip) {
                continue;
            }
            samples.push_back(i->getExposureInfo()->calibrateFlux(flux, fluxErr, fluxScale));
        }
        FluxAndErr mean = weightedMeanFlux(samples);
        Flux::MeasKey const measKey = clusterSchema[filter + "." + fluxDef];
        cluster.set(measKey, mean.first);
        Flux::ErrKey const errKey = clusterSchema[filter + "." + fluxDef + ".err"];
        cluster.set(errKey, mean.second);
        Key<int> const countKey = clusterSchema[filter + "." + fluxDef + ".count"];
        cluster.set(countKey, static_cast<int>(samples.size()));
    }
}

void lsst::ap::cluster::computeShapeMean	(	SourceClusterRecord &	cluster,
		std::vector< SourceAndExposure > const &	sources,
		std::string const &	shapeDef,
		std::vector< lsst::afw::table::Key< lsst::afw::table::Flag > > const &	skipFlags
	)

Compute per-filter shape means for a cluster.

The input source vector is assumed to be grouped by filter. Shape measurments and errors are obtained from fields named '<shapeDef>' and '<shapeDef>.err'. Results are stored in fields named '<filter>.<shapeDef>' and '<filter>.<shapeDef>.err', with a sample count stored in '<filter>.<shapeDef>.count' field. For an input source to be included in the mean, the following must hold:

• the source must have finite moments
• the source must have positive entries in the covariance matrix diagonals
• off diagonal elements of the covariance matrix must not be NaN
• the covariance matrix must be symmetric.
• none of the flag fields in skipFlags must be set.

If, after applying these filtering critera, 2 or more shape samples are available, then the inverse variance weighted mean of the samples (after transformation to a coordinate system S centered on the cluster position with the standard N,E basis) is computed and stored. If only one sample is available, the input shape and error are transformed to S and stored directly. If no samples are available, the output moments and covariance matrix are set to NaNs, and the sample count is zeroed.

Parameters

[out]	cluster	Cluster to store results in.
[in]	sources	List of sources and their exposures.
[in]	shapeDef	Name of shape field to operate on.
[in]	skipFlags	Flags marking sources to skip.

Definition at line 535 of file attributes.cc.

{
    typedef std::vector<SourceAndExposure>::const_iterator Iter;
    typedef std::vector<Key<Flag > >::const_iterator FlagIter;

    if (sources.empty()) {
        throw LSST_EXCEPT(InvalidParameterError, "No sources in cluster");
    }
    Schema const sourceSchema = sources[0].getSource()->getSchema();
    Schema const clusterSchema = cluster.getSchema();
    Shape::MeasKey const sourceShapeKey = sourceSchema[shapeDef];
    Shape::ErrKey const sourceShapeErrKey = sourceSchema[shapeDef + ".err"];

    Iter i = sources.begin();
    Iter const e = sources.end();
    while (i != e) {
        // Determine range [i, j) containing sources from the same filter
        int const filterId = i->getExposureInfo()->getFilter().getId();
        Iter j = i;
        for (++j; j != e && filterId == j->getExposureInfo()->getFilter().getId(); ++j) { }
        // iterate over [i, j) storing all valid shape/error pairs in samples
        std::vector<ShapeAndErr> samples;
        samples.reserve(j - i);
        std::string const filter = i->getExposureInfo()->getFilter().getName();
        for (; i < j; ++i) {
            SourceRecord const * r = i->getSource().get();
            bool skip = false;
            for (FlagIter f = skipFlags.begin(), fe = skipFlags.end(); f != fe; ++f) {
                if (r->get(*f)) {
                    skip = true;
                    break;
                }
            }
            if (skip) {
                continue;
            }
            Quadrupole q = r->get(sourceShapeKey);
            if (!lsst::utils::isfinite(q.getIxx()) ||
                !lsst::utils::isfinite(q.getIyy()) ||
                !lsst::utils::isfinite(q.getIxy())) {
                continue; // shape contains NaNs
            }
            Eigen::Matrix3d cov = r->get(sourceShapeErrKey).cast<double>();
            if (lsst::utils::isnan(cov(0,0)) ||
                lsst::utils::isnan(cov(1,1)) ||
                lsst::utils::isnan(cov(2,2))) {
                continue; // covariance matrix contains NaNs
            } else if (cov(0,0) <= 0.0 ||
                       cov(1,1) <= 0.0 ||
                       cov(2,2) <= 0.0) {
                continue; // negative variance
            }
            // FIXME: for now, no measurement algorithms actually
            //        compute full covariance matrixes, and the
            //        off diagonal matrix elements are always NaN.
            //        Longer term, we will need some algorithm metadata
            //        to decide whether an off-diagonal NaN means a
            //        sample should be ignored because a computation failed,
            //        or whether it should be zeroed because the algorithm
            //        never computes it.
            if (cov(0,1) != cov(1,0) || cov(0,2) != cov(2,0) || cov(1,2) != cov(2,1)) {
                if (lsst::utils::isnan(cov(0,1)) && lsst::utils::isnan(cov(1,0)) &&
                    lsst::utils::isnan(cov(0,2)) && lsst::utils::isnan(cov(2,0)) &&
                    lsst::utils::isnan(cov(1,2)) && lsst::utils::isnan(cov(2,1))) {
                    cov(0,1) = 0.0; cov(1,0) = 0.0;
                    cov(0,2) = 0.0; cov(2,0) = 0.0;
                    cov(1,2) = 0.0; cov(2,1) = 0.0;
                } else {
                    continue; // covariance matrix not symmetric
                }
            }
            // transform moments and covariance matrix to N,E basis
            LinearTransform const * xform = &(i->getTransform().getLinear());
            Eigen::Matrix3d j = q.transform(*xform).d();
            q.transform(*xform).inPlace();
            cov = j * cov * j.transpose();
            samples.push_back(ShapeAndErr(q, cov));
        }
        ShapeAndErr mean;
        weightedMeanShape(mean, samples);
        Shape::MeasKey const measKey = clusterSchema[filter + "." + shapeDef];
        cluster.set(measKey, mean.first);
        Shape::ErrKey const errKey = clusterSchema[filter + "." + shapeDef + ".err"];
        cluster.set(errKey, mean.second);
        Key<int> const countKey = clusterSchema[filter + "." + shapeDef + ".count"];
        cluster.set(countKey, static_cast<int>(samples.size()));
    }
}

std::pair< boost::shared_ptr< lsst::afw::table::SourceTable >, lsst::afw::table::SchemaMapper > const lsst::ap::cluster::makeOutputSourceTable	(	lsst::afw::table::SourceTable const &	prototype,
		SourceProcessingControl const &	control
	)

Create an output source table based on a prototypical input source table.
The schema and slot mappings from the input table are used as is, and exposure
related fields, required for database ingestion / clustering, are added.

In particular, the following additions are performed:

    - "coord.err"             (CovPointD) : Source sky-coordinate error
    - "<exposure>.id"         (L)         : ID of exposure source was measured on.

i * - "<cluster>.id" (L) : ID of cluster containing source.

• "<cluster>.coord" (Coord) : ICRS coordinates of cluster containing source.

If the control.multiBand config parameter is false:

• "<exposure>.filter.id" (I) : ID of filter for exposure.

If the control.coadd config parameter is false:

• "<exposure>.time" (F) : Exposure time (s).
• "<exposure>.time.mid" (D) : Middle of exposure time (MJD).

"<exposure>" and "<cluster>" are field name prefixes obtained from SourceProcessingControl. If either prefix is empty, the corresponding fields will not be added. Note however that the exposure related fields are necessary for attribute computation and database ingest into the LSST schema. The cluster related fields are required only for database ingest.

Parameters

[in]	prototype	Prototypical schema/slot mappings.
[in]	control	Source processing parameters.

Returns

An output SourceTable and a SchemaMapper that can be used to copy common field values between input and output records.

Definition at line 135 of file clustering.cc.

{
    if (!prototype.getCoordKey().isValid() ||
        !prototype.getCentroidKey().isValid()) {
        throw LSST_EXCEPT(InvalidParameterError, "Prototypical "
            "SourceCatalog must have valid Coord and Centroid fields");
    }
    SchemaMapper mapper(prototype.getSchema());
    // add all fields in the input schema to the output schema.
    mapper.addMappingsWhere(True());
    // add coordinate error
    mapper.addOutputField(Field<Covariance<lsst::afw::table::Point<float> > >(
        "coord.err",
        "covariance matrix for source sky-coordinates",
        "rad^2"));
    if (!control.exposurePrefix.empty()) {
        // add exposure related fields
        mapper.addOutputField(Field<int64_t>(
            control.exposurePrefix + ".id",
            "ID of exposure source was measured on"));
        if (!control.multiBand) {
            mapper.addOutputField(Field<int>(
                control.exposurePrefix + ".filter.id",
                "ID of filter for exposure"));
        }
        if (!control.coadd) {
            mapper.addOutputField(Field<float>(
                control.exposurePrefix + ".time",
                "exposure time",
                "s"));
            mapper.addOutputField(Field<double>(
                control.exposurePrefix + ".time.mid",
                "middle of exposure time",
                "mjd"));
        }
    }
    if (!control.clusterPrefix.empty()) {
        // add cluster related fields
        mapper.addOutputField(Field<int64_t>(
            control.clusterPrefix + ".id",
            "ID of cluster containing source; 0 if source is not assigned to any cluster."));
        mapper.addOutputField(Field<lsst::afw::table::Coord>(
            control.clusterPrefix + ".coord",
            "ICRS sky-coordinates of cluster containing source",
            "rad"));
    }
    boost::shared_ptr<SourceTable> table =
        SourceTable::make(mapper.getOutputSchema(), boost::shared_ptr<IdFactory>());
    // copy slot mappings from prototype
    table->definePsfFlux(prototype.getPsfFluxDefinition());
    table->defineModelFlux(prototype.getModelFluxDefinition());
    table->defineApFlux(prototype.getApFluxDefinition());
    table->defineInstFlux(prototype.getInstFluxDefinition());
    table->defineCentroid(prototype.getCentroidDefinition());
    table->defineShape(prototype.getShapeDefinition());

    return std::make_pair(table, mapper);
}

boost::shared_ptr<SourceClusterTable> const lsst::ap::cluster::makeSourceClusterTable	(	lsst::afw::table::SourceTable const &	prototype,
		boost::shared_ptr< lsst::afw::table::IdFactory > const &	idFactory,
		SourceProcessingControl const &	control
	)

Create a source cluster table based on a prototypical input source table. The names of flux/shape fields from SourceProcessingControl which correspond to fields in the prototype are used to create flux/shape fields for every filter that has been defined via the lsst::afw::image::Filter API. Slots from the prototype are used to setup corresponding filter specific slots in the source cluster table.

Note that filters are typically defined by the data butler's mapper class (see e.g. the obs_lsstSim package). The names of flux/shape fields in the cluster schema are prefixed with filter names - for example, if the PSF flux slot in the prototype is mapped to a field named "flux.naive" and filters "u" and "g" are defined, then the cluster schema will contain fields "u.flux.naive" and "g.flux.naive", with filter specific PSF flux slots mapped accordingly.

The following additional fields will be created (names not configurable at the moment):

- "id"                  (L)         : cluster ID (from minimal schema)
- "coord"               (Coord)     : mean sky-coordinates (from minimal schema)
- "coord.err"           (CovPointD) : uncertainty of coord
- "obs.count"           (I)         : number of sources in cluster
- "flag.noise"          (Flag)      : was cluster created from a single noise source?
- "<filter>.obs.count"  (I)         : numer of sources in a specific filter

If the input table has valid centroid and centroid error slot keys, then the following fields are added:

- "coord.weightedmean"       (Coord)     : inverse variance weighted mean sky-coordinates
- "coord.weightedmean.err"   (CovPointD) : covariance matrix for coord.weighted
- "coord.weightedmean.count" (I)         : number of samples included in "coord.weighted"

Finally, if the "<exposure>.time.mid" field exists in the input table, then:

- "obs.time.min"          (D)       : earliest observation time of sources in cluster
- "obs.time.mean"         (D)       : mean observation time of sources in cluster
- "obs.time.max"          (D)       : latest observation time of sources in cluster
- "<filter>.obs.time.min" (D)       : earliest observation time in a specific filter
- "<filter>.obs.time.max" (D)       : latest observation time in a specific filter

are added (where the "<exposure>" prefix is obtained from SourceProcessingControl).

Parameters

[in]	prototype	Prototypical schema/slot mappings.
[in]	idFactory	ID generator.
[in]	control	Source processing parameters.

Returns

A SourceClusterTable corresponding to clusters of sources having the given prototypical schema and slots.

boost::shared_ptr<SourceClusterTable> const lsst::ap::cluster::makeSourceClusterTable	(	SourceTable const &	prototype,
		boost::shared_ptr< IdFactory > const &	idFactory,
		SourceProcessingControl const &	control
	)

Definition at line 197 of file clustering.cc.

{
    typedef std::vector<std::string>::const_iterator Iter;

    Schema schema = SourceClusterTable::makeMinimalSchema();
    schema.setVersion(0);
    // Add basic fields
    schema.addField<Covariance<lsst::afw::table::Point<float> > >(
        "coord.err",
        "sample covariance matrix for coord field (unweighted mean sky-coordinates)",
        "rad^2");
    schema.addField<int>(
        "obs.count",
        "number of sources in cluster");
    schema.addField<Flag>(
        "flag.noise",
        "set if cluster was created from a single noise source");
    if (prototype.getCentroidKey().isValid() &&
        prototype.getCentroidErrKey().isValid()) {
        schema.addField<lsst::afw::table::Coord>(
            "coord.weightedmean",
            "inverse variance weighted mean sky-coordinates (ICRS)",
            "rad");
        schema.addField<Covariance<lsst::afw::table::Point<float> > >(
            "coord.weightedmean.err",
            "covariance matrix for coord.weightedmean field",
            "rad^2");
        schema.addField<int>(
            "coord.weightedmean.count",
            "number of samples used to compute coord.weightedmean");
    }
    if (!control.exposurePrefix.empty()) {
        schema.addField<double>(
            "obs.time.min",
            "earliest observation time of sources in cluster (TAI)",
            "mjd");
        schema.addField<double>(
            "obs.time.mean",
            "mean observation time of sources in cluster (TAI)",
            "mjd");
        schema.addField<double>(
            "obs.time.max",
            "latest observation time of sources in cluster (TAI)",
            "mjd");
    }
    std::vector<std::string> filters = Filter::getNames();
    // Sort filter names by ID so that fields are added in a consistent order
    std::sort(filters.begin(), filters.end(), FilterNameCmp());
    // Add filter specific fields
    for (Iter filt = filters.begin(), eFilt = filters.end(); filt != eFilt; ++filt) {
        schema.addField<int>(
            *filt + ".obs.count",
            "number of " + *filt + "-filter sources in cluster");
        if (!control.exposurePrefix.empty()) {
            schema.addField<double>(
                *filt + ".obs.time.min",
                "earliest observation time of " + *filt + "-filter sources in cluster (TAI)",
                "mjd");
            schema.addField<double>(
                *filt + ".obs.time.max",
                "latest observation time of " + *filt + "-filter sources in cluster (TAI)",
                "mjd");
        }
        for (Iter flux = control.fluxFields.begin(), eFlux = control.fluxFields.end();
             flux != eFlux; ++flux) {
            if (hasFluxField(prototype.getSchema(), *flux)) {
                 addFlux(schema, prototype.getSchema(), *filt, *flux, control.fluxUnit);
            }
        }
        for (Iter shape = control.shapeFields.begin(), eShape = control.shapeFields.end();
             shape != eShape; ++shape) {
            if (hasShapeField(prototype.getSchema(), *shape)) {
                 addShape(schema, prototype.getSchema(), *filt, *shape);
            }
        }
    }

    // Create table
    boost::shared_ptr<SourceClusterTable> table =
        SourceClusterTable::make(schema, idFactory);

    // setup slot mappings
    table->defineCoordErr("coord.err");
    table->defineNumSources("obs.count");
    if (prototype.getCentroidKey().isValid() &&
        prototype.getCentroidErrKey().isValid()) {
        table->defineWeightedMeanCoord("coord.weightedmean");
        table->defineWeightedMeanCoordErr("coord.weightedmean.err");
        table->defineWeightedMeanCoordCount("coord.weightedmean.count");
    }
    if (!control.exposurePrefix.empty()) {
        table->defineTimeMin("obs.time.min");
        table->defineTimeMean("obs.time.mean");
        table->defineTimeMax("obs.time.max");
    }
    for (Iter filt = filters.begin(), eFilt = filters.end();
         filt != eFilt; ++filt) {
        table->defineNumSources(*filt, "obs.count");
        if (!control.exposurePrefix.empty()) {
            table->defineTimeMin(*filt, "obs.time.min");
            table->defineTimeMax(*filt, "obs.time.max");
        }
        Iter flux = control.fluxFields.begin(), eFlux = control.fluxFields.end();
        if (prototype.getPsfFluxKey().isValid() &&
            prototype.getPsfFluxErrKey().isValid()) {
            std::string def = prototype.getPsfFluxDefinition();
            if (std::find(flux, eFlux, def) != eFlux) {
                table->definePsfFlux(*filt, def);
            }
        }
        if (prototype.getModelFluxKey().isValid() &&
            prototype.getModelFluxErrKey().isValid()) {
            std::string def = prototype.getModelFluxDefinition();
            if (std::find(flux, eFlux, def) != eFlux) {
                table->defineModelFlux(*filt, def);
            }
        }
        if (prototype.getApFluxKey().isValid() &&
            prototype.getApFluxErrKey().isValid()) {
            std::string def = prototype.getApFluxDefinition();
            if (std::find(flux, eFlux, def) != eFlux) {
                table->defineApFlux(*filt, def);
            }
        }
        if (prototype.getInstFluxKey().isValid() &&
            prototype.getInstFluxErrKey().isValid()) {
            std::string def = prototype.getInstFluxDefinition();
            if (std::find(flux, eFlux, def) != eFlux) {
                table->defineInstFlux(*filt, def);
            }
        }
        if (prototype.getShapeKey().isValid() &&
            prototype.getShapeErrKey().isValid()) {
            std::string def = prototype.getShapeDefinition();
            Iter shape = control.shapeFields.begin(), eShape = control.shapeFields.end();
            if (std::find(shape, eShape, def) != eShape) {
                table->defineShape(*filt, def);
            }
        }
    }
    return table;
}

void lsst::ap::cluster::processSources	(	lsst::afw::table::SourceCatalog const &	expSources,
		lsst::ap::match::ExposureInfo const &	expInfo,
		lsst::ap::utils::PT1SkyTile const *	skyTile,
		SourceProcessingControl const &	control,
		lsst::afw::table::SchemaMapper const &	mapper,
		lsst::afw::table::SourceCatalog &	sources,
		lsst::afw::table::SourceCatalog &	badSources,
		lsst::afw::table::SourceCatalog &	invalidSources
	)

Process input sourcees, distributing them to one of 3 output catalogs.

• sources: sources suitable for spatial clustering.
• badSources: sources ignored for the purposes of clustering; identified by one or more flags.
• invalidSources: sources which do not have sky-coordinates or centroids.

Sources not lying in the given sky-tile are discarded. Note that one can pass in a NULL sky-tile pointer, in which case no sources are discarded.

Existing fields in the input sources are modified as follows:

• source footprints are discarded
• source sky-coordinate errors are computed and added
• exposure information (ID, filter ID, etc...) is added to each source.
• the cluster coordinates for a source are set to the source coordinates.

The output catalogs will typically have been constructed from tables obtained via makeOutputSourceTable() - their schemas and slot mappings must all be identical. The input table slots must match output table slots, and the input schema must be fully contained in the output schema.

Parameters

[in]	expSources	Single exposure sources to process.
[in]	expInfo	Exposure information.
[in]	skyTile	Sky-tile being processed - may be NULL.
[in]	control	Source processing parameters.
[in]	mapper	Maps between input and output source records.
[out]	sources	Catalog for sources that will be clustered.
[out]	badSources	Catalog for sources with bad measurement flags.
[out]	invalidSources	Catalog for sources with invalid measurements.

void lsst::ap::cluster::processSources	(	SourceCatalog const &	expSources,
		ExposureInfo const &	expInfo,
		lsst::ap::utils::PT1SkyTile const *	skyTile,
		SourceProcessingControl const &	control,
		SchemaMapper const &	mapper,
		SourceCatalog &	sources,
		SourceCatalog &	badSources,
		SourceCatalog &	invalidSources
	)

Definition at line 375 of file clustering.cc.

{
    typedef std::vector<std::string>::const_iterator StringIter;
    typedef std::vector<Key<Flag> >::const_iterator FlagKeyIter;
    typedef SourceCatalog::const_iterator SourceIter;

    Log log(Log::getDefaultLog(), "lsst.ap.cluster");
    log.format(Log::INFO, "Processing sources from exposure %lld",
               static_cast<long long>(expInfo.getId()));

    // Validate that schemas and slots match as expected
    if (sources.getSchema() != badSources.getSchema() ||
        sources.getSchema() != invalidSources.getSchema()) {
        throw LSST_EXCEPT(InvalidParameterError, "output "
            "SourceCatalog schema mismatch");
    }
    if (!sources.getSchema().contains(expSources.getSchema())) {
        throw LSST_EXCEPT(InvalidParameterError, "output SourceCatalog "
            "schema does not contain input SourceCatalog schema");
    }
    if (mapper.getInputSchema() != expSources.getSchema()) {
        throw LSST_EXCEPT(InvalidParameterError, "SchemaMapper and "
            "input source catalog disagree on input schema.");
    }
    if (mapper.getOutputSchema() != sources.getSchema()) {
        throw LSST_EXCEPT(InvalidParameterError, "SchemaMapper and "
            "output source catalogs disagree on output schema.");
    }
    if (!compareSlots(*expSources.getTable(), *sources.getTable()) ||
        !compareSlots(*sources.getTable(), *badSources.getTable()) ||
        !compareSlots(*badSources.getTable(), *invalidSources.getTable())) {
        throw LSST_EXCEPT(InvalidParameterError, "slot mappings "
            "for input and output SourceCatalogs differ");
    }

    // Construct bad source flag key vector
    std::vector<Key<Flag> > badKeys;
    for (StringIter i = control.badFlagFields.begin(),
                    e = control.badFlagFields.end(); i != e; ++i) {
        badKeys.push_back(expSources.getSchema().find<Flag>(*i).key);
    }
    // Extract exposure X/Y bounds relative to parent exposure (if there is one).
    // This is because centroids of coadd-sources are in the pixel coordinate system
    // of the tract exposure containing the patch they were measured on, NOT in
    // the pixel coordinate system of the patch.
    double const xMin = indexToPosition(expInfo.getX0()) - 0.5;
    double const yMin = indexToPosition(expInfo.getY0()) - 0.5;
    double const xMax = indexToPosition(expInfo.getX0() + expInfo.getWidth() - 1) + 0.5;
    double const yMax = indexToPosition(expInfo.getY0() + expInfo.getHeight() - 1) + 0.5;
    
    // Set up keys to additional fields
    Key<int64_t> expIdKey;
    Key<int> expFilterIdKey;
    Key<float> expTimeKey;
    Key<double> expTimeMidKey;
    afw::table::CovarianceMatrixKey<float,2> centroidErrKey =
        sources.getTable()->getCentroidErrKey();
    Key<Covariance<lsst::afw::table::Point<float> > > coordErrKey;
    Key<lsst::afw::coord::Coord> clusterCoordKey;
    try {
        coordErrKey = sources.getSchema()["coord.err"];
    } catch (NotFoundError &) {
        // no easy way to ask whether a field exists by name?
    }
    if (!control.exposurePrefix.empty()) {
        Schema schema = mapper.getOutputSchema();
        expIdKey = schema[control.exposurePrefix + ".id"];
        if (!control.multiBand) {
            expFilterIdKey = schema[control.exposurePrefix + ".filter.id"];
        }
        if (!control.coadd) {
            expTimeKey = schema[control.exposurePrefix + ".time"];
            expTimeMidKey = schema[control.exposurePrefix + ".time.mid"];
        }
    }
    if (!control.clusterPrefix.empty()) {
        Schema schema = mapper.getOutputSchema();
        clusterCoordKey = schema[control.clusterPrefix + ".coord"];
    }
    size_t ngood = 0, nbad = 0, ninvalid = 0, noutside = 0;
    // Loop over input sources
    for (SourceIter s = expSources.begin(), es = expSources.end(); s != es; ++s) {
        // Check source validity
        bool invalid = false;
        double x = s->getX();
        double y = s->getY();
        if (lsst::utils::isnan(x) || lsst::utils::isnan(y)) {
            log.format(Log::WARN, "Centroid of source %lld contains NaNs",
                       static_cast<long long>(s->getId()));
            invalid = true;
        } else if (x < xMin || x > xMax || y < yMin || y > yMax) {
            log.format(Log::WARN, "Centroid of source %lld does not lie on exposure",
                       static_cast<long long>(s->getId()));
            invalid = true;
        } else if (lsst::utils::isnan(s->getRa().asRadians()) ||
                   lsst::utils::isnan(s->getDec().asRadians())) {
            log.format(Log::WARN, "Sky-coord of source %lld contains NaNs",
                       static_cast<long long>(s->getId()));
            invalid = true;
        }
        if (!invalid && skyTile && !skyTile->contains(s->getCoord())) {
            // skip sources outside the sky-tilea
            noutside += 1;
            continue;
        }
        bool bad = false;
        // Check whether source is flagged as bad
        for (FlagKeyIter f = badKeys.begin(), ef = badKeys.end(); f != ef; ++f) {
            if (s->get(*f)) {
                bad = true;
                break;
            }
        }
        // Copy input source to appropriate output catalog
        boost::shared_ptr<SourceRecord> os;
        if (invalid) {
            ninvalid += 1;
            os = invalidSources.addNew();
        } else if (bad) {
            nbad += 1;
            os = badSources.addNew();
        } else {
            ngood += 1;
            os = sources.addNew();
        }
        os->assign(*s, mapper);
        // Add exposure parameters
        if (expIdKey.isValid()) {
            os->set(expIdKey, expInfo.getId());
        }
        if (expFilterIdKey.isValid()) {
            os->set(expFilterIdKey, expInfo.getFilter().getId());
        }
        if (expTimeKey.isValid()) {
            os->set(expTimeKey, expInfo.getExposureTime());
        }
        if (expTimeMidKey.isValid()) {
            os->set(expTimeMidKey, expInfo.getEpoch());
        }
        // Source cluster coords default to source coordinates
        if (clusterCoordKey.isValid()) {
            os->set(clusterCoordKey, s->getCoord());
        }
        // Compute sky-coordinate errors
        if (!invalid && coordErrKey.isValid() && centroidErrKey.isValid()) {
            Eigen::Matrix2d cov = s->getCentroidErr().cast<double>();
            bool computeErr = true;
            if (lsst::utils::isnan(cov(0,0)) || lsst::utils::isnan(cov(1,1))) {
                computeErr = false;
            } else if (cov(0,1) != cov(1,0)) {
                // FIXME: for now, no measurement algorithms actually
                //        compute full covariance matrixes, and the
                //        off diagonal matrix elements are always NaN.
                //        Longer term, we will need some algorithm metadata
                //        to decide whether an off-diagonal NaN means a
                //        sample should be ignored because a computation failed,
                //        or whether it should be zeroed because the algorithm
                //        never computes it.
                if (lsst::utils::isnan(cov(0,1)) && lsst::utils::isnan(cov(1,0))) {
                    cov(0,1) = 0.0; cov(1,0) = 0.0;
                } else {
                    computeErr = false;
                }
            }
            if (computeErr) {
                Eigen::Matrix2d m = expInfo.getWcs()->linearizePixelToSky(
                    s->getCentroid(), radians).getLinear().getMatrix();
                os->set(coordErrKey, (m * cov * m.transpose()).cast<float>());
            }
        }
    }
    log.format(Log::INFO, "processed %lld sources (invalid: %lld, "
               "outside sky-tile: %lld, bad: %lld, good: %lld)",
               static_cast<long long>(expSources.size()),
               static_cast<long long>(ninvalid),
               static_cast<long long>(noutside),
               static_cast<long long>(nbad),
               static_cast<long long>(ngood));
}

void lsst::ap::cluster::setClusterFields	(	lsst::afw::table::SourceCatalog &	sources,
		SourceClusterRecord const &	record,
		SourceProcessingControl const &	control
	)

Set the "<cluster>.id" and "<cluster>.coord" fields of each source in the given catalog to the ID and sky-coordinates of the given cluster. The "<cluster>" field name prefix is obtained from SourceProcessingControl. If those fields are not setup in the given catalog, this function is a no-op.

This is intended to support database ingest via qserv, where sources are partitioned by their associated astrophysical object, and objects (for which clusters are a temporary stand-in) are partitioned by position. Denormalizing the output source schema by appending cluster position avoids a potentially very expensive join during database ingest.

Parameters

[out]	sources	Sources to update.
[in]	record	Cluster to obtain ID/sky-coordinates from.
[in]	control	Supplies cluster field name prefix.

void lsst::ap::cluster::setClusterFields	(	SourceCatalog &	sources,
		SourceClusterRecord const &	record,
		SourceProcessingControl const &	control
	)

Definition at line 615 of file clustering.cc.

{
    typedef SourceCatalog::iterator Iter;

    if (control.clusterPrefix.empty()) {
        return;
    }
    Key<int64_t> idKey;
    Key<lsst::afw::coord::Coord> coordKey;
    try {
        idKey = sources.getSchema().find<int64_t>(control.clusterPrefix + ".id").key;
        coordKey = sources.getSchema().find<lsst::afw::coord::Coord>(
            control.clusterPrefix + ".coord").key;
    } catch (NotFoundError &) {
        // could not find fields with the expected name
        return;
    }
    int64_t const id = record.getId();
    IcrsCoord const coord = record.getCoord();
    for (Iter i = sources.begin(), e = sources.end(); i != e; ++i) {
        i->set(idKey, id);
        i->set(coordKey, coord);
    }
}