LSSTApplications  1.1.2+25,10.0+13,10.0+132,10.0+133,10.0+224,10.0+41,10.0+8,10.0-1-g0f53050+14,10.0-1-g4b7b172+19,10.0-1-g61a5bae+98,10.0-1-g7408a83+3,10.0-1-gc1e0f5a+19,10.0-1-gdb4482e+14,10.0-11-g3947115+2,10.0-12-g8719d8b+2,10.0-15-ga3f480f+1,10.0-2-g4f67435,10.0-2-gcb4bc6c+26,10.0-28-gf7f57a9+1,10.0-3-g1bbe32c+14,10.0-3-g5b46d21,10.0-4-g027f45f+5,10.0-4-g86f66b5+2,10.0-4-gc4fccf3+24,10.0-40-g4349866+2,10.0-5-g766159b,10.0-5-gca2295e+25,10.0-6-g462a451+1
LSSTDataManagementBasePackage
Namespaces | Classes | Typedefs | Functions
lsst::ap Namespace Reference

Namespaces

 cluster
 
 detail
 @ endcond
 
 io
 
 match
 
 pipeline
 
 SqlAppendTables
 
 SqlDropTables
 
 SqlStoreOutputs
 
 tasks
 
 utils
 
 version
 

Classes

class  Bitset
 A fixed size set of bits. More...
 
struct  BinChunkHeader
 Simple header for binary chunk files – allows some sanity checking at read time. More...
 
class  ChunkDescriptor
 A generic descriptor containing state for different kinds of chunks. More...
 
class  SharedObjectChunkManager
 A manager for Object chunks that exist in shared memory. More...
 
class  CircularRegion
 A circular region of the unit sphere (sky). More...
 
class  Condition
 Encapsulates a POSIX condition variable. More...
 
struct  DataTraits
 Provides basic chunk parameters at compile time. More...
 
struct  DataTraits< Object >
 
class  Ellipse
 Contains spatial information for a single ellipse on the unit sphere (sky). More...
 
struct  EllipsePtrLessThan
 Comparison functor for Ellipse pointers that orders ellipses by minimum overlapping zone. More...
 
class  EllipseList
 A list of ellipses, implemented using std::vector. More...
 
class  Fifo
 A First In, First Out (FIFO) queue of fixed capacity. More...
 
struct  PassthroughFilter
 A default "let everything through" filter implementation. More...
 
struct  MatchWithDistance
 Contains a pointer to a match and an associated distance. More...
 
struct  MatchWithoutDistance
 Contains a pointer to a match. More...
 
struct  EmptyMatchListProcessor
 A default "do nothing" match list processing implementation. More...
 
struct  EmptyMatchPairProcessor
 A default "do nothing" match pair processing implementation. More...
 
class  ScopedLock
 Grants access to a mutex, enforcing the RAII principle. More...
 
class  Mutex
 A wrapper for a process private POSIX mutual exclusion lock. More...
 
class  SharedMutex
 A wrapper for a POSIX process shared mutual exclusion lock. More...
 
struct  Object
 A partial representation of a full LSST Object containing only id, position, proper motions, and per-filter variability probabilities. More...
 
struct  Point
 A point on the unit sphere (sky), specified in spherical polar coordinates. More...
 
class  RectangularRegion
 A rectangular region (in right ascension and declination) of the unit sphere. More...
 
class  MatchPair
 Holds a pair of ids and the distance between the corresponding positions on the unit sphere. More...
 
class  PersistableMatchPairVector
 A persistable wrapper for a MatchPairVector. More...
 
class  PersistableIdPairVector
 A persistable wrapper for an IdPairVector. More...
 
class  PersistableIdVector
 A persistable wrapper for an IdVector. More...
 
class  ScopeGuard
 Utility class for automatically invoking a function when leaving a scope. More...
 
class  ZoneStripeChunkDecomposition
 A decomposition of the unit sphere into zones, stripes, and chunks. More...
 
struct  DiaSourceChunk
 
class  VisitProcessingContext
 Container for inter-stage association pipeline state. More...
 
class  TimeSpec
 Wraps the C library timespec struct. More...
 
class  Stopwatch
 Utility class for profiling. More...
 
struct  ZoneEntry
 Contains spatial information for a single point used during cross-matching. More...
 
struct  ZoneEntryArray
 Stores entries inside a single zone (a narrow declination stripe) in a sorted array. More...
 
class  ZoneIndex
 Container for a sequence of adjacent zones. More...
 

Typedefs

typedef unsigned char ChunkEntryFlag
 
typedef std::pair
< boost::int64_t,
boost::int64_t > 		IdPair
 Holds a pair of ids. More...
 
typedef std::vector< MatchPair > MatchPairVector
 A list of MatchPair instances. More...
 
typedef std::vector< IdPair > IdPairVector
 A list of IdPair instances. More...
 
typedef std::vector
< boost::int64_t > 		IdVector
 A list of integer ids. More...
 

Functions

void initialize (std::string const &runId)
 
void registerVisit (VisitProcessingContext &context)
 
void loadSliceObjects (VisitProcessingContext &context)
 
void buildObjectIndex (VisitProcessingContext &context)
 
void matchDiaSources (MatchPairVector &matches, VisitProcessingContext &context)
 
void matchMops (MatchPairVector &matches, IdPairVector &newObjects, VisitProcessingContext &context, lsst::mops::MovingObjectPredictionVector &predictions)
 
void storeSliceObjects (VisitProcessingContext &context)
 
void failVisit (VisitProcessingContext &context)
 
bool endVisit (VisitProcessingContext &context, bool const rollback)
 
std::ostream & operator<< (std::ostream &os, Stopwatch const &watch)
 
template<typename DataT >
void swap (Ellipse< DataT > &a, Ellipse< DataT > &b)
 
template<typename DataT >
bool operator< (Ellipse< DataT > const &a, Ellipse< DataT > const &b)
 
template<typename DataT >
bool operator== (Ellipse< DataT > const &a, Ellipse< DataT > const &b)
 
template<typename DataT >
bool operator< (int32_t const a, Ellipse< DataT > const &b)
 
template<typename DataT >
bool operator< (Ellipse< DataT > const &a, int32_t const b)
 
template<typename DataT >
bool operator== (int32_t const a, Ellipse< DataT > const &b)
 
template<typename DataT >
bool operator== (Ellipse< DataT > const &a, int32_t const b)
 
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchListProcessorT >
std::size_t distanceMatch (ZoneIndex< FirstEntryT > &first, ZoneIndex< SecondEntryT > &second, double const radius, FirstFilterT &firstFilter, SecondFilterT &secondFilter, MatchListProcessorT &matchListProcessor)
 
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchPairProcessorT >
std::size_t ellipseMatch (EllipseList< FirstEntryT > &first, ZoneIndex< SecondEntryT > &second, FirstFilterT &firstFilter, SecondFilterT &secondFilter, MatchPairProcessorT &matchPairProcessor)
 
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchListProcessorT >
std::size_t ellipseGroupedMatch (EllipseList< FirstEntryT > &first, ZoneIndex< SecondEntryT > &second, FirstFilterT &firstFilter, SecondFilterT &secondFilter, MatchListProcessorT &matchListProcessor)
 
bool operator== (Object const &o1, Object const &o2)
 
bool operator!= (Object const &o1, Object const &o2)
 
void swap (ZoneStripeChunkDecomposition &a, ZoneStripeChunkDecomposition &b)
 
double alpha (double const theta, double const centerDec, double const dec)
 
double maxAlpha (double const theta, double const centerDec)
 
void computeChunkIds (std::vector< int > &chunkIds, CircularRegion const &region, ZoneStripeChunkDecomposition const &zsc, int const workerId=0, int const numWorkers=1)
 
void computeChunkIds (std::vector< int > &chunkIds, RectangularRegion const &region, ZoneStripeChunkDecomposition const &zsc, int const workerId=0, int const numWorkers=1)
 
boost::uint32_t raToScaledInteger (double const ra)
 
boost::uint32_t deltaRaToScaledInteger (double const delta)
 
boost::int32_t decToScaledInteger (double const dec)
 
boost::int32_t deltaDecToScaledInteger (double const delta)
 
double clampDec (double const dec)
 
void verifyPathName (std::string const &name)
 
template<typename ChunkT >
bool operator< (ZoneEntry< ChunkT > const &a, ZoneEntry< ChunkT > const &b)
 
