#include <Histogram.h>

Public Member Functions
	Histogram (double binSize, double minValue, double maxValue)

	~Histogram ()

void	add (double value, const T &ref)

double	findPeak (double minVal, double maxVal) const

bool	hasSinglePeak (double minVal, double maxVal) const

double	findValley (double minVal, double maxVal) const

double	findFirstValueAfter (double start) const

double	findFirstValleyAfter (double val1, bool hasPoissonNoise=false) const

double	findFirstValleyBefore (double val1, bool hasPoissonNoise=false) const

double	findFirstPeakAfter (double val1, bool hasPoissonNoise=false) const

double	findFirstPeakBefore (double val1, bool hasPoissonNoise=false) const

int	getTotalCountBefore (double val1) const

int	getTotalCountAfter (double val1) const

int	operator[] (double value) const

double	findThresh (double minVal, double maxVal) const

std::vector< T >	getRefsInRange (double min, double max) const

std::vector< double >	getValuesInRange (double min, double max) const

int	getRefinedPeakCount (double *peak) const

int	getRefinedValleyCount (double *valley) const

void	print (std::ostream &fout, double val1=-1.e10, double val2=1.e10) const

Private Member Functions
int	getCount (int i) const

int	index (double value) const

double	value (int i) const

Private Attributes
double	_binSize

double	_minValue

double	_maxValue

std::vector< std::vector< T > >	_refs

std::vector< std::vector < double > >	_values

Detailed Description

template<typename T>
class lsst::meas::algorithms::shapelet::Histogram< T >

Definition at line 13 of file Histogram.h.

Constructor & Destructor Documentation

template<typename T >

lsst::meas::algorithms::shapelet::Histogram< T >::Histogram	(	double	binSize,
		double	minValue,
		double	maxValue
	)

Definition at line 38 of file Histogram.cc.

                                                                             :
         _binSize(binSize)
     {
         // Add an extra bin to prevent rounding errors from giving you a
         // value outsize the allowed range.
         int nBins = int(ceil((maxValue-minValue)/_binSize))+1;
         _minValue = minValue - binSize/2.;
         _maxValue = maxValue + binSize/2.;
         _refs = std::vector<std::vector<T> >(nBins);
         _values = std::vector<std::vector<double> >(nBins);
         if (value(nBins-1) > _maxValue) _maxValue = value(nBins-1)+binSize/4.;
 
         dbg<<"made histogram:\n";
         dbg<<"minvalue="<<_minValue<<
             " index(minvalue)="<<index(_minValue)<<std::endl;
         dbg<<"maxvalue="<<_maxValue<<
             " index(maxvalue)="<<index(_maxValue)<<std::endl;
         dbg<<"mini=0  value(mini)="<<value(0)<<std::endl;
         dbg<<"maxi="<<_refs.size()-1;
         dbg<<" value(maxi)="<<value(_refs.size()-1)<<std::endl;
     }

template<typename T>

lsst::meas::algorithms::shapelet::Histogram< T >::~Histogram ( )

inline

Definition at line 18 of file Histogram.h.

18 {}

Member Function Documentation

template<typename T >

void lsst::meas::algorithms::shapelet::Histogram< T >::add	(	double	value,
		const T &	ref
	)

