Index: /issm/trunk-jpl/src/c/classes/kriging/Covertree.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/kriging/Covertree.cpp (revision 19152) +++ /issm/trunk-jpl/src/c/classes/kriging/Covertree.cpp (revision 19153) @@ -1,10 +1,8 @@ #include "../classes.h" -//#include #include -//#include #include /*Constructors/Destructors*/ -Covertree::Covertree(int maxLevel,const std::vector& points){ +Covertree::Covertree(int maxLevel,const std::vector& points){/*{{{*/ this->base = 2.; _root=NULL; @@ -16,7 +14,6 @@ this->insert(*it);//adds data to the covertree object } -} - -Covertree::~Covertree(){ +}/*}}}*/ +Covertree::~Covertree(){/*{{{*/ if(_root==NULL) return; //Get all of the root's children (from any level), @@ -31,8 +28,8 @@ delete byeNode; } -} +}/*}}}*/ -/*Methods*/ -std::vector Covertree::kNearestNodes(const Observation& p, const unsigned int& k) const{ + /*Methods*/ +std::vector Covertree::kNearestNodes(const Observation& p, const unsigned int& k) const{/*{{{*/ if(_root==NULL) return std::vector(); //maxDist is the kth nearest known point to p, and also the farthest @@ -80,6 +77,6 @@ } return kNN; -} -bool Covertree::insert_rec(const Observation& p, const std::vector& Qi, const int& level){ +}/*}}}*/ +bool Covertree::insert_rec(const Observation& p, const std::vector& Qi, const int& level){/*{{{*/ std::vector > Qj; double sep = pow(base,level); @@ -119,7 +116,6 @@ } } -} - -void Covertree::remove_rec(const Observation& p, std::map >& coverSets, int level, bool& multi){ +}/*}}}*/ +void Covertree::remove_rec(const Observation& p, std::map >& coverSets, int level, bool& multi){/*{{{*/ std::vector& Qi = coverSets[level]; std::vector& Qj = coverSets[level-1]; @@ -222,14 +218,9 @@ } } -} - -int Covertree::get_numberofobs(){ +}/*}}}*/ +int Covertree::get_numberofobs(){/*{{{*/ return _numNodes; -} - -std::pair -Covertree::distance(const Observation& p, - const std::vector& Q) -{ +}/*}}}*/ +std::pair Covertree::distance(const Observation& p, const std::vector& Q){/*{{{*/ double minDist = 1.e+50; CoverTreeNode* minNode; @@ -243,8 +234,6 @@ } return std::make_pair(minDist,minNode); -} - -void Covertree::insert(const Observation& newObservation) -{ +}/*}}}*/ +void Covertree::insert(const Observation& newObservation){/*{{{*/ if(_root==NULL) { _root = new CoverTreeNode(newObservation); @@ -265,8 +254,6 @@ _maxLevel); } -} - -void Covertree::remove(const Observation& p) -{ +}/*}}}*/ +void Covertree::remove(const Observation& p){/*{{{*/ //Most of this function's code is for the special case of removing the root if(_root==NULL) return; @@ -305,7 +292,6 @@ _root=newRoot; } -} - -std::vector Covertree::kNearestNeighbors(const Observation& p, const unsigned int& k) const{ +}/*}}}*/ +std::vector Covertree::kNearestNeighbors(const Observation& p, const unsigned int& k) const{/*{{{*/ if(_root==NULL) return std::vector(); std::vector v = kNearestNodes(p, k); @@ -318,7 +304,6 @@ } return kNN; -} - -void Covertree::print() const{ +}/*}}}*/ +void Covertree::print() const{/*{{{*/ int d = _maxLevel-_minLevel+1; std::vector Q; @@ -346,18 +331,13 @@ std::cout << "\n\n"; } -} +}/*}}}*/ -Covertree::CoverTreeNode* Covertree::getRoot() const -{ +Covertree::CoverTreeNode* Covertree::getRoot() const{/*{{{*/ return _root; -} - -Covertree::CoverTreeNode::CoverTreeNode(const Observation& p) { +}/*}}}*/ +Covertree::CoverTreeNode::CoverTreeNode(const Observation& p) {/*{{{*/ _observations.push_back(p); -} - -std::vector -Covertree::CoverTreeNode::getChildren(int level) const -{ +}/*}}}*/ +std::vector Covertree::CoverTreeNode::getChildren(int level) const{/*{{{*/ std::map >::const_iterator it = _childMap.find(level); @@ -366,13 +346,9 @@ } return std::vector(); -} - -void Covertree::CoverTreeNode::addChild(int level, CoverTreeNode* p) -{ +}/*}}}*/ +void Covertree::CoverTreeNode::addChild(int level, CoverTreeNode* p){/*{{{*/ _childMap[level].push_back(p); -} - -void Covertree::CoverTreeNode::removeChild(int level, CoverTreeNode* p) -{ +}/*}}}*/ +void Covertree::CoverTreeNode::removeChild(int level, CoverTreeNode* p){/*{{{*/ std::vector& v = _childMap[level]; for(unsigned int i=0;i::iterator it = find(_observations.begin(), _observations.end(), p); if(it != _observations.end()) _observations.erase(it); -} - -double Covertree::CoverTreeNode::distance(const CoverTreeNode& p) const -{ +}/*}}}*/ +double Covertree::CoverTreeNode::distance(const CoverTreeNode& p) const{/*{{{*/ return _observations[0].distance(p.getObservation()); -} - -bool Covertree::CoverTreeNode::isSingle() const -{ +}/*}}}*/ +bool Covertree::CoverTreeNode::isSingle() const{/*{{{*/ return _observations.size() == 1; -} - -bool Covertree::CoverTreeNode::hasObservation(const Observation& p) const -{ +}/*}}}*/ +bool Covertree::CoverTreeNode::hasObservation(const Observation& p) const{/*{{{*/ return find(_observations.begin(), _observations.end(), p) != _observations.end(); -} - -const Observation& Covertree::CoverTreeNode::getObservation() const { return _observations[0]; } - -std::vector -Covertree::CoverTreeNode::getAllChildren() const -{ +}/*}}}*/ +const Observation& Covertree::CoverTreeNode::getObservation() const{/*{{{*/ + return _observations[0]; +}/*}}}*/ +std::vector Covertree::CoverTreeNode::getAllChildren() const{/*{{{*/ std::vector children; std::map >::const_iterator it; @@ -425,7 +389,6 @@ } return children; -} - -bool Covertree::isValidTree() const { +}/*}}}*/ +bool Covertree::isValidTree() const {/*{{{*/ if(_numNodes==0) return _root==NULL; @@ -465,3 +428,3 @@ } return true; -} +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/kriging/Covertree.h =================================================================== --- /issm/trunk-jpl/src/c/classes/kriging/Covertree.h (revision 19152) +++ /issm/trunk-jpl/src/c/classes/kriging/Covertree.h (revision 19153) @@ -8,8 +8,6 @@ class Covertree{ - /** - * Cover tree node. Consists of arbitrarily many points P, as long as - * they have distance 0 to each other. Keeps track of its children. - */ + /* Cover tree node. Consists of arbitrarily many points P, as long as they + * have distance 0 to each other. Keeps track of its children. */ class CoverTreeNode{ private: @@ -34,8 +32,6 @@ const std::vector& getObservations() { return _observations; } double distance(const CoverTreeNode& p) const; - - bool isSingle() const; - bool hasObservation(const Observation& o) const; - + bool isSingle() const; + bool hasObservation(const Observation& o) const; const Observation& getObservation() const; @@ -49,31 +45,22 @@ typedef std::pair distNodePair; - CoverTreeNode* _root; - unsigned int _numNodes; - int _maxLevel;//base^_maxLevel should be the max distance + CoverTreeNode *_root; + unsigned int _numNodes; + int _maxLevel; //base^_maxLevel should be the max distance //between any 2 points - int _minLevel;//A level beneath which there are no more new nodes. + int _minLevel; //A level beneath which there are no more new nodes. - std::vector - kNearestNodes(const Observation& o, const unsigned int& k) const; + std::vector kNearestNodes(const Observation& o, const unsigned int& k) const; /** * Recursive implementation of the insert algorithm (see paper). */ - bool insert_rec(const Observation& p, - const std::vector& Qi, - const int& level); + bool insert_rec(const Observation& p, const std::vector& Qi,const int& level); - /** - * Finds the node in Q with the minimum distance to p. Returns a - * pair consisting of this node and the distance. - */ - distNodePair distance(const Observation& p, - const std::vector& Q); + /* Finds the node in Q with the minimum distance to p. Returns a pair + * consisting of this node and the distance. */ + distNodePair distance(const Observation& p,const std::vector& Q); - void remove_rec(const Observation& p, - std::map >& coverSets, - int level, - bool& multi); + void remove_rec(const Observation& p, std::map >& coverSets, int level, bool& multi); public: @@ -89,6 +76,5 @@ */ - Covertree(int maxDist, - const std::vector& points=std::vector()); + Covertree(int maxDist,const std::vector& points=std::vector()); ~Covertree(); Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 19152) +++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 19153) @@ -59,7 +59,4 @@ _error_("Tree type "<treetype<<" not supported yet (1: quadtree, 2: covertree)"); } - - - } /*}}}*/ @@ -131,5 +128,5 @@ if(observations_list[i]quadtree->ClosestObs(&index,x[i],y[i]); if(index>=0){ @@ -157,10 +154,36 @@ void Observations::InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/ - int maxdepth = 20; + /*Intermediaries*/ + int maxdepth = 20; + IssmPDouble minspacing,mintrimming,maxtrimming; + + /*Checks*/ + _assert_(n); + + /*Get trimming limits*/ + options->Get(&mintrimming,"mintrimming",-1.e+21); + options->Get(&maxtrimming,"maxtrimming",+1.e+21); + options->Get(&minspacing,"minspacing",0.01); + if(minspacing<=0) _error_("minspacing must > 0"); + _printf0_("Generating covertree with a maximum depth " << maxdepth <<"... "); this->covertree=new Covertree(maxdepth); for(int i=0;icovertree->insert(Observation(x[i],y[i],observations_list[i])); + + /*First check limits*/ + if(observations_list[i]>maxtrimming) continue; + if(observations_list[i]0 && this->covertree->getRoot()){ + /*Get closest obs and see if it is too close*/ + std::vector kNN=(this->covertree->kNearestNeighbors(newobs,1)); + Observation oldobs = (*kNN.begin()); + if(oldobs.distance(newobs)covertree->insert(newobs); } _printf0_("done\n"); @@ -170,4 +193,18 @@ /*Methods*/ void Observations::ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/ + + switch(this->treetype){ + case 1: + this->ClosestObservationQuadtree(px,py,pobs,x_interp,y_interp,radius); + break; + case 2: + this->ClosestObservationCovertree(px,py,pobs,x_interp,y_interp,radius); + break; + default: + _error_("Tree type "<treetype<<" not supported yet (1: quadtree, 2: covertree)"); + } + +}/*}}}*/ +void Observations::ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/ /*Output and Intermediaries*/ @@ -220,4 +257,26 @@ }/*}}}*/ +void Observations::ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/ + + IssmPDouble hmin = UNDEF; + + if(this->covertree->getRoot()){ + /*Get closest obs and see if it is too close*/ + Observation newobs = Observation(x_interp,y_interp,0.); + std::vector kNN=(this->covertree->kNearestNeighbors(newobs,1)); + Observation observation = (*kNN.begin()); + hmin = observation.distance(newobs); + if(hmin<=radius){ + *px = observation.x; + *py = observation.y; + *pobs = observation.value; + return; + } + } + + *px = UNDEF; + *py = UNDEF; + *pobs = UNDEF; +}/*}}}*/ void Observations::Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius){/*{{{*/ @@ -234,10 +293,45 @@ } }/*}}}*/ +void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){/*{{{*/ + + /*Output and Intermediaries*/ + int nobs; + IssmPDouble *x = NULL; + IssmPDouble *y = NULL; + IssmPDouble *obs = NULL; + Observation *observation = NULL; + + nobs = this->Size(); + + if(nobs){ + x = xNew(nobs); + y = xNew(nobs); + obs = xNew(nobs); + for(int i=0;iSize();i++){ + observation=xDynamicCast(this->GetObjectByOffset(i)); + observation->WriteXYObs(&x[i],&y[i],&obs[i]); + } + } + + /*Assign output pointer*/ + *px=x; + *py=y; + *pobs=obs; + *pnobs=nobs; +}/*}}}*/ void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/ - if(this->treetype==2){ - this->ObservationList2(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata); - return; - } + switch(this->treetype){ + case 1: + this->ObservationListQuadtree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata); + break; + case 2: + this->ObservationListCovertree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata); + break; + default: + _error_("Tree type "<treetype<<" not supported yet (1: quadtree, 2: covertree)"); + } +}/*}}}*/ +void Observations::ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/ /*Output and Intermediaries*/ @@ -321,30 +415,46 @@ *pnobs=nobs; }/*}}}*/ -void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){/*{{{*/ - - /*Output and Intermediaries*/ - int nobs; - IssmPDouble *x = NULL; - IssmPDouble *y = NULL; - IssmPDouble *obs = NULL; - Observation *observation = NULL; - - nobs = this->Size(); - - if(nobs){ - x = xNew(nobs); - y = xNew(nobs); - obs = xNew(nobs); - for(int i=0;iSize();i++){ - observation=xDynamicCast(this->GetObjectByOffset(i)); - observation->WriteXYObs(&x[i],&y[i],&obs[i]); - } - } - - /*Assign output pointer*/ +void Observations::ObservationListCovertree(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){/*{{{*/ + + double *x = NULL; + double *y = NULL; + double *obs = NULL; + Observation observation=Observation(x_interp,y_interp,0.); + std::vector kNN; + + kNN=(this->covertree->kNearestNeighbors(observation, maxdata)); + //cout << "kNN's size: " << kNN.size() << " (maxdata = " <::iterator it; + if(radius>0.){ + for (it = kNN.begin(); it != kNN.end(); ++it) { + //(*it).print(); + //cout << "\n" << (*it).distance(observation) << endl; + if ((*it).distance(observation) > radius) { + break; + } + } + kNN.erase(it, kNN.end()); + } + + /*Allocate vectors*/ + x = new double[kNN.size()]; + y = new double[kNN.size()]; + obs = new double[kNN.size()]; + + /*Loop over all observations and fill in x, y and obs*/ + int i = 0; + for(it = kNN.begin(); it != kNN.end(); ++it) { + (*it).WriteXYObs((*it), &x[i], &y[i], &obs[i]); + i++; + } + *px=x; *py=y; *pobs=obs; - *pnobs=nobs; + *pnobs = kNN.size(); }/*}}}*/ void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){/*{{{*/ @@ -362,7 +472,4 @@ _assert_(pprediction); _assert_(power>0); - - /*If radius is not provided or is 0, return all observations*/ - if(radius==0) radius=this->quadtree->root->length; /*Get list of observations for current point*/ @@ -410,7 +517,4 @@ _assert_(pprediction && perror); - /*If radius is not provided or is 0, return all observations*/ - if(radius==0) radius=this->quadtree->root->length; - /*Get list of observations for current point*/ this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata); @@ -470,5 +574,4 @@ xDelete(B); xDelete(Lambda); - }/*}}}*/ void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/ @@ -617,46 +720,2 @@ }/*}}}*/ -void Observations::ObservationList2(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){/*{{{*/ - - - double *x = NULL; - double *y = NULL; - double *obs = NULL; - Observation observation=Observation(x_interp,y_interp,0.); - std::vector kNN; - - kNN=(this->covertree->kNearestNeighbors(observation, maxdata)); - //cout << "kNN's size: " << kNN.size() << " (maxdata = " <::iterator it; - if(radius>0.){ - for (it = kNN.begin(); it != kNN.end(); ++it) { - //(*it).print(); - //cout << "\n" << (*it).distance(observation) << endl; - if ((*it).distance(observation) > radius) { - break; - } - } - kNN.erase(it, kNN.end()); - } - - /*Allocate vectors*/ - x = new double[kNN.size()]; - y = new double[kNN.size()]; - obs = new double[kNN.size()]; - - /*Loop over all observations and fill in x, y and obs*/ - int i = 0; - for(it = kNN.begin(); it != kNN.end(); ++it) { - (*it).WriteXYObs((*it), &x[i], &y[i], &obs[i]); - i++; - } - - *px=x; - *py=y; - *pobs=obs; - *pnobs = kNN.size(); -}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.h =================================================================== --- /issm/trunk-jpl/src/c/classes/kriging/Observations.h (revision 19152) +++ /issm/trunk-jpl/src/c/classes/kriging/Observations.h (revision 19153) @@ -34,4 +34,6 @@ /*Methods*/ void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius); + void ClosestObservationQuadtree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius); + void ClosestObservationCovertree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius); void Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius); void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power); @@ -41,5 +43,6 @@ void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs); void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata); - void ObservationList2(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata); + void ObservationListQuadtree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata); + void ObservationListCovertree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata); void QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n); void Variomap(IssmDouble* gamma,IssmDouble *x,int n);