template<typename ChunkT >
bool operator== (ZoneEntry< ChunkT > const &a, ZoneEntry< ChunkT > const &b)
 
template<typename ChunkT >
bool operator< (boost::uint32_t const a, ZoneEntry< ChunkT > const &b)
 
template<typename ChunkT >
bool operator< (ZoneEntry< ChunkT > const &a, boost::uint32_t const b)
 
template<typename ChunkT >
bool operator== (boost::uint32_t const a, ZoneEntry< ChunkT > const &b)
 
template<typename ChunkT >
bool operator== (ZoneEntry< ChunkT > const &a, boost::uint32_t const b)
 

Typedef Documentation

typedef unsigned char lsst::ap::ChunkEntryFlag

Definition at line 142 of file Chunk.h.

typedef std::pair<boost::int64_t, boost::int64_t> lsst::ap::IdPair

Holds a pair of ids.

Definition at line 121 of file Results.h.

typedef std::vector<IdPair> lsst::ap::IdPairVector

A list of IdPair instances.

Definition at line 125 of file Results.h.

typedef std::vector<boost::int64_t> lsst::ap::IdVector

A list of integer ids.

Definition at line 127 of file Results.h.

typedef std::vector<MatchPair> lsst::ap::MatchPairVector

A list of MatchPair instances.

Definition at line 123 of file Results.h.

Function Documentation

double lsst::ap::alpha ( double const  theta,
double const  centerDec,
double const  dec 
)

Intersects the plane given by z = sin(dec) with the circle of radius theta and center (0, centerDec) on the unit sphere, then returns the right ascension alpha giving the two resulting points: (-alpha, dec) and (alpha, dec).

Precondition
theta > 0.0 && theta < 10.0
centerDec >= -90.0 && centerDec <= 90.0
dec >= -90.0 && dec <= 90.0
Parameters
[in]thetathe radius of the input circle
[in]centerDecthe declination of the circle center
[in]decthe declination of the horizontal plane to intersect with the input circle
Returns
alpha, the largest right ascension of the two points in the intersection of the input circle with the input plane.

Definition at line 248 of file SpatialUtil.cc.

252  {
253  assert(theta > 0.0 && theta < 10.0 && "radius out of range");
254  assert(centerDec >= -90.0 && centerDec <= 90.0 && "declination out of range");
255  assert(dec >= -90.0 && dec <= 90.0 && "declination out of range");
256 
257  if (std::fabs(centerDec) + theta > 89.9) {
258  return 180.0;
259  }
260  double x = std::cos(afwGeom::degToRad(theta)) - std::sin(afwGeom::degToRad(centerDec))*std::sin(afwGeom::degToRad(dec));
261  double u = std::cos(afwGeom::degToRad(centerDec))*std::cos(afwGeom::degToRad(dec));
262  double y = std::sqrt(std::fabs(u*u - x*x));
263  return afwGeom::radToDeg(std::fabs(std::atan2(y,x)));
264 }
y
int y
Definition: GaussianCentroid.cc:36
radToDeg
double radToDeg(long double angleInRadians)
Definition: RaDecStr.cc:61
x
int x
Definition: GaussianCentroid.cc:36
degToRad
double degToRad(long double angleInDegrees)
Definition: RaDecStr.cc:67
void lsst::ap::buildObjectIndex ( VisitProcessingContext &  context)

Checks to see if all parallel workers succeeded in loading their chunks and, if so, builds a zone index for the objects just loaded.

Parameters
[in,out]contextState involved in processing a single visit.

Definition at line 775 of file Stages.cc.

775  {
776  if (!context.debugSharedMemory()) {
777  // if the shared memory object used for chunk storage hasn't yet been unlinked, do so now
778  SharedObjectChunkManager::destroyInstance(context.getRunId());
779  }
780  SharedObjectChunkManager manager(context.getRunId());
781  if (manager.isVisitInFlight(context.getVisitId())) {
782  try {
783  // Build zone index on objects
784  manager.getChunks(context.getChunks(), context.getChunkIds());
785  context.buildObjectIndex();
786  } catch(...) {
787  manager.endVisit(context.getVisitId(), true);
788  throw;
789  }
790  } else {
791  // One or more workers failed in the load phase - rollback the visit
792  manager.endVisit(context.getVisitId(), true);
793  throw LSST_EXCEPT(ex::RuntimeError,
794  "Association pipeline failed to read Object data for FOV");
795  }
796 }
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
double lsst::ap::clampDec ( double const  dec)
inline

Clamps the given declination value to [-90,90].

Definition at line 276 of file SpatialUtil.h.

276  {
277  if (dec <= -90.0) {
278  return -90.0;
279  } else if (dec >= 90.0) {
280  return 90.0;
281  }
282  return dec;
283 }
void lsst::ap::computeChunkIds ( std::vector< int > &  chunkIds,
CircularRegion const &  region,
ZoneStripeChunkDecomposition const &  zsc,
int const  workerId = 0,
int const  numWorkers = 1 
)

Computes identifiers for all chunks in the given ZoneStripeChunkDecomposition that overlap the given region and belong to the specified worker. Chunks belonging to a stripe s such that the euclidian remainder of s/numWorkers is workerId belong to the worker identified by workerId.

Parameters
[out]chunkIdsThe list in which to store the computed chunk identifiers.
[in]regionThe region for which overlapping chunks are to be computed.
[in]zscA decomposition of the unit sphere into stripes, chunks, and zones.
[in]workerIdThe integer id of the current worker (in a set of numWorkers parallel workers).
[in]numWorkersThe number of parallel workers.

Definition at line 308 of file SpatialUtil.cc.

314  {
315  if (numWorkers < 1) {
316  throw LSST_EXCEPT(ex::InvalidParameterError, "number of workers must be positive");
317  }
318  if (workerId < 0 || workerId >= numWorkers) {
319  throw LSST_EXCEPT(ex::InvalidParameterError,
320  "Worker id must be between 0 and N - 1, where N is the total number of workers");
321  }
322 
323  for (int s = zsc.decToStripe(region.getMinDec()); s <= zsc.decToStripe(region.getMaxDec()); ++s) {
324 
325  // round-robin stripes to workers
326  int rem = s % numWorkers;
327  if (rem < 0) {
328  rem += numWorkers;
329  }
330  if (rem != workerId) {
331  continue;
332  }
333 
334  int const fc = zsc.getFirstChunkForStripe(s);
335  int const nc = zsc.getNumChunksPerStripe(s);
336  double a = zsc.stripeAndCircleToRaRange(s,
337  region.getCenterRa(), region.getCenterDec(), region.getRadius());
338  double raMin = region.getCenterRa() - a;
339  double raMax = region.getCenterRa() + a;
340  bool wrap = false;
341 
342  if (raMin < 0.0) {
343  raMin += 360.0;
344  wrap = true;
345  }
346  if (raMax >= 360.0) {
347  raMax -= 360.0;
348  wrap = true;
349  }
350  int maxChunk = static_cast<int>(std::floor((raMax/360.0)*nc));
351  if (maxChunk == nc) {
352  --maxChunk;
353  }
354  int chunk = static_cast<int>(std::floor((raMin/360.0)*nc));
355  if (chunk == nc) {
356  --chunk;
357  }
358  chunk += fc;
359  maxChunk += fc;
360 
361  if (raMax < raMin || (raMax == raMin && wrap)) {
362  if (chunk == maxChunk) {
363  --maxChunk; // avoid adding the same chunk twice
364  }
365  for ( ; chunk < fc + nc; ++chunk) {
366  chunkIds.push_back(chunk);
367  }
368  for (chunk = 0; chunk <= maxChunk; ++chunk) {
369  chunkIds.push_back(chunk);
370  }
371  } else {
372  for ( ; chunk <= maxChunk; ++chunk) {
373  chunkIds.push_back(chunk);
374  }
375  }
376  }
377 }
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
lsst.afw.geom::floor
Extent< int, N > floor(Extent< double, N > const &input)
void lsst::ap::computeChunkIds ( std::vector< int > &  chunkIds,
RectangularRegion const &  region,
ZoneStripeChunkDecomposition const &  zsc,
int const  workerId = 0,
int const  numWorkers = 1 
)

Computes identifiers for all chunks in the given ZoneStripeChunkDecomposition that overlap the given region and belong to the specified worker. Chunks belonging to a stripe s such that the euclidian remainder of s/numWorkers is workerId belong to the worker identified by workerId.

