lsst::ap::detail Namespace Reference

@ endcond More...

Classes
struct	BootstrapLock
	Mutual exclusion lock suitable for initializing shared memory objects. More...
class	ObjectMatchProcessor
	Processor for lists of objects matching a difference source. More...
struct	MovingObjectPredictionMatchProcessor
	Processor for matches between moving object predictions and difference sources. More...
struct	DiscardKnownVariableFilter
	Filter which discards difference sources matching known variable objects. More...
struct	DiscardLargeEllipseFilter
	Filter which discards predicted moving objects with large position error ellipses. More...
struct	NewObjectCreator
	Records ids of difference sources with no matches. More...
struct	BitTraits
struct	BitTraits< uint8_t >
struct	BitTraits< uint16_t >
struct	BitTraits< uint32_t >
struct	BitTraits< uint64_t >
class	HashedSet
	A set of up to NumEntries elements of type EntryT, hashed by an integer identifier. More...
class	BlockAllocator
	A thread-safe memory block allocator that uses a Bitset to track which blocks (out of a fixed size pool of blocks) are in-use/free. More...
class	Visit
	State for a single visit to a field of view. More...
class	VisitTracker
	Tracks a set of visits. More...
class	SubManager
	Helper class for managing chunks of a particular type. More...
class	ChunkManagerImpl
	A manager for a set of chunks of a single type. More...

Typedefs
typedef BlockAllocator < SharedMutex, Object >	ObjAllocator
typedef ChunkManagerImpl < SharedMutex, Object >	ObjChunkMgr
typedef SharedObjectChunkManager::ObjectChunk	ObjectChunk
typedef std::vector< ObjectChunk >	ObjectChunkVector
typedef ZoneEntry< ObjectChunk >	ObjectEntry
typedef ZoneEntry< DiaSourceChunk >	DiaSourceEntry
typedef Ellipse < MovingObjectPrediction >	MovingObjectEllipse

Functions
template<typename ManagerT >
ManagerT *	getSingleton (char const *const shmObjName, char const *const shmLockName)
std::pair< double, double >	correctProperMotion (Object const &obj, double const epoch)
template<typename EntryT >
void	buildZoneIndex (ZoneIndex< EntryT > &index, std::vector< typename EntryT::Chunk > const &chunks, double const epoch)
template<typename WordT >
int	wordForBit (int const i)
template<typename WordT >
WordT	maskForBit (int const i)
template<typename WordT >
bool	setBits (int *const indexes, WordT *const words, int const numBitsToSet, int const numBits)
template<typename WordT >
void	resetBits (WordT *const words, int const *const indexes, int const numBitsToReset, int const numBits)
int	hash (boost::uint32_t key)
int	hash (int key)

Detailed Description

@ endcond

Typedef Documentation

typedef ZoneEntry<DiaSourceChunk> lsst::ap::detail::DiaSourceEntry

Definition at line 116 of file Stages.cc.

typedef Ellipse<MovingObjectPrediction> lsst::ap::detail::MovingObjectEllipse

Definition at line 118 of file Stages.cc.

typedef BlockAllocator<SharedMutex, Object> lsst::ap::detail::ObjAllocator

Definition at line 59 of file ChunkManager.cc.

typedef ChunkManagerImpl<SharedMutex, Object> lsst::ap::detail::ObjChunkMgr

Definition at line 67 of file ChunkManager.cc.

typedef SharedObjectChunkManager::ObjectChunk lsst::ap::detail::ObjectChunk

Definition at line 111 of file Stages.cc.

typedef std::vector<ObjectChunk> lsst::ap::detail::ObjectChunkVector

Definition at line 113 of file Stages.cc.

typedef ZoneEntry<ObjectChunk> lsst::ap::detail::ObjectEntry

Definition at line 115 of file Stages.cc.

Function Documentation

template<typename EntryT >

void lsst::ap::detail::buildZoneIndex	(	ZoneIndex< EntryT > &	index,
		std::vector< typename EntryT::Chunk > const &	chunks,
		double const	epoch
	)