Definition at line 61 of file Histogram.cc.

     {
         _refs[index(value)].push_back(ref);
         _values[index(value)].push_back(value);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findFirstPeakAfter	(	double	val1,
		bool	hasPoissonNoise = `false`
	)		const

Definition at line 204 of file Histogram.cc.

     {
         if (val1<_minValue) val1 = _minValue;
         int i1=index(val1);
         const int nRefs = _refs.size();
         Assert(i1 < nRefs);
         // Keep going until it stops rising
         int iPeak=i1;
         double maxCount = getCount(i1);
         int iCheck=i1+1;
         while(iCheck < nRefs && maxCount <= getCount(iCheck) +
               (hasPoissonNoise ? int(sqrt(float(getCount(iCheck)))) : 0)) {
             if (getCount(iCheck) > maxCount) {
                 maxCount = getCount(iCheck);
                 iPeak = iCheck;
             }
             if (++iCheck == nRefs) break;
         }
         return value(iPeak);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findFirstPeakBefore	(	double	val1,
		bool	hasPoissonNoise = `false`
	)		const

Definition at line 226 of file Histogram.cc.

     {
         if (val1>_maxValue) val1 = _maxValue;
         int i1=index(val1);
         Assert(i1 < int(_refs.size()));
         // Keep going until it stops rising
         int iPeak=i1;
         double maxCount = getCount(i1);
         if (i1 > 0) {
             int iCheck=i1-1;
             while(maxCount <= getCount(iCheck) +
                   (hasPoissonNoise ? int(sqrt(float(getCount(iCheck)))) : 0)) {
                 if (getCount(iCheck) > maxCount) {
                     maxCount = getCount(iCheck);
                     iPeak = iCheck;
                 }
                 if (iCheck == 0) break;
                 --iCheck;
             }
         }
         return value(iPeak);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findFirstValleyAfter	(	double	val1,
		bool	hasPoissonNoise = `false`
	)		const

Definition at line 140 of file Histogram.cc.

     {
         dbg<<"Start FindFirstValleyAfter "<<val1<<std::endl;
         if (val1<_minValue) val1 = _minValue;
         int i1=index(val1);
         Assert(i1<int(_refs.size()));
         // Keep going until it stops dropping
         // The sqrt bit allows for Poisson noise in the counts
         // Otherwise you can get caught early by ledges.
         // eg. 15 9 10 4 1 3 3 ... would stop at the 9 instead of at the 1
         int iValley=i1;
         double minCount = getCount(i1);
         int iCheck=i1+1;
         const int nRefs = _refs.size();
         while(iCheck < nRefs && getCount(iCheck) <= minCount +
               (hasPoissonNoise ? int(sqrt(float(getCount(iCheck)))) : 0)) {
             if (getCount(iCheck) < minCount) {
                 minCount = getCount(iCheck);
                 iValley = iCheck;
             }
             if (++iCheck == nRefs) break;
         }
         dbg<<"valley = "<<value(iValley)<<std::endl;
         return value(iValley);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findFirstValleyBefore	(	double	val1,
		bool	hasPoissonNoise = `false`
	)		const

Definition at line 167 of file Histogram.cc.

     {
         if (val1>_maxValue) val1 = _maxValue;
         int i1=index(val1);
         Assert(i1<int(_refs.size()));
         // Keep going until it stops dropping
         // The sqrt bit allows for Poisson noise in the counts
         // Otherwise you can get caught early by ledges.
         // eg. 15 9 10 4 1 3 3 ... would stop at the 9 instead of at the 1
         int iValley=i1;
         double minCount = getCount(i1);
         if (i1 > 0) {
             int iCheck=i1-1;
             while(getCount(iCheck) <= minCount +
                   (hasPoissonNoise ? int(sqrt(float(getCount(iCheck)))) : 0)) {
                 if (getCount(iCheck) < minCount) {
                     minCount = getCount(iCheck);
                     iValley = iCheck;
                 }
                 if (iCheck == 0) break;
                 --iCheck;
             }
         }
         return value(iValley);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findFirstValueAfter ( double start ) const

Definition at line 194 of file Histogram.cc.

     {
         int i = (start<_minValue) ? index(_minValue) : index(start)+1;
         const int nValues = _values.size();
         while (i<nValues && _values[i].size() == 0) ++i;
         if (i == nValues) return _maxValue;
         else return *std::min_element(_values[i].begin(),_values[i].end());
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findPeak	(	double	minVal,
		double	maxVal
	)		const

Definition at line 76 of file Histogram.cc.

     {
         if (val1<_minValue) val1 = _minValue;
         if (val2>_maxValue) val2 = _maxValue;
         int i1=index(val1), i2=index(val2);
         Assert(i1<int(_refs.size()));
         Assert(i2<int(_refs.size()));
         int iPeak=i1;
         int maxCount = getCount(i1);
 
         for(int i=i1+1;i<=i2;++i) {
             if (getCount(i) > maxCount) {
                 maxCount = getCount(i);
                 iPeak = i;
             }
         }
         return value(iPeak);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findThresh	(	double	minVal,
		double	maxVal
	)		const

Definition at line 347 of file Histogram.cc.

     {
         const double EPSILON = 1.e-6;
 
         dbg<<"starting FindThresh\n";
         if (dbgout) print(*dbgout);
 
         double nTot=0.,meanTot=0.,meanSqTot=0.;
 
         if (minVal < _minValue) minVal = _minValue;
         if (maxVal > _maxValue) maxVal = _maxValue;
         int i1 = index(minVal);
         int i2 = index(maxVal);
         Assert(i1 < int(_refs.size()));
         Assert(i2 < int(_refs.size()));
 
         // First calculate the "tot" values:
         for(int i=i1;i<=i2;++i) {
             double dcount = getCount(i);
             nTot += dcount;
             meanTot += dcount*i;
             meanSqTot += dcount*i*i;
         }
         meanTot /= nTot;
         meanSqTot /= nTot;
         dbg<<"ntot="<<nTot<<", meantot="<<meanTot<<
             ", meansqtot="<<meanSqTot<<std::endl;
 
         double sumN=0.,sumNi=0.,sumNii=0.,varBest=-1.;
         int iBest=-1;
 
         const double minVariance = 0.25;
 
         for(int i=i1;i<=i2;++i) {
             double dcount = getCount(i);
             sumN += dcount;
             sumNi += dcount*i;
             sumNii += dcount*i*i;
             double fL = sumN/nTot;
             double fR = 1.-fL;
             Assert(fL >=0. && fR >=0.);
             if (fL == 0. || fR == 0.) continue; // var = 0
             double xL = sumNi/nTot;  // actuall fL*xL now
             double xR = (meanTot-xL)/fR;
             xL /= fL;
             double vL = sumNii/nTot; // actually fL*(vL+xL^2)
             double vR = (meanSqTot-vL)/fR - xR*xR;
             vL = vL/fL - xL*xL;
             if (vL<minVariance) vL = minVariance;
             if (vR<minVariance) vR = minVariance;
             double var = fL*fR*std::pow(xR-xL,2)/sqrt(vL*vR);
             Assert(var >= 0.);
             if (var > varBest+EPSILON) {
                 iBest = i;
                 varBest = var;
             }
         }
 
         dbg<<"besti = "<<iBest<<", bestvar = "<<varBest<<std::endl;
         Assert(iBest >= 0);
         Assert(varBest > 0.);
         // +1 below forces thresh to be at least 1 bin away from first peak.
         // otherwise if peak is only 1 bin, can get iBest = that bin.
         dbg<<"returning thresh = "<<value(iBest)<<std::endl;
         return value(iBest);
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::findValley	(	double	minVal,
		double	maxVal
	)		const

Definition at line 113 of file Histogram.cc.

     {
         if (val1<_minValue) val1 = _minValue;
         if (val2>_maxValue) val2 = _maxValue;
         int i1=index(val1), i2=index(val2);
         Assert(i1<int(_refs.size()));
         Assert(i2<int(_refs.size()));
         int iValley=i1;
         int minCount = getCount(i1);
 
         for(int i=i1+1;i<=i2;++i) {
             if (getCount(i) < minCount) {
                 minCount = getCount(i);
                 iValley = i;
             }
         }
         // Make sure you're not still sloping down.  If so, keep going.
         if (iValley==i2) {
             int nRefs = _refs.size();
             while(iValley+1 < nRefs && getCount(iValley+1) < getCount(iValley)) 
                 ++iValley;
         }
 
         return value(iValley);
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::getCount ( int i ) const

private

Definition at line 68 of file Histogram.cc.

     {
         Assert(i<int(_refs.size()));
         Assert(_values[i].size() == _refs[i].size());
         return _refs[i].size();
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::getRefinedPeakCount ( double * peak ) const

Definition at line 271 of file Histogram.cc.

     {
         std::vector<double> vals = 
             getValuesInRange(*peak-1.5*_binSize,
                              *peak+1.5*_binSize);
         Assert(vals.size()>=1);
         std::sort(vals.begin(),vals.end());
         int bestCount=0;
         double bestVal=*peak;
         // i is the object at the left side of the virtual bin
         // j is the first object after the end of the virtual bin
         const int nVals = vals.size();
         for(int i=0,j=0;j<nVals;++i) {
             // Determine j for this value of i
             while (j<nVals && vals[j] < vals[i]+_binSize) ++j;
             // If the corresponding count is better than bestCount, update it
             Assert(j>i);
             if (j-i > bestCount) {
                 bestCount = j-i;
                 bestVal = (vals[i]+vals[j-1])/2.;
             }
         }
         *peak = bestVal;
         return bestCount;
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::getRefinedValleyCount ( double * valley ) const

Definition at line 302 of file Histogram.cc.

     {
         std::vector<double> vals = 
             getValuesInRange(*valley-1.5*_binSize,
                              *valley+1.5*_binSize);
         const int nVals = vals.size();
         if (nVals == 0) return 0;
         std::sort(vals.begin(),vals.end());
         // start high since want to find the lowest count
         int bestCount=nVals;
         double bestVal=*valley;
         // i is the last object before the start of the virtual bin
         // j is the first object after the end of the virtual bin
         for(int i=0,j=0; i<nVals && vals[i]<*valley+0.5*_binSize; ++i) {
             // Determine j for this value of i
             while (j<nVals && vals[j] < vals[i]+_binSize) ++j;
             // If we can get rid of any i values, given this j, do so.
             double nextBinStart = j<nVals ? vals[j] : *valley+1.5*_binSize;
             while (i+1<nVals && vals[i+1]+_binSize < nextBinStart) ++i;
             // If the corresponding count is better than bestCount, update it
             Assert(j>i);
             if (j-i-1 < bestCount) {
                 bestCount = j-i-1;
                 bestVal = vals[i]+_binSize/2.;
             }
         }
         *valley = bestVal;
         return bestCount;
     }

template<typename T >

std::vector< T > lsst::meas::algorithms::shapelet::Histogram< T >::getRefsInRange	(	double	min,
		double	max
	)		const

Definition at line 448 of file Histogram.cc.

     {
         if (min < _minValue) min = _minValue;
         if (max > _maxValue) max = _maxValue;
         int i1=index(min);
         int i2=index(max);
 
         dbg<<"in getrefs: min,max = "<<min<<','<<max<<std::endl;
 
         std::vector<T> temp;
         const int nRefs1 = _refs[i1].size();
         for(int k=0; k<nRefs1; ++k) {
             if(_values[i1][k]>=min) {
                 dbg<<"i1 - add ref @ "<<_values[i1][k]<<std::endl;
                 temp.push_back(_refs[i1][k]);
             }
         }
         for(int i=i1+1;i<i2;++i) {
             if (_refs[i].size() > 0) 
                 dbg<<"i - add all ("<<_refs[i].size()<<") refs near "<<
                     _values[i].front()<<std::endl;
             temp.insert(temp.end(),_refs[i].begin(),_refs[i].end());
         }
         const int nRefs2 = _refs[i2].size();
         for(int k=0; k<nRefs2; ++k) {
             if(_values[i2][k] <= max) {
                 dbg<<"i2 - add ref @ "<<_values[i2][k]<<std::endl;
                 temp.push_back(_refs[i2][k]);
             }
         }
         return temp;
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::getTotalCountAfter ( double val1 ) const

Definition at line 260 of file Histogram.cc.

     {
         int count=0;
         int i1 = index(val1);
         const int nRefs = _refs.size();
         Assert(i1 < nRefs);
         for(int i=i1+1;i<nRefs;++i) count += getCount(i);
         return count;
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::getTotalCountBefore ( double val1 ) const

Definition at line 250 of file Histogram.cc.

     {
         int count=0;
         int i1 = index(val1);
         Assert(i1 < int(_refs.size()));
         for(int i=0;i<i1;++i) count += getCount(i);
         return count;
     }

template<typename T >

std::vector< double > lsst::meas::algorithms::shapelet::Histogram< T >::getValuesInRange	(	double	min,
		double	max
	)		const

Definition at line 482 of file Histogram.cc.

     {
         if (min < _minValue) min = _minValue;
         if (max > _maxValue) max = _maxValue;
         int i1=index(min);
         int i2=index(max);
 
         std::vector<double> temp;
         const int nVals1 = _values[i1].size();
         for(int k=0;k<nVals1;++k)
             if(_values[i1][k]>=min) temp.push_back(_values[i1][k]);
         for(int i=i1+1;i<i2;++i)
             temp.insert(temp.end(),_values[i].begin(),_values[i].end());
         const int nVals2 = _values[i2].size();
         for(int k=0;k<nVals2;++k)
             if(_values[i2][k]<=max) temp.push_back(_values[i2][k]);
         return temp;
     }

template<typename T >

bool lsst::meas::algorithms::shapelet::Histogram< T >::hasSinglePeak	(	double	minVal,
		double	maxVal
	)		const

Definition at line 96 of file Histogram.cc.

     {
         if (val1<_minValue) val1 = _minValue;
         if (val2>_maxValue) val2 = _maxValue;
         int i1=index(val1), i2=index(val2);
         Assert(i1<int(_refs.size()));
         Assert(i2<int(_refs.size()));
         int iPeak1=i1,iPeak2=i2;
         while(iPeak1+1<=i2 && getCount(iPeak1+1)>=getCount(iPeak1)) ++iPeak1;
         while(iPeak2-1>=i1 && getCount(iPeak2-1)>=getCount(iPeak2)) --iPeak2;
         return iPeak1 >= iPeak2; // Usually ==, but if plateau, they swap
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::index ( double value ) const

private

Definition at line 503 of file Histogram.cc.

     {
         Assert (value >= _minValue && value <= _maxValue);
         return int(floor((value - _minValue)/_binSize));
     }

template<typename T >

int lsst::meas::algorithms::shapelet::Histogram< T >::operator[] ( double value ) const

Definition at line 339 of file Histogram.cc.

     {
         if (index(value) == int(_refs.size())) return 0;
         Assert(index(value) < int(_refs.size()));
         return getCount(index(value));
     }

template<typename T >

void lsst::meas::algorithms::shapelet::Histogram< T >::print	(	std::ostream &	fout,
		double	val1 = `-1.e10`,
		double	val2 = `1.e10`
	)		const

Definition at line 517 of file Histogram.cc.

     {
         if (val1 < _minValue) val1 = _minValue;
         if (val2 > _maxValue) val2 = _maxValue;
         int i1=index(val1);
         int i2=index(val2);
         Assert(i1 < int(_refs.size()));
         Assert(i2 < int(_refs.size()));
         Form sci2; sci2.sci().width(10).prec(2);
 
         for(int i=i1;i<=i2;++i) {
             //if ((i-i1)%5 == 0) 
             fout << sci2(value(i));
             //else fout << "           ";
             for(int j=0;j<getCount(i);++j) fout<<'*';
             fout<<std::endl;
         }
     }

template<typename T >

double lsst::meas::algorithms::shapelet::Histogram< T >::value ( int i ) const

private

Definition at line 510 of file Histogram.cc.

     {
         Assert(i<int(_refs.size()));
         return _binSize*(i+0.5)+_minValue;
     }

Member Data Documentation

template<typename T>

double lsst::meas::algorithms::shapelet::Histogram< T >::_binSize

private

Definition at line 45 of file Histogram.h.

template<typename T>

double lsst::meas::algorithms::shapelet::Histogram< T >::_maxValue

private

Definition at line 45 of file Histogram.h.

template<typename T>

double lsst::meas::algorithms::shapelet::Histogram< T >::_minValue

private

Definition at line 45 of file Histogram.h.

template<typename T>

std::vector<std::vector<T> > lsst::meas::algorithms::shapelet::Histogram< T >::_refs

private

Definition at line 46 of file Histogram.h.

template<typename T>

std::vector<std::vector<double> > lsst::meas::algorithms::shapelet::Histogram< T >::_values

private

Definition at line 47 of file Histogram.h.

The documentation for this class was generated from the following files:

/home/lsstsw/stack/Linux64/meas_algorithms/11.0-2-gb8b8ce7/include/lsst/meas/algorithms/shapelet/Histogram.h
/home/lsstsw/stack/Linux64/meas_algorithms/11.0-2-gb8b8ce7/src/shapelet/Histogram.cc

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

template<typename T> class lsst::meas::algorithms::shapelet::Histogram< T >

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<typename T>
class lsst::meas::algorithms::shapelet::Histogram< T >