Parameters
[out]chunkIdsThe list in which to store the computed chunk identifiers.
[in]regionThe region for which overlapping chunks are to be computed.
[in]zscA decomposition of the unit sphere into stripes, chunks, and zones.
[in]workerIdThe integer id of the current worker (in a set of numWorkers parallel workers).
[in]numWorkersThe number of parallel workers.

Definition at line 393 of file SpatialUtil.cc.

399  {
400  if (numWorkers < 1) {
401  throw LSST_EXCEPT(ex::InvalidParameterError,
402  "number of workers must be positive");
403  }
404  if (workerId < 0 || workerId >= numWorkers) {
405  throw LSST_EXCEPT(ex::InvalidParameterError,
406  "Worker id must be between 0 and N - 1, where N is the total number of workers");
407  }
408 
409  double const raMin = region.getMinRa();
410  double const raMax = region.getMaxRa();
411 
412  for (int s = zsc.decToStripe(region.getMinDec()); s <= zsc.decToStripe(region.getMaxDec()); ++s) {
413 
414  // round-robin stripes to workers
415  int rem = s % numWorkers;
416  if (rem < 0) {
417  rem += numWorkers;
418  }
419  if (rem != workerId) {
420  continue;
421  }
422 
423  int const fc = zsc.getFirstChunkForStripe(s);
424  int const nc = zsc.getNumChunksPerStripe(s);
425 
426  int maxChunk = fc + static_cast<int>(std::floor((raMax*nc)/360.0)) % nc;
427  int chunk = fc + static_cast<int>(std::floor((raMin*nc)/360.0)) % nc;
428 
429  if (raMax < raMin) {
430  if (chunk == maxChunk) {
431  --maxChunk; // avoid adding the same chunk twice
432  }
433  for ( ; chunk < fc + nc; ++chunk) {
434  chunkIds.push_back(chunk);
435  }
436  for (chunk = 0; chunk <= maxChunk; ++chunk) {
437  chunkIds.push_back(chunk);
438  }
439  } else {
440  for (; chunk <= maxChunk; ++chunk) {
441  chunkIds.push_back(chunk);
442  }
443  }
444  }
445 }
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
lsst.afw.geom::floor
Extent< int, N > floor(Extent< double, N > const &input)
boost::int32_t lsst::ap::decToScaledInteger ( double const  dec)
inline

Converts a declination (in degrees) to an integer between -2^30 and 2^30

Definition at line 264 of file SpatialUtil.h.

264  {
265  assert(dec >= -90.0 && dec <= 90.0);
266  return static_cast<int32_t>(std::floor(dec*RA_DEC_SCALE));
267 }
lsst.afw.geom::floor
Extent< int, N > floor(Extent< double, N > const &input)
boost::int32_t lsst::ap::deltaDecToScaledInteger ( double const  delta)
inline

Converts a declination delta (in degrees) to an integer

Definition at line 270 of file SpatialUtil.h.

270  {
271  assert(delta >= 0.0 && delta <= 90.0);
272  return static_cast<int32_t>(std::ceil(delta*RA_DEC_SCALE));
273 }
lsst.afw.geom::ceil
Extent< int, N > ceil(Extent< double, N > const &input)
boost::uint32_t lsst::ap::deltaRaToScaledInteger ( double const  delta)
inline

Converts a right ascension delta (in degrees) to an integer

Definition at line 258 of file SpatialUtil.h.

258  {
259  assert(delta >= 0.0 && delta <= 180.0);
260  return static_cast<int32_t>(std::ceil(delta*RA_DEC_SCALE));
261 }
lsst.afw.geom::ceil
Extent< int, N > ceil(Extent< double, N > const &input)
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchListProcessorT >
std::size_t lsst::ap::distanceMatch ( ZoneIndex< FirstEntryT > &  first,
ZoneIndex< SecondEntryT > &  second,
double const  radius,
FirstFilterT &  firstFilter,
SecondFilterT &  secondFilter,
MatchListProcessorT &  matchListProcessor 
)

Spatial cross-match routine – finds match pairs in a first and a second set of entities (both subject to filtering), where both sets consist of points. An entity from the second set is deemed a match to an entity from the first set if the two are within the given angle of eachother. All matches for a given entity in the first set are found at once, and sent off to a match list processor for further inspection.

This routine is optimized for the case where few matches are expected for any given entity.

Precondition
radius >= 0
Parameters
[in]firstA first set of entities.
[in]secondA second set of entities.
[in]radiusThe match-radius (in degrees).
[in]firstFilterA filter on the first set of entities.
[in]secondFilterA filter on the second set of entities.
[in]matchListProcessorA processor for match lists.
Returns
The number of match pairs found.

Definition at line 153 of file Match.h.

160  {
161  typedef typename ZoneIndex<FirstEntryT>::Zone FirstZone;
162  typedef typename ZoneIndex<SecondEntryT>::Zone SecondZone;
163  typedef typename MatchListProcessorT::Match Match;
164 
165  if (radius < 0) {
166  throw LSST_EXCEPT(lsst::pex::exceptions::RangeError,
167  "match radius must be greater than or equal to zero degrees");
168  }
169 
170  double const shr = std::sin(radians(radius*0.5));
171  double const d2Limit = 4.0*shr*shr;
172  int const minZone = first.getMinZone();
173  int const maxZone = first.getMaxZone();
174  boost::int32_t const deltaDec = deltaDecToScaledInteger(radius);
175 
176  std::size_t numMatchPairs = 0;
177 
178  second.computeMatchParams(radius);
179 
180  // loop over first set of zones in parallel
181 #if LSST_AP_HAVE_OPEN_MP
182 # pragma omp parallel default(shared)
183 #endif
184  {
185  // allocate per-thread data structures
186  SecondZone * zones[2048];
187  int limits[2048];
188  std::vector<Match> matches;
189  matches.reserve(32);
190 
191  // loop over the first set of zones in parallel, assigning batches of
192  // adjacent zones to threads (for cache coherency)
193 #if LSST_AP_HAVE_OPEN_MP
194 # pragma omp for \
195  reduction(+:numMatchPairs) \
196  schedule(static,8)
197 #endif
198  for (int fzi = minZone; fzi <= maxZone; ++fzi) {
199 
200  FirstZone * const __restrict fz = first.getZone(fzi);
201  FirstEntryT * const __restrict fze = fz->_entries;
202  int const nfze = fz->_size;
203  if (nfze <= 0) {
204  continue; // no entries in first zone
205  }
206 
207  // populate secondary zone array with potentially matching zones
208  int nsz = 0;
209  {
210  double d = first.getDecomposition().getZoneDecMin(fz->_zone) - radius;
211  int const minz = second.getDecomposition().decToZone(d <= -90.0 ? -90.0 : d);
212  d = first.getDecomposition().getZoneDecMax(fz->_zone) + radius;
213  int const maxz = second.getDecomposition().decToZone(d >= 90.0 ? 90.0 : d);
214  // A search circle should never cover more than 2048 zones
215  assert(maxz - minz + 1 <= 2048 && "match radius too large");
216 
217  SecondZone * __restrict sz = second.firstZone(minz, maxz);
218  SecondZone * const __restrict szend = second.endZone(minz, maxz);
219 
220  for ( ; sz < szend; ++sz) {
221  int const nsze = sz->_size;
222  if (nsze > 0) {
223  zones[nsz] = sz;
224  if ((nsze >> 4) > nfze) {
225  // second set much larger than first, use binary search in inner loop
226  limits[nsz] = -1;
227  } else {
228  // use linear walk in inner loop, find starting point now
229  limits[nsz] = sz->findGte(fze[0]._ra - sz->_deltaRa);
230  }
231  ++nsz;
232  }
233  }
234  }
235 
236  if (nsz == 0) {
237  // no entries in any potentially matching zones
238  continue;
239  }
240 
241  // loop over entries in first zone
242  for (int fe = 0; fe < nfze; ++fe) {
243 
244  if (!firstFilter(fze[fe])) {
245  continue; // entry was filtered out
246  }
247 
248  matches.clear();
249 
250  boost::uint32_t const ra = fze[fe]._ra;
251  boost::int32_t const dec = fze[fe]._dec;
252  double const fx = fze[fe]._x;
253  double const fy = fze[fe]._y;
254  double const fz = fze[fe]._z;
255 
256  // loop over all potentially matching zones.
257  for (int szi = 0; szi < nsz; ++szi) {
258 
259  SecondZone * const __restrict sz = zones[szi];
260  SecondEntryT * const __restrict sze = sz->_entries;
261 
262  int const seWrap = sz->_size;
263  boost::uint32_t const deltaRa = sz->_deltaRa;
264  boost::uint32_t const deltaRaWrap = -deltaRa;
265 
266  int start = limits[szi];
267  int se;
268  boost::uint32_t dra;
269 
270  if (start >= 0) {
271 
272  // use linear walk from last starting point
273  // to get to first point in ra range
274  se = start;
275  dra = ra - sze[se]._ra;
276  if (dra > deltaRa && dra < deltaRaWrap) {
277 
278  bool cont = false;
279  do {
280  ++se;
281  if (se == seWrap) {
282  se = 0; // ra wrap around
283  }
284  if (se == start) {
285  cont = true;
286  break; // avoid infinite loops
287  }
288  dra = ra - sze[se]._ra;
289  } while (dra > deltaRa && dra < deltaRaWrap);
290 
291  if (cont) {
292  continue; // no starting point found
293  }
294  // found starting point -- remember it
295  limits[szi] = se;
296  start = se;
297  }
298 
299  } else {
300 
301  // use binary search to find starting point
302  start = sz->findGte(ra - deltaRa);
303  se = start;
304  dra = ra - sze[se]._ra;
305  if (dra > deltaRa && dra < deltaRaWrap) {
306  continue; // no starting point found
307  }
308 
309  }
310 
311  // At this point, entry se is within ra range of fe --
312  // loop over all zone entries within ra range
313  do {
314 
315  // test whether entry se is within dec range
316  if (abs(dec - sze[se]._dec) <= deltaDec) {
317  // yes -- perform detailed distance test
318  double xd = (fx - sze[se]._x);
319  double yd = (fy - sze[se]._y);
320  double zd = (fz - sze[se]._z);
321  double d2 = xd*xd + yd*yd + zd*zd;
322  if (d2 < d2Limit) {
323  // Note: this isn't necessarily the best place for the second filter test...
324  if (secondFilter(sze[se])) {
325  // found a match, record it
326  matches.push_back(Match(&sze[se], d2));
327  }
328  }
329  }
330  ++se;
331  if (se == seWrap) {
332  se = 0; // ra wrap around
333  }
334  if (se == start) {
335  break; // avoid infinite loops
336  }
337  dra = ra - sze[se]._ra;
338 
339  } while (dra <= deltaRa || dra >= deltaRaWrap);
340 
341  } // end of loop over potentially matching zones
342 
343  // All matches (if any) for fe are found
344  std::size_t nm = matches.size();
345  if (nm > 0) {
346  // pass them on to the match processor
347  numMatchPairs = numMatchPairs + nm;
348  matchListProcessor(fze[fe], matches.begin(), matches.end());
349  }
350 
351  } // end of loop over entries in first zone
352 
353  } // end of omp for
354  } // end of omp parallel
355 
356  return numMatchPairs;
357 }
lsst::ap::deltaDecToScaledInteger
boost::int32_t deltaDecToScaledInteger(double const delta)
Definition: SpatialUtil.h:270
lsst.afw.geom::radians
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
d
int d
Definition: KDTree.cc:89
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchListProcessorT >
std::size_t lsst::ap::ellipseGroupedMatch ( EllipseList< FirstEntryT > &  first,
ZoneIndex< SecondEntryT > &  second,
FirstFilterT &  firstFilter,
SecondFilterT &  secondFilter,
MatchListProcessorT &  matchListProcessor 
)