Definition at line 383 of file Stages.cc.

  {
    typedef typename EntryT::Data Data;
    typedef typename EntryT::Chunk Chunk;
    typedef std::vector<Chunk> ChunkVector;
    typedef typename ChunkVector::const_iterator ChunkIterator;
    typedef typename ChunkVector::size_type Size;

    index.clear();
    if (chunks.empty()) {
        return;
    }

    Stopwatch watch(true);

    // determine stripe bounds for the input chunks
    ZoneStripeChunkDecomposition const & zsc = index.getDecomposition();
    int minStripe = 0x7FFFFFFF;
    int maxStripe = -1 - minStripe;
    ChunkIterator const end(chunks.end());
    for (ChunkIterator c(chunks.begin()); c != end; ++c) {
        int const stripeId = ZoneStripeChunkDecomposition::chunkToStripe(c->getId());
        if (stripeId > maxStripe) {
            maxStripe = stripeId;
        }
        if (stripeId < minStripe) {
            minStripe = stripeId;
        }
    }
    assert(maxStripe >= minStripe && "invalid stripe bounds for chunk list");

    // Partition input chunks into stripes
    int const numStripes = maxStripe - minStripe + 1;
    boost::scoped_array<ChunkVector> stripes(new ChunkVector[numStripes]);
    for (ChunkIterator c(chunks.begin()); c != end; ++c) {
        int const stripeId = ZoneStripeChunkDecomposition::chunkToStripe(c->getId());
        assert(stripeId >= minStripe && stripeId <= maxStripe && "stripe id out of bounds");
        stripes[stripeId - minStripe].push_back(*c);
    }

    double minDec = zsc.getStripeDecMin(minStripe) - 0.001;
    double maxDec = zsc.getStripeDecMax(maxStripe) + 0.001;
    index.setDecBounds(std::max(minDec, -90.0), std::min(maxDec, 90.0));

    // Loop over stripes. Note that due to proper motion, objects from
    // different stripes can end up in the same zone. Objects are never
    // migrated across chunks/stripes; rather, their J2000 coordinates
    // determine the chunks that they will belong to. As time goes on,
    // the spatial extent of a chunk will grow, where the rate of growth
    // is bounded by Barnard's star with a proper motion of ~10.5 arcsec/year.
    // When intersecting the bounding circle of a FOV with chunk boundaries
    // to determine which chunks must be loaded for association, the bounding
    // circle must be padded to ensure all relevant objects are loaded.
    try {
        for (int s = 0; s < numStripes; ++s) {
            ChunkVector &  vec = stripes[s];
            Size const numChunks = vec.size();

            // Loop over chunks in stripe
            for (Size c = 0; c < numChunks; ++c) {
                Chunk * const ch = &vec[c];
                int const numBlocks  = ch->blocks();
                int i = 0;

                // loop over blocks in chunk
                for (int b = 0; b < numBlocks; ++b) {
                    int const numEntries = ch->entries(b);
                    Data * const block = ch->getBlock(b);
                    ChunkEntryFlag const * const flags = ch->getFlagBlock(b);

                    // loop over entries in block
                    for (int e = 0; e < numEntries; ++e, ++i) {
                        if ((flags[e] & Chunk::DELETED) == 0) {
                            std::pair<double, double> pos = correctProperMotion(block[e], epoch);
                            index.insert(pos.first, pos.second, &block[e], ch, i);
                        }
                    }
                }
            }
        }
    } catch(...) {
       index.clear();
       throw LSST_EXCEPT(ex::RuntimeError, "Failed to build zone index"); 
    }

    watch.stop();
    Log log(Log::getDefaultLog(), "lsst.ap");
    int numElements = static_cast<int>(index.size());
    Rec(log, Log::INFO) << "inserted elements into zone index" <<
        Prop<int>("numElements", numElements) <<
        Prop<double>("time", watch.seconds()) << Rec::endr;

    // zone structure is filled, sort individual zones (on right ascension)
    watch.start();
    index.sort();
    watch.stop();
    Rec(log, Log::INFO) << "sorted zone index" <<
        Prop<int>("numElements", numElements) <<
        Prop<double>("time", watch.seconds()) << Rec::endr;
}

std::pair<double, double> lsst::ap::detail::correctProperMotion	(	Object const &	obj,
		double const	epoch
	)

Returns

the proper motion corrected position of the given object

Definition at line 149 of file Stages.cc.

                                                                                   {
    static double const RAD_PER_MAS = (RADIANS_PER_DEGREE/360000.0);
    // (rad/mas)*(sec/year)/(km/AU)
    static double const SCALE = RAD_PER_MAS*(365.25*86400/149597870.691);

    double ra   = radians(obj.getRa());   // rad
    double decl = radians(obj.getDec()); // rad

    // Convert ra, dec to unit vector in cartesian coordinate system
    double const sinDecl = std::sin(decl);
    double const cosDecl = std::cos(decl);
    double const sinRa   = std::sin(ra);
    double const cosRa   = std::cos(ra);
    double x = cosRa * cosDecl;
    double y = sinRa * cosDecl;
    double z = sinDecl;

    // compute space motion vector (radians per year)
    double const pmRa   = obj.getMuRa()*RAD_PER_MAS/cosDecl;
    double const pmDecl = obj.getMuDecl()*RAD_PER_MAS;
    // divide radial velociy by distance to source
    double const w = obj.getParallax()*obj.getRadialVelocity()*SCALE;

    double const mx = - pmRa*y - pmDecl*cosRa   + x*w;
    double const my =   pmRa*x - pmDecl*sinRa   + y*w;
    double const mz =            pmDecl*cosDecl + z*w;

    // Linear interpolation of position
    double const dt = (epoch - obj.getEpoch()) * (1/365.25); // julian years
    x += mx*dt;
    y += my*dt;
    z += mz*dt;

    // Store unit vector for corrected position, convert back to
    // spherical coords and store scaled integer ra/dec
    double d2 = x*x + y*y;
    ra   = (d2 == 0.0) ? 0 : degrees(std::atan2(y, x));
    decl = (z  == 0.0) ? 0 : degrees(std::atan2(z, std::sqrt(d2)));
    if (ra < 0.0) {
        ra += 360.0;
    }
    return std::make_pair(ra, decl);
}