Spatial cross-match routine – finds match pairs in a first and a second set of entities (both subject to filtering), where the first set consists of ellipses and the second of points. An entity in the second set is deemed a match for an entity in the first set if it is within the ellipse defined by the first entity. All matches for a given ellipse are found at once and sent off to a match list processor for further inspection.

Parameters
[in]firstA first set of entities.
[in]secondA second set of entities.
[in]firstFilterA filter on the first set of entities.
[in]secondFilterA filter on the second set of entities.
[in]matchListProcessorA processor for match lists.
Returns
The number of match pairs found.

Definition at line 535 of file Match.h.

541  {
542  typedef typename ZoneIndex<SecondEntryT>::Zone SecondZone;
543  typedef typename MatchListProcessorT::Match Match;
544 
545  std::size_t const numEllipses = first.size();
546  std::size_t numMatchPairs = 0;
547 
548  first.prepareForMatch(second.getDecomposition());
549 
550 #if LSST_AP_HAVE_OPEN_MP
551 # pragma omp parallel default(shared)
552 #endif
553  {
554  // allocate per-thread match list
555  std::vector<Match> matches;
556  matches.reserve(2048);
557 
558  // loop over the list of ellipses in parallel
559 #if LSST_AP_HAVE_OPEN_MP
560 # pragma omp for \
561  reduction(+:numMatchPairs) \
562  schedule(static,128)
563 #endif
564  for (std::size_t i = 0; i < numEllipses; ++i) {
565 
566  if (!firstFilter(first[i])) {
567  continue; // ellipse was filtered out
568  }
569 
570  Ellipse<FirstEntryT> * const __restrict ell = &first[i];
571  SecondZone * __restrict sz = second.firstZone(ell->_minZone, ell->_maxZone);
572  SecondZone * const __restrict szend = second.endZone(ell->_minZone, ell->_maxZone);
573 
574  boost::uint32_t const ra = ell->_ra;
575  boost::uint32_t const deltaRa = ell->_deltaRa;
576  boost::uint32_t const deltaRaWrap = -deltaRa;
577 
578  matches.clear();
579 
580  // loop over zones covered by the ellipse
581  for ( ; sz < szend; ++sz) {
582 
583  int const seWrap = sz->_size;
584  if (seWrap == 0) {
585  continue;
586  }
587 
588  // find starting point within ra range of the ellipse
589  SecondEntryT * const __restrict sze = sz->_entries;
590  int const start = sz->findGte(ra - deltaRa);
591  int se = start;
592  boost::uint32_t dra = ra - sze[se]._ra;
593 
594  while (dra <= deltaRa || dra >= deltaRaWrap) {
595 
596  // check whether entry is within dec range
597  boost::int32_t const dec = sze[se]._dec;
598  if (dec >= ell->_minDec && dec <= ell->_maxDec) {
599  // perform detailed in ellipse test
600  if (ell->contains(sze[se]._x, sze[se]._y, sze[se]._z)) {
601  if (secondFilter(sze[se])) {
602  // record match
603  matches.push_back(Match(&sze[se]));
604  }
605  }
606  }
607  ++se;
608  if (se == seWrap) {
609  se = 0; // ra wrap around
610  }
611  if (se == start) {
612  break; // avoid infinite loops
613  }
614  dra = ra - sze[se]._ra;
615  }
616  }
617 
618  // All matches (if any) for ell are found
619  std::size_t nm = matches.size();
620  if (nm > 0) {
621  // pass them on to the match processor
622  numMatchPairs = numMatchPairs + nm;
623  matchListProcessor(*ell, matches.begin(), matches.end());
624  }
625 
626  } // end of omp for
627  } // end of omp parallel
628 
629  return numMatchPairs;
630 }
template<typename FirstEntryT , typename SecondEntryT , typename FirstFilterT , typename SecondFilterT , typename MatchPairProcessorT >
std::size_t lsst::ap::ellipseMatch ( EllipseList< FirstEntryT > &  first,
ZoneIndex< SecondEntryT > &  second,
FirstFilterT &  firstFilter,
SecondFilterT &  secondFilter,
MatchPairProcessorT &  matchPairProcessor 
)

Spatial cross-match routine – finds match pairs in a first and a second set of entities (both subject to filtering), where the first set consists of ellipses and the second of points. An entity in the second set is deemed a match for an entity in the first set if it is within the ellipse defined by the first entity. Note – this particular routine does not package up all matches for a given ellipse before sending them off to a match list processor. Instead, it reports match pairs to a match pair processor in no particular order. The routine is also currently single threaded. For a parallel routine that uses a match list processor, see ellipseGroupedMatch(). The advantage of this routine is that it should have much better memory access patterns (smaller cache footprint, more sequential access) in the presence of ellipses with wildly varying size.

Parameters
[in]firstA first set of entities.
[in]secondA second set of entities.
[in]firstFilterA filter on the first set of entities.
[in]secondFilterA filter on the second set of entities.
[in]matchPairProcessorA processor for match pairs.
Returns
The number of match pairs found.

Definition at line 386 of file Match.h.

392  {
393  typedef typename ZoneIndex<SecondEntryT>::Zone SecondZone;
394 
395  int const minZone = second.getMinZone();
396  int const maxZone = second.getMaxZone();
397  std::size_t numMatchPairs = 0;
398 
399  first.prepareForMatch(second.getDecomposition());
400 
401  Ellipse<FirstEntryT> * activeHead = 0;
402  Ellipse<FirstEntryT> * activeTail = 0;
403  Ellipse<FirstEntryT> * searchHead = &(*first.begin());
404  Ellipse<FirstEntryT> * const end = &(*first.end());
405 
406  // initialize the linked list of active ellipses (those that intersect minZone)
407  while (searchHead < end && searchHead->_minZone <= minZone) {
408  if (searchHead->_maxZone >= minZone) {
409  if (firstFilter(*searchHead)) {
410  if (activeTail == 0) {
411  activeHead = searchHead;
412  } else {
413  activeTail->_next = searchHead;
414  }
415  activeTail = searchHead;
416  }
417  }
418  ++searchHead;
419  }
420 
421  // loop over the set of zones in the second entity set
422  int szi = minZone;
423  while (true) {
424 
425  SecondZone * const __restrict sz = second.getZone(szi);
426  SecondEntryT * const __restrict sze = sz->_entries;
427  int const seWrap = sz->_size;
428 
429  ++szi;
430 
431  // Traverse the active ellipse list
432  Ellipse<FirstEntryT> * active = activeHead;
433  Ellipse<FirstEntryT> * prev = 0;
434 
435  while (active != 0) {
436 
437 #if defined(__GNUC__)
438  __builtin_prefetch(active->_next, 0, 3);
439 #endif
440  if (seWrap > 0) {
441  // find entries within ra range of the active ellipse
442  boost::uint32_t const ra = active->_ra;
443  boost::uint32_t const deltaRa = active->_deltaRa;
444  boost::uint32_t const deltaRaWrap = -deltaRa;
445  int const start = sz->findGte(ra - deltaRa);
446 
447  int se = start;
448  boost::uint32_t dra = ra - sze[se]._ra;
449  while (dra <= deltaRa || dra >= deltaRaWrap) {
450 
451  // check whether entry is within dec range
452  boost::int32_t const dec = sze[se]._dec;
453  if (dec >= active->_minDec && dec <= active->_maxDec) {
454  // perform detailed in ellipse test
455  if (active->contains(sze[se]._x, sze[se]._y, sze[se]._z)) {
456  if (secondFilter(sze[se])) {
457  // process match pair
458  matchPairProcessor(*active, sze[se]);
459  ++numMatchPairs;
460  }
461  }
462  }
463  ++se;
464  if (se == seWrap) {
465  se = 0; // ra wrap around
466  }
467  if (se == start) {
468  break; // avoid infinite loops
469  }
470  dra = ra - sze[se]._ra;
471  }
472  }
473 
474  if (active->_maxZone < szi) {
475  // ellipse no longer active in following zone -- remove it from the active list
476  active = active->_next;
477  if (prev == 0) {
478  activeHead = active;
479  } else {
480  prev->_next = active;
481  }
482  } else {
483  // keep ellipse in the active list
484  prev = active;
485  active = active->_next;
486  }
487 
488  } // end of loop over active ellipses
489 
490  if (szi > maxZone) {
491  break;
492  }
493 
494  // Add ellipses with minimum zone equal to szi (the zone that will be searched
495  // in the next loop iteration) to the active list
496  while (searchHead < end && searchHead->_minZone <= szi) {
497  if (firstFilter(*searchHead)) {
498  if (prev == 0) {
499  activeHead = searchHead;
500  } else {
501  prev->_next = searchHead;
502  }
503  prev = searchHead;
504  }
505  ++searchHead;
506  }
507 
508  } // end of loop over second set of zones
509 
510  return numMatchPairs;
511 }
bool lsst::ap::endVisit ( VisitProcessingContext &  context,
bool const  rollback 
)

Ends visit processing - in memory changes are rolled back if the given visit failed in any way, or if rollback is true .

Parameters
[in,out]contextState involved in processing a single visit.
[in]rollbackIndicates whether or not results for the visit should be rolled back
Returns
true if the visit existed, was not marked as a failure and was committed, false otherwise.

Definition at line 1028 of file Stages.cc.

1028  {
1029  SharedObjectChunkManager manager(context.getRunId());
1030  bool committed = manager.endVisit(context.getVisitId(), rollback);
1031  Log log(Log::getDefaultLog(), "lsst.ap");
1032  if (committed) {
1033  log.log(Log::INFO, "Committed visit");
1034  } else {
1035  log.log(Log::FATAL, "Rolled back visit");
1036  }
1037  return committed;
1038 }
lsst.pex::logging::Log
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
lsst.log.log.FATAL
int FATAL
Definition: log.py:40
lsst.log.log.log
def log
Definition: log.py:85
lsst.log.log.INFO
int INFO
Definition: log.py:37
void lsst::ap::failVisit ( VisitProcessingContext &  context)

Marks processing for the given visit as a failure.

Parameters
[in,out]contextState involved in processing a single visit.

Definition at line 1012 of file Stages.cc.

1012  {
1013  SharedObjectChunkManager manager(context.getRunId());
1014  manager.failVisit(context.getVisitId());
1015 }
void lsst::ap::initialize ( std::string const &  runId)

Sets up all fundamental visit processing parameters using a policy and ensure that a reference to the shared memory object used for chunk storage exists.

Definition at line 644 of file Stages.cc.

644  {
645  // create shared memory object if it doesn't already exist
646  volatile SharedObjectChunkManager manager(runId);
647 }
void lsst::ap::loadSliceObjects ( VisitProcessingContext &  context)

Ensures that object data for the chunks assigned to the calling slice has been read in or is owned by the given visit.

Definition at line 666 of file Stages.cc.