template<typename ManagerT >

ManagerT* lsst::ap::detail::getSingleton	(	char const *const	shmObjName,
		char const *const	shmLockName
	)

Definition at line 129 of file ChunkManager.cc.

                                                                                       {

    static Mutex mutex;
    static ManagerT * singleton = 0;

    ScopedLock<Mutex> lck(mutex);
    if (singleton != 0) {
        return singleton;
    }

    // Block interference from other processes
    BootstrapLock blck(shmLockName);

    std::size_t const pageSize  = static_cast<std::size_t>(::getpagesize());
    std::size_t const numBytes  = (ManagerT::size() + pageSize - 1) & ~(pageSize - 1);

    // 1. try to create shared memory object
    int fd = ::shm_open(shmObjName, O_RDWR | O_CREAT | O_EXCL, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);

    if (fd == -1) {

        // failed ...
        if (errno != EEXIST) {
            throw LSST_EXCEPT(ex::RuntimeError,
                (boost::format("shm_open(): failed to create shared memory object %1%, errno: %2%")
                    % shmObjName % errno).str());
        }
        // because the shared memory object already exists -- so try to open existing object instead
        fd = ::shm_open(shmObjName, O_RDWR, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
        if (fd == -1) {
            throw LSST_EXCEPT(ex::RuntimeError,
                (boost::format("shm_open(): failed to open existing shared memory object %1%, errno: %2%")
                    % shmObjName % errno).str());
        }
        ScopeGuard g(boost::bind(::close, fd));

        void * mem = ::mmap(0, numBytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
        if (mem == 0) {
            throw LSST_EXCEPT(ex::RuntimeError,
                (boost::format("mmap(): failed to map %1% bytes of shared memory object %2%, errno: %3%")
                    % numBytes % shmObjName % errno).str());
        }
        singleton = static_cast<ManagerT *>(mem);
        g.dismiss();

    } else {

        // shared memory object was created (initially of zero size)
        ScopeGuard g1(boost::bind(::shm_unlink, shmObjName));
        ScopeGuard g2(boost::bind(::close, fd));

        // set size of shared memory object
        if (::ftruncate(fd, numBytes) != 0) {
            throw LSST_EXCEPT(ex::RuntimeError,
                (boost::format(
                    "ftruncate(): failed to set size of shared memory object %1% to %2% bytes, errno: %3%")
                    % shmObjName % numBytes % errno).str());
        }

        // map shared memory object
        void * mem = ::mmap(0, numBytes, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
        if (mem == 0) {
            throw LSST_EXCEPT(ex::RuntimeError,
                (boost::format("mmap(): failed to map %1% bytes of shared memory object %2%, errno: %3%")
                    % numBytes % shmObjName % errno).str());
        }
        ScopeGuard g3(boost::bind(::munmap, mem, numBytes));

        new (mem) ManagerT();
        singleton = static_cast<ManagerT *>(mem);
        g3.dismiss();
        g2.dismiss();
        g1.dismiss();
    }

    // Try to lock chunk pages into memory to avoid swapping, but ignore failure to do so
    // - Solaris      : process must be run as root
    // - Linux/Darwin : process must be run as root or RLIMIT_MEMLOCK should be set to a large value
    ::mlock(static_cast<void *>(singleton), numBytes);

    // Note: we rely on the system to munmap and close the file descriptor above at process exit.

    return singleton;
}

int lsst::ap::detail::hash ( boost::uint32_t  key )

inline

Returns the 32 bit hash of a 32 bit value using Thomas Wang's mixing function.

Definition at line 55 of file ChunkManagerImpl.cc.

                                   {
    key = (~key) + (key << 15); // key = (key << 15) - key - 1;
    key = key ^ (key >> 12);
    key = key + (key << 2);
    key = key ^ (key >> 4);
    key = key * 2057;           // key = (key + (key << 3)) + (key << 11);
    key = key ^ (key >> 16);
    return static_cast<int>(key);
}

int lsst::ap::detail::hash ( int  key )

inline

Definition at line 65 of file ChunkManagerImpl.cc.

                         {
    return hash(static_cast<boost::uint32_t>(key));
}

template<typename WordT >

WordT lsst::ap::detail::maskForBit ( int const  i )

inline

Definition at line 84 of file Bitset.h.

                                     {
    return static_cast<WordT>(1) << (i & (BitTraits<WordT>::BITS_PER_WORD - 1));
}

template<typename WordT >

void lsst::ap::detail::resetBits	(	WordT *const	words,
		int const *const	indexes,
		int const	numBitsToReset,
		int const	numBits
	)

Given the array words containing numBits bits and a list of numBitsToReset indexes, resets (to zero) the bit at each index.

Parameters

[in,out]	words	A memory buffer containing at least numBits bits
[in]	indexes	An array of at least numBitsToReset integers, each corresponding to the index of a single bit to reset (to zero).
[in]	numBitsToReset	The number of bits to reset (to zero)
[in]	numBits	The number of bits in words

Definition at line 212 of file Bitset.cc.

  {
    assert(words != 0 && numBits > 0 && "null or empty bitset");
    assert(indexes != 0 && numBitsToReset >= 0 && "null or empty index array");

    for (int i = 0; i < numBitsToReset; ++i) {
        int const j = indexes[i];
        assert(j >= 0 && j < numBits);
        words[wordForBit<WordT>(j)] &= ~ maskForBit<WordT>(j);
    }
}

template<typename WordT >

bool lsst::ap::detail::setBits	(	int *const	indexes,
		WordT *const	words,
		int const	numBitsToSet,
		int const	numBits
	)

Attempts to find at least numBitsToSet zero bits in the array words (containing at least numBits bits). If numBitsToSet zero bits cannot be found, false is returned. Otherwise, the first numBitsToSet zero bits in words are set to one and true is returned. The indexes of the bits which were set to one are stored in the indexes array.

Parameters

[out]	indexes	An array of at least numBitsToSet integers, which is filled with the indexes of zero bits that were set (to one).
[in,out]	words	A memory buffer containing at least numBits bits
[in]	numBitsToSet	The number of bits to reset (to zero)
[in]	numBits	The number of bits in words

Returns

true if numBitsToSet zero bits were found and set to one in words, false otherwise.

Definition at line 149 of file Bitset.cc.

  {
    assert(words   != 0 && numBits > 0 && "null or empty bitset");
    assert(indexes != 0 && numBitsToSet > 0 && "null or empty index array");

    bool  const special = (numBits & (BitTraits<WordT>::BITS_PER_WORD - 1)) > 0;
    WordT const last    = ~(maskForBit<WordT>(numBits) - 1);
    int   const nwords  = (numBits + (BitTraits<WordT>::BITS_PER_WORD - 1)) >>
                         BitTraits<WordT>::BITS_PER_WORD_LOG2;

    int zeroes = numBitsToSet;
    int i;
    for (i = 0; i < nwords - special && zeroes > 0; ++i) {
        zeroes -= BitTraits<WordT>::BITS_PER_WORD -
                  populationCount(static_cast<WordT>(words[i] & BitTraits<WordT>::WORD_MASK));
    }
    if (zeroes > 0 && special) {
        zeroes -= BitTraits<WordT>::BITS_PER_WORD -
                  populationCount(static_cast<WordT>((words[i] & BitTraits<WordT>::WORD_MASK) | last));
    }
    if (zeroes > 0) {
        // didn't find enough zeroes
        return false;
    }

    zeroes = 0;
    for (int i = 0; zeroes < numBitsToSet; ++i) {
        WordT w = words[i] & BitTraits<WordT>::WORD_MASK;
        if (w == BitTraits<WordT>::WORD_MASK) {
            continue;
        }
        WordT mask = 1;
        for (int j = 0; j < BitTraits<WordT>::BITS_PER_WORD; ++j, mask <<= 1) {
            if ((w & mask) == 0) {
                indexes[zeroes++] = (i << BitTraits<WordT>::BITS_PER_WORD_LOG2) + j;
                w |= mask;
                if (zeroes == numBitsToSet) {
                    break;
                }
            }
        }
        words[i] = w;
    }

    return true;
}

template<typename WordT >

int lsst::ap::detail::wordForBit ( int const  i )

inline

Definition at line 79 of file Bitset.h.

                                   {
    return (i >> BitTraits<WordT>::BITS_PER_WORD_LOG2);
}