666  {
667 
668  typedef VisitProcessingContext::ObjectChunk Chunk;
669  typedef std::vector<Chunk> ChunkVector;
670  typedef std::vector<Chunk>::iterator ChunkIterator;
671 
672  SharedObjectChunkManager manager(context.getRunId());
673  Log log(Log::getDefaultLog(), "lsst.ap");
674 
675  try {
676 
677  // From FOV, get chunk ids to be handled by this worker
678  Stopwatch watch(true);
679  computeChunkIds(
680  context.getChunkIds(),
681  context.getFov(),
682  context.getDecomposition(),
683  context.getWorkerId(),
684  context.getNumWorkers()
685  );
686  watch.stop();
687  Rec(log, Log::INFO) << "computed chunk ids in FOV for worker " <<
688  Prop<int>("numChunks", static_cast<int>(context.getChunkIds().size())) <<
689  Prop<double>("time", watch.seconds()) << Rec::endr;
690 
691  // Register interest in or create chunks via the chunk manager
692  std::vector<Chunk> toRead;
693  std::vector<Chunk> toWaitFor;
694  watch.start();
695  manager.startVisit(toRead, toWaitFor, context.getVisitId(), context.getChunkIds());
696  watch.stop();
697  Rec(log, Log::INFO) << "started processing visit" <<
698  Prop<double>("time", watch.seconds()) << Rec::endr;
699 
700  // record pointers to all chunks being handled by the slice
701  ChunkVector & chunks = context.getChunks();
702  chunks.clear();
703  chunks.insert(chunks.end(), toRead.begin(), toRead.end());
704  chunks.insert(chunks.end(), toWaitFor.begin(), toWaitFor.end());
705 
706  // Read data files
707  watch.start();
708  std::string refNamePattern(context.getPipelinePolicy()->getString("objectChunkFileNamePattern"));
709  std::string deltaNamePattern(context.getPipelinePolicy()->getString("objectDeltaChunkFileNamePattern"));
710  PropertySet::Ptr ps(new PropertySet);
711  ps->set<std::string>("runId", context.getRunId());
712  ChunkVector::size_type numToRead(toRead.size());
713  for (ChunkIterator i(toRead.begin()), end(toRead.end()); i != end; ++i) {
714  Chunk & c = *i;
715  ps->set<int>("chunkId", c.getId());
716  ps->set<int>("stripeId", ZoneStripeChunkDecomposition::chunkToStripe(c.getId()));
717  ps->set<int>("chunkSeqNum", ZoneStripeChunkDecomposition::chunkToSequence(c.getId()));
718  c.read(LogicalLocation(refNamePattern, ps).locString(), false);
719  c.readDelta(LogicalLocation(deltaNamePattern, ps).locString(), false);
720  c.setUsable();
721  }
722  watch.stop();
723  Rec(log, Log::INFO) << "read chunk files" <<
724  Prop<int>("numChunks", static_cast<int>(numToRead)) <<
725  Prop<double>("time", watch.seconds()) << Rec::endr;
726 
727  toRead.clear();
728  ChunkVector::size_type numToWaitFor = toWaitFor.size();
729  if (numToWaitFor > 0) {
730  watch.start();
731  // Wait for chunks that are owned by another visit
732  manager.waitForOwnership(toRead, toWaitFor, context.getVisitId(), context.getDeadline());
733  watch.stop();
734  Rec(log, Log::INFO) << "acquired ownership of pre-existing chunks" <<
735  Prop<int>("numChunks", static_cast<int>(numToWaitFor)) <<
736  Prop<double>("time", watch.seconds()) << Rec::endr;
737 
738  // Read in chunks that were not successfully read by the previous owner
739  watch.start();
740  numToRead = toRead.size();
741  for (ChunkIterator i(toRead.begin()), end(toRead.end()); i != end; ++i) {
742  Chunk & c = *i;
743  ps->set<int>("chunkId", c.getId());
744  ps->set<int>("stripeId", ZoneStripeChunkDecomposition::chunkToStripe(c.getId()));
745  ps->set<int>("chunkSeqNum", ZoneStripeChunkDecomposition::chunkToSequence(c.getId()));
746  c.read(LogicalLocation(refNamePattern, ps).locString(), false);
747  c.readDelta(LogicalLocation(deltaNamePattern, ps).locString(), false);
748  c.setUsable();
749  }
750  watch.stop();
751  Rec(log, Log::INFO) << "read straggling chunks" <<
752  Prop<int>("numChunks", static_cast<int>(numToRead)) <<
753  Prop<double>("time", watch.seconds()) << Rec::endr;
754  }
755 
756  } catch (ex::Exception & except) {
757  Rec(log, Log::FATAL) << except.what() << Rec::endr;
758  manager.failVisit(context.getVisitId());
759  } catch (std::exception & except) {
760  log.log(Log::FATAL, except.what());
761  manager.failVisit(context.getVisitId());
762  } catch (...) {
763  log.log(Log::FATAL, "caught unknown exception");
764  manager.failVisit(context.getVisitId());
765  }
766 }
lsst::ap::computeChunkIds
void computeChunkIds(std::vector< int > &chunkIds, CircularRegion const &region, ZoneStripeChunkDecomposition const &zsc, int const workerId=0, int const numWorkers=1)
Definition: SpatialUtil.cc:308
lsst::ap::utils::endr
CsvWriter & endr(CsvWriter &w)
Definition: Csv.cc:330
lsst.pex::logging::Rec
LogRec Rec
Definition: Log.h:838
lsst.daf.persistence::LogicalLocation
Class for logical location of a persisted Persistable instance.
Definition: LogicalLocation.h:60
lsst.pex::logging::Log
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
lsst.log.log.FATAL
int FATAL
Definition: log.py:40
lsst.log.log.log
def log
Definition: log.py:85
lsst::ap::detail::ObjectChunk
SharedObjectChunkManager::ObjectChunk ObjectChunk
Definition: Stages.cc:111
lsst.pex.exceptions::Exception::what
virtual char const * what(void) const
Definition: Exception.cc:112
lsst.log.log.INFO
int INFO
Definition: log.py:37
lsst.daf.base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:82
lsst.daf.base::PropertySet::Ptr
boost::shared_ptr< PropertySet > Ptr
Definition: PropertySet.h:90
void lsst::ap::matchDiaSources ( MatchPairVector &  matches,
VisitProcessingContext &  context 
)

Matches difference sources for a visit (obtained from the detection pipeline) against the objects in the visit FOV.

Parameters
[out]matchesSet to a list of difference source to object match pairs.
[in,out]contextState involved in processing a single visit.

Definition at line 808 of file Stages.cc.

811  {
812  SharedObjectChunkManager manager(context.getRunId());
813 
814  try {
815 
816  matches.clear();
817  matches.reserve(65536);
818 
819  detail::ObjectMatchProcessor<detail::ObjectEntry> mlp(context, matches, context.getFilter());
820  PassthroughFilter<detail::DiaSourceEntry> pdf;
821  PassthroughFilter<detail::ObjectEntry> pof;
822 
823  Stopwatch watch(true);
824  std::size_t nm = distanceMatch<
825  detail::DiaSourceEntry,
826  detail::ObjectEntry,
827  PassthroughFilter<detail::DiaSourceEntry>,
828  PassthroughFilter<detail::ObjectEntry>,
829  detail::ObjectMatchProcessor<detail::ObjectEntry>
830  >(
831  context.getDiaSourceIndex(),
832  context.getObjectIndex(),
833  context.getMatchRadius()/3600.0, // match routine expects degrees, not arc-seconds
834  pdf,
835  pof,
836  mlp
837  );
838  watch.stop();
839  Log log(Log::getDefaultLog(), "lsst.ap");
840  Rec(log, Log::INFO) << "matched difference sources to objects" <<
841  Prop<int>("numDiaSources", context.getDiaSourceIndex().size()) <<
842  Prop<int>("numObjects", context.getObjectIndex().size()) <<
843  Prop<int>("numMatches", static_cast<int>(nm)) <<
844  Prop<double>("time", watch.seconds()) << Rec::endr;
845 
846  } catch (...) {
847  manager.endVisit(context.getVisitId(), true);
848  throw;
849  }
850 }
lsst::ap::detail::DiaSourceEntry
ZoneEntry< DiaSourceChunk > DiaSourceEntry
Definition: Stages.cc:116
lsst::ap::utils::endr
CsvWriter & endr(CsvWriter &w)
Definition: Csv.cc:330
lsst.pex::logging::Rec
LogRec Rec
Definition: Log.h:838
lsst::ap::distanceMatch
std::size_t distanceMatch(ZoneIndex< FirstEntryT > &first, ZoneIndex< SecondEntryT > &second, double const radius, FirstFilterT &firstFilter, SecondFilterT &secondFilter, MatchListProcessorT &matchListProcessor)
Definition: Match.h:153
lsst.pex::logging::Log
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
lsst.log.log.log
def log
Definition: log.py:85
lsst::ap::detail::ObjectEntry
ZoneEntry< ObjectChunk > ObjectEntry
Definition: Stages.cc:115
lsst.log.log.INFO
int INFO
Definition: log.py:37
void lsst::ap::matchMops ( MatchPairVector &  matches,
IdPairVector &  newObjects,
VisitProcessingContext &  context,
lsst::mops::MovingObjectPredictionVector &  predictions 
)

Matches moving object predictions falling within the FOV of a visit against the difference sources for that visit.

Parameters
[out]matchesSet to a list of moving object prediction to difference source match pairs.
[out]newObjectsSet to a list of (difference source id, object id) pairs that specifies which difference sources should be used to create new object and what the id of each new object should be set to.
[in,out]contextState involved in processing a single visit.
[in]predictionsThe list of moving object predictions to match against difference sources.

Definition at line 866 of file Stages.cc.

871  {
872  SharedObjectChunkManager manager(context.getRunId());
873 
874  try {
875  double const ellScale = context.getPipelinePolicy()->getDouble("ellipseScalingFactor");
876  double const smaaClamp = context.getPipelinePolicy()->getDouble("semiMajorAxisClamp");
877  double const smiaClamp = context.getPipelinePolicy()->getDouble("semiMinorAxisClamp");
878 
879  matches.clear();
880  matches.reserve(8192);
881  newObjects.clear();
882  newObjects.reserve(8192);
883 
884  detail::MovingObjectPredictionMatchProcessor mpp(matches);
885 
886  // discard difference sources with known variable matches
887  Stopwatch watch(true);
888  detail::DiscardKnownVariableFilter dvf;
889  int nr = context.getDiaSourceIndex().pack(dvf);
890  watch.stop();
891  Log log(Log::getDefaultLog(), "lsst.ap");
892  Rec(log, Log::INFO) << "removed difference sources matching known variables from index" <<
893  Prop<int>("numRemoved", nr) <<
894  Prop<double>("time", watch.seconds()) << Rec::endr;
895 
896  // build ellipses required for matching from predictions
897  watch.start();
898  EllipseList<lsst::mops::MovingObjectPrediction::MovingObjectPrediction> ellipses;
899  ellipses.reserve(predictions.size());
900  detail::DiscardLargeEllipseFilter elf(context.getPipelinePolicy());
901  lsst::mops::MovingObjectPredictionVector::iterator const end = predictions.end();
902  for (lsst::mops::MovingObjectPredictionVector::iterator i = predictions.begin(); i != end; ++i) {
903  i->setSemiMajorAxisLength(i->getSemiMajorAxisLength() * ellScale);
904  i->setSemiMinorAxisLength(i->getSemiMinorAxisLength() * ellScale);
905  if (elf(*i)) {
906  // clamp error ellipses if necessary
907  if (smaaClamp > 0.0 && i->getSemiMajorAxisLength() > smaaClamp) {
908  i->setSemiMajorAxisLength(smaaClamp);
909  }
910  if (smiaClamp > 0.0 && i->getSemiMinorAxisLength() > smiaClamp) {
911  i->setSemiMinorAxisLength(smiaClamp);
912  }
913  ellipses.push_back(Ellipse<lsst::mops::MovingObjectPrediction>(*i));
914  }
915  }
916  watch.stop();
917  Rec(log, Log::INFO) << "built list of match parameters for moving object predictions" <<
918  Prop<int>("numPredictions", static_cast<int>(ellipses.size())) <<
919  Prop<double>("time", watch.seconds()) << Rec::endr;
920 
921  // match them against difference sources
922  PassthroughFilter<detail::MovingObjectEllipse> pef;
923  PassthroughFilter<detail::DiaSourceEntry> pdf;
924  watch.start();
925  std::size_t nm = ellipseMatch<
926  MovingObjectPrediction,
927  detail::DiaSourceEntry,
928  PassthroughFilter<detail::MovingObjectEllipse>,
929  PassthroughFilter<detail::DiaSourceEntry>,
930  detail::MovingObjectPredictionMatchProcessor
931  >(
932  ellipses,
933  context.getDiaSourceIndex(),
934  pef,
935  pdf,
936  mpp
937  );
938  watch.stop();
939  Rec(log, Log::INFO) << "matched moving object predictions to difference sources" <<
940  Prop<int>("numPredictions", static_cast<int>(ellipses.size())) <<
941  Prop<int>("numDiaSources", context.getDiaSourceIndex().size()) <<
942  Prop<int>("numMatches", static_cast<int>(nm)) <<
943  Prop<double>("time", watch.seconds()) << Rec::endr;
944 
945  // Create new objects from difference sources with no matches
946  watch.start();
947  detail::NewObjectCreator createObjects(newObjects, context);
948  context.getDiaSourceIndex().apply(createObjects);
949  watch.stop();
950  Rec(log, Log::INFO) << "created new objects" <<
951  Prop<int>("numObjects", static_cast<int>(newObjects.size())) <<
952  Prop<double>("time", watch.seconds()) << Rec::endr;
953  } catch (...) {
954  manager.endVisit(context.getVisitId(), true);
955  throw;
956  }
957 }
lsst::ap::detail::DiaSourceEntry
ZoneEntry< DiaSourceChunk > DiaSourceEntry
Definition: Stages.cc:116
lsst::ap::utils::endr
CsvWriter & endr(CsvWriter &w)
Definition: Csv.cc:330
lsst.pex::logging::Rec
LogRec Rec
Definition: Log.h:838
lsst::ap::ellipseMatch
std::size_t ellipseMatch(EllipseList< FirstEntryT > &first, ZoneIndex< SecondEntryT > &second, FirstFilterT &firstFilter, SecondFilterT &secondFilter, MatchPairProcessorT &matchPairProcessor)
Definition: Match.h:386
lsst.pex::logging::Log
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
lsst.log.log.log
def log
Definition: log.py:85
lsst.log.log.INFO
int INFO
Definition: log.py:37
double lsst::ap::maxAlpha ( double const  theta,
double const  centerDec 
)

Computes the extent in right ascension [-alpha, alpha] of the circle with radius theta and center (0, centerDec) on the unit sphere.

Precondition
theta > 0.0 && theta < 10.0
centerDec >= -90.0 && centerDec <= 90.0
Parameters
[in]thetathe radius of the circle to find ra extents for
[in]centerDecthe declination of the circle center (in degrees)
Returns
the largest right ascension of any point on the input circle

Definition at line 279 of file SpatialUtil.cc.

282  {
283  assert(theta > 0.0 && theta < 10.0 && "radius out of range");
284  assert(centerDec >= -90.0 && centerDec <= 90.0 && "declination out of range");
285 
286  if (std::fabs(centerDec) + theta > 89.9) {
287  return 180.0;
288  }
289  double y = std::sin(afwGeom::degToRad(theta));
290  double x = std::sqrt(std::fabs(cos(afwGeom::degToRad(centerDec - theta))*cos(afwGeom::degToRad(centerDec + theta))));
291  return afwGeom::radToDeg(std::fabs(std::atan(y/x)));
292 }
y
int y
Definition: GaussianCentroid.cc:36
radToDeg
double radToDeg(long double angleInRadians)
Definition: RaDecStr.cc:61
x
int x
Definition: GaussianCentroid.cc:36
degToRad
double degToRad(long double angleInDegrees)
Definition: RaDecStr.cc:67
bool lsst::ap::operator!= ( Object const &  o1,
Object const &  o2 
)
inline

Definition at line 112 of file Object.h.

112  {
113  return !(o1 == o2);
114 }
template<typename ChunkT >
bool lsst::ap::operator< ( ZoneEntry< ChunkT > const &  a,
ZoneEntry< ChunkT > const &  b 
)
inline

Definition at line 78 of file ZoneTypes.h.

78  {
79  return a._ra < b._ra;
80 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename ChunkT >
bool lsst::ap::operator< ( boost::uint32_t const  a,
ZoneEntry< ChunkT > const &  b 
)
inline

Definition at line 88 of file ZoneTypes.h.

88  {
89  return a < b._ra;
90 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename ChunkT >
bool lsst::ap::operator< ( ZoneEntry< ChunkT > const &  a,
boost::uint32_t const  b 
)
inline

Definition at line 93 of file ZoneTypes.h.

93  {
94  return a._ra < b;
95 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename DataT >
bool lsst::ap::operator< ( Ellipse< DataT > const &  a,
Ellipse< DataT > const &  b 
)
inline

Definition at line 95 of file EllipseTypes.h.

95  {
96  return a._minZone < b._minZone;
97 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename DataT >
bool lsst::ap::operator< ( int32_t const  a,
Ellipse< DataT > const &  b 
)
inline

Definition at line 105 of file EllipseTypes.h.

105  {
106  return a < b._minZone;
107 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename DataT >
bool lsst::ap::operator< ( Ellipse< DataT > const &  a,
int32_t const  b 
)
inline

Definition at line 110 of file EllipseTypes.h.

110  {
111  return a._minZone < b;
112 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
std::ostream& lsst::ap::operator<< ( std::ostream &  os,
Stopwatch const &  watch 
)

Definition at line 168 of file Time.cc.

168  {
169  TimeSpec t;
170  if (watch._stopped) {
171  t = watch._ts;
172  } else {
173  t.now();
174  t -= watch._ts;
175  }
176  long nsec = t.tv_nsec % 1000000000l;
177  long sec = t.tv_sec + t.tv_nsec/1000000000l;
178  long min = (sec/60) % 60;
179  long hour = sec/3600;
180 
181  double s = static_cast<double>(sec % 60) + 1e-9 * static_cast<double>(nsec);
182 
183  if (hour > 0) {
184  os << hour << "hr ";
185  }
186  if (min > 0) {
187  os << min << "min ";
188  }
189  os << s << "sec";
190  return os;
191 }
min
double min
Definition: attributes.cc:216
template<typename ChunkT >
bool lsst::ap::operator== ( ZoneEntry< ChunkT > const &  a,
ZoneEntry< ChunkT > const &  b 
)
inline

Definition at line 83 of file ZoneTypes.h.

83  {
84  return a._ra == b._ra;
85 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename ChunkT >
bool lsst::ap::operator== ( boost::uint32_t const  a,
ZoneEntry< ChunkT > const &  b 
)
inline

Definition at line 98 of file ZoneTypes.h.

98  {
99  return a == b._ra;
100 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename DataT >
bool lsst::ap::operator== ( Ellipse< DataT > const &  a,
Ellipse< DataT > const &  b 
)
inline

Definition at line 100 of file EllipseTypes.h.

100  {
101  return a._minZone == b._minZone;
102 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename ChunkT >
bool lsst::ap::operator== ( ZoneEntry< ChunkT > const &  a,
boost::uint32_t const  b 
)
inline

Definition at line 103 of file ZoneTypes.h.

103  {
104  return a._ra == b;
105 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
bool lsst::ap::operator== ( Object const &  o1,
Object const &  o2 
)

Definition at line 36 of file Object.cc.

36  {
37  if (o1._objectId == o2._objectId && o1._ra == o2._ra && o1._decl == o2._decl) {
38  for (int i = 0; i < lsst::ap::Object::NUM_FILTERS; ++i) {
39  if (o1._varProb[i] != o2._varProb[i]) {
40  return false;
41  }
42  }
43  return true;
44  }
45  return false;
46 }
lsst::ap::Object::NUM_FILTERS
static int const NUM_FILTERS
Definition: Object.h:54
template<typename DataT >
bool lsst::ap::operator== ( int32_t const  a,
Ellipse< DataT > const &  b 
)
inline

Definition at line 115 of file EllipseTypes.h.

115  {
116  return a == b._minZone;
117 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
template<typename DataT >
bool lsst::ap::operator== ( Ellipse< DataT > const &  a,
int32_t const  b 
)
inline

Definition at line 120 of file EllipseTypes.h.

120  {
121  return a._minZone == b;
122 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
boost::uint32_t lsst::ap::raToScaledInteger ( double const  ra)
inline

Converts a right ascension (in degrees) to an integer in the range [0, 2^32)

Definition at line 251 of file SpatialUtil.h.

251  {
252  assert(ra >= 0.0 && ra < 360.0);
253  double d = std::floor(ra*RA_DEC_SCALE);
254  return d >= 4294967295.0 ? 4294967295u : static_cast<uint32_t>(d);
255 }
d
int d
Definition: KDTree.cc:89
lsst.afw.geom::floor
Extent< int, N > floor(Extent< double, N > const &input)
void lsst::ap::registerVisit ( VisitProcessingContext &  context)

Computes ids for all object chunks covering the visit FOV and registers the visit with the shared memory chunk manager.

Definition at line 654 of file Stages.cc.

654  {
655  context.getChunkIds().clear();
656  computeChunkIds(context.getChunkIds(), context.getFov(), context.getDecomposition(), 0, 1);
657  SharedObjectChunkManager manager(context.getRunId());
658  manager.registerVisit(context.getVisitId());
659 }
lsst::ap::computeChunkIds
void computeChunkIds(std::vector< int > &chunkIds, CircularRegion const &region, ZoneStripeChunkDecomposition const &zsc, int const workerId=0, int const numWorkers=1)
Definition: SpatialUtil.cc:308
void lsst::ap::storeSliceObjects ( VisitProcessingContext &  context)

Stores any new objects that have been added to the FOV of the visit.

Parameters
[in,out]contextState involved in processing a single visit.

Definition at line 967 of file Stages.cc.

967  {
968 
969  typedef VisitProcessingContext::ObjectChunk Chunk;
970  typedef std::vector<Chunk> ChunkVector;
971  typedef std::vector<Chunk>::iterator ChunkIterator;
972 
973  SharedObjectChunkManager manager(context.getRunId());
974  Log log(Log::getDefaultLog(), "lsst.ap");
975  try {
976  Stopwatch watch(true);
977  std::string deltaNamePattern = context.getPipelinePolicy()->getString("objectDeltaChunkFileNamePattern");
978  PropertySet::Ptr ps(new PropertySet);
979  ps->set<std::string>("runId", context.getRunId());
980  ChunkVector & chunks = context.getChunks();
981  for (ChunkIterator i(chunks.begin()), end(chunks.end()); i != end; ++i) {
982  Chunk & c = *i;
983  ps->set<int>("chunkId", c.getId());
984  ps->set<int>("stripeId", ZoneStripeChunkDecomposition::chunkToStripe(c.getId()));
985  ps->set<int>("chunkSeqNum", ZoneStripeChunkDecomposition::chunkToSequence(c.getId()));
986  std::string file = LogicalLocation(deltaNamePattern, ps).locString();
987  verifyPathName(file);
988  c.writeDelta(file, true, false);
989  }
990  watch.stop();
991  Rec(log, Log::INFO) << "wrote chunk delta files" <<
992  Prop<int>("numChunks", static_cast<int>(chunks.size())) <<
993  Prop<double>("time", watch.seconds()) << Rec::endr;
994  } catch (ex::Exception & except) {
995  Rec(log, Log::FATAL) << except.what() << Rec::endr;
996  manager.failVisit(context.getVisitId());
997  } catch (std::exception & except) {
998  log.log(Log::FATAL, except.what());
999  manager.failVisit(context.getVisitId());
1000  } catch (...) {
1001  log.log(Log::FATAL, "caught unknown exception");
1002  manager.failVisit(context.getVisitId());
1003  }
1004 }
lsst::ap::utils::endr
CsvWriter & endr(CsvWriter &w)
Definition: Csv.cc:330
lsst.pex::logging::Rec
LogRec Rec
Definition: Log.h:838
lsst.daf.persistence::LogicalLocation
Class for logical location of a persisted Persistable instance.
Definition: LogicalLocation.h:60
lsst.daf.persistence::LogicalLocation::locString
std::string const & locString(void) const
Definition: LogicalLocation.cc:132
lsst.pex::logging::Log
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
lsst.log.log.FATAL
int FATAL
Definition: log.py:40
lsst.log.log.log
def log
Definition: log.py:85
lsst::ap::detail::ObjectChunk
SharedObjectChunkManager::ObjectChunk ObjectChunk
Definition: Stages.cc:111
lsst::ap::verifyPathName
void verifyPathName(std::string const &name)
Definition: Utils.cc:49
lsst.pex.exceptions::Exception::what
virtual char const * what(void) const
Definition: Exception.cc:112
lsst.log.log.INFO
int INFO
Definition: log.py:37
lsst.daf.base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:82
lsst.daf.base::PropertySet::Ptr
boost::shared_ptr< PropertySet > Ptr
Definition: PropertySet.h:90
template<typename DataT >
void lsst::ap::swap ( Ellipse< DataT > &  a,
Ellipse< DataT > &  b 
)
inline

Definition at line 90 of file EllipseTypes.h.

90  {
91  a.swap(b);
92 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
void lsst::ap::swap ( ZoneStripeChunkDecomposition &  a,
ZoneStripeChunkDecomposition &  b 
)
inline

Definition at line 225 of file SpatialUtil.h.

225  {
226  a.swap(b);
227 }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:48
void lsst::ap::verifyPathName ( std::string const &  name)

Ensure that all directories along a path exist, creating them if necessary.

Parameters
[in]namePathname to file to be created

Definition at line 49 of file Utils.cc.

49  {
50  // Get the directory by stripping off anything after the last slash.
51  std::string::size_type pos = name.find_last_of('/');
52  if (pos == std::string::npos) return;
53  std::string dirName = name.substr(0, pos);
54 
55  // Check to see if the directory exists.
56  struct stat buf;
57  int ret = ::stat(dirName.c_str(), &buf);
58 
59  if (ret == -1 && errno == ENOENT) {
60  // It doesn't; check its parent and then create it.
61  verifyPathName(dirName);
62 
63  ret = ::mkdir(dirName.c_str(), 0777);
64  if (ret == -1) {
65  throw LSST_EXCEPT(ex::IoError, "Error creating directory " + dirName);
66  }
67  }
68  else if (ret == -1) {
69  // We couldn't read the (existing) directory for some reason.
70  throw LSST_EXCEPT(ex::IoError, "Error searching for directory " + dirName);
71  }
72  else if (!S_ISDIR(buf.st_mode)) {
73  // It's not a directory.
74  throw LSST_EXCEPT(ex::IoError, "Non-directory in path " + dirName);
75  }
76 }
lsst::ap::verifyPathName
void verifyPathName(std::string const &name)
Definition: Utils.cc:49
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Generated on Tue Apr 21 2015 10:05:04 for LSSTApplications by   doxygen 1.8.5