From 3fca927dbb900be9d9e3bb0884080e89f22a6f60 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Sat, 7 Sep 2019 10:27:26 -0400 Subject: [PATCH 01/31] Added a way to stride through spatial index values at a level. Improved the increment/decrement striding. Requires a function call. Would be nice to add something that could be inlined or a macro. --- VERSION | 2 +- include/STARE.h.in | 8 +++ src/EmbeddedLevelNameEncoding.C | 68 ++++++++++++++++++--- src/STARE.C | 29 +++++++++ tests/CUTE/STARE_test.cpp | 102 +++++++++++++++++++++++++++++++- tests/CUTE/Test.cpp | 67 ++++++++++++++++++--- 6 files changed, 258 insertions(+), 18 deletions(-) diff --git a/VERSION b/VERSION index bcaffe1..5db90e1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.7.0 \ No newline at end of file +0.8.0-candidate \ No newline at end of file diff --git a/include/STARE.h.in b/include/STARE.h.in index 1be7eaa..c37e0e3 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -145,6 +145,14 @@ bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId); /// Check if the HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); +STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( + STARE_ArrayIndexSpatialValue spatialStareId, + int resolution, + int shiftAmount + ); + +uint64 spatialLevelMask(); + void STARE_test(); // I added this to simplify finding the library with autoconf. jhrg 5/22/19 diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index 5d97ce3..120ebfe 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -22,11 +22,13 @@ EmbeddedLevelNameEncoding::~EmbeddedLevelNameEncoding() {} * @return */ char* EmbeddedLevelNameEncoding::nameById(uint64 id) { +#ifdef DIAG if(id == 0) { // Throw an exception? std::cout << "EmbeddedLevelNameEncoding::nameById WARNING ID == 0 -> 'S0'..."; // throw SpatialFailure("EmbeddedLevelNameEncoding::nameById-INVALID_ID_0"); } +#endif int nameSize = levelById(id)+3; ///< levelById is local to the encoding char *returnedName = new char[nameSize]; @@ -248,7 +250,7 @@ void EmbeddedLevelNameEncoding::setIdFromSciDBLeftJustifiedFormat( int64 id_scid uint64 level = id_scidb & levelMaskSciDB; iTmp = iTmp << 1; iTmp = iTmp | level; - iTmp = iTmp | TopBit; + iTmp = iTmp | TopBit; // This says we have a valid ID. How does this relate to HTM-ID? // this->id = iTmp; // hacking... this->setId(iTmp); } @@ -365,6 +367,9 @@ bool EmbeddedLevelNameEncoding::terminatorp(uint64 terminatorCandidate) { uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int n) const { /// TODO Error checking of overflow not trustworthy here. using namespace std; + + bool topBitSet = (lowerBound & TopBit) == TopBit; // Some left-justified may not have a topbit set. + uint64 successor = lowerBound; // Bump up one, but we still need the level. // Should clean up successor just in case terminator non-3 prefix is not consistent with level. uint64 one_mask_to_level = 0; @@ -378,17 +383,56 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int } // one_at_level = one_at_level >> 2; - successor = successor & (~one_mask_to_level); +#ifdef DIAG + cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; + cout << "successor: " << setw(23) << hex << successor << endl << flush; +#endif + + // Remove TopBit to test for overflow. + successor = successor & (~TopBit); + + // Truncate to level + successor = successor & (~one_mask_to_level); + + // Increment successor += n*one_at_level; +#ifdef DIAG + cout << "successor: " << setw(23) << hex << successor << endl << flush; + cout << "n: " << setw(23) << dec << n << endl << flush; + cout << "level: " << setw(23) << dec << level << endl << flush; + cout << "n*one_at_level: " << setw(23) << hex << n*one_at_level << endl << flush; + cout << "one_at_level: " << setw(23) << hex << one_at_level << endl << flush; + cout << "successor&TB: " << setw(23) << hex << (successor & TopBit) << endl << flush; + cout << "TopBit: " << setw(23) << hex << TopBit << endl << flush; +#endif + // cout << "one_at_level: "<< hex << one_at_level << dec << endl << flush; // cout << "one_mask_to_: "<< hex << one_mask_to_level << dec << endl << flush; - // Check for overflow. - if( successor == TopBit ) { - return 0; // It's invalid! + // Check for overflow. Not a completely accurate check. + if( (successor & TopBit) == TopBit ) { +#ifdef DIAG + cout << "n: " << setw(23) << dec << n << endl << flush; + cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; + cout << "level: " << setw(23) << dec << level << endl << flush; + cout << "successor: " << setw(23) << hex << successor << endl << flush; + cout << "n*one_at_level: " << setw(23) << hex << n*one_at_level << endl << flush; + cout << "one_at_level: " << setw(23) << hex << one_at_level << endl << flush; + cout << "TopBit: " << setw(23) << hex << TopBit << endl << flush; +#endif + + throw SpatialFailure("EmbeddedLevelNameEncoding::error-increment-overflow"); + // TODO THROW EXCEPTION INSTEAD! + // cout << "ERROR 1" << endl << flush; + // cout << "successor " << successor << ", topbit " << TopBit << endl << flush; + // return 0; // It's invalid! } + // Add TopBit back in if necessary. + if( topBitSet ) { + successor = successor | TopBit; + } successor += level; @@ -403,9 +447,11 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int // if( successor < (( lowerBound & stripMask ) + level) ) { // if( (successor & ~levelMask) < (( lowerBound & ~levelMask )) ) { if( (successor & stripMask) < ( lowerBound & stripMask ) ) { - return 0; // It's invalid! Wrap around! + throw SpatialFailure("EmbeddedLevelNameEncoding::error-increment-wrap-around"); + // TODO THROW EXCEPTION INSTEAD! + // cout << "ERROR 2" << endl << flush; + // return 0; // It's invalid! Wrap around! } - return successor; } uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int n) const { @@ -432,7 +478,9 @@ uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int successor -= n*one_at_level; if( successor == 0 ) { - return 0; // It's invalid! + throw SpatialFailure("EmbeddedLevelNameEncoding::error-decrement-underflow"); + // TODO THROW EXCEPTION INSTEAD! + // return 0; // It's invalid! } successor += level; @@ -453,8 +501,10 @@ uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int // Check for underflow if( (successor & stripMask) > (lowerBound & stripMask)) { + throw SpatialFailure("EmbeddedLevelNameEncoding::error-decrement-wrap-around"); // if( successor > (lowerBound & ~levelMask)+level) { - return 0; // It's invalid! Wrap around! + // TODO THROW EXCEPTION INSTEAD! + // return 0; // It's invalid! Wrap around! } return successor; } diff --git a/src/STARE.C b/src/STARE.C index 3a25943..d4691bb 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -348,6 +348,9 @@ STARE_SpatialIntervals STARE::CoverBoundingBoxFromLatLonDegrees( id1 = EmbeddedLevelNameEncoding(BitShiftNameEncoding(hi).leftJustifiedId()).getSciDBTerminatorLeftJustifiedFormat(); intervals.push_back(id1); } + + + } while( r.getNext(lo,hi) ); } return intervals; @@ -621,3 +624,29 @@ HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals) { return range; } +STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( + STARE_ArrayIndexSpatialValue spatialStareId, + int resolution, + int shiftAmount + ) { +// if( shiftAmount == 0 ) { +// return spatialStareId +// } else { + EmbeddedLevelNameEncoding leftJustifiedWithResolution; + leftJustifiedWithResolution.setIdFromSciDBLeftJustifiedFormat(spatialStareId); + + if( shiftAmount >= 0) { + leftJustifiedWithResolution.setId( + leftJustifiedWithResolution.increment(leftJustifiedWithResolution.getId(), resolution, shiftAmount)); + } else { + leftJustifiedWithResolution.setId( + leftJustifiedWithResolution.decrement(leftJustifiedWithResolution.getId(), resolution, -shiftAmount)); + } + return leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); +// } +} + +uint64 spatialLevelMask() { + EmbeddedLevelNameEncoding leftJustified; + return leftJustified.levelMaskSciDB; +} diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index 8522420..57db842 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -908,6 +908,106 @@ void STARE_test() { // cout << " idx=0, latlon0: " << latlon0.lat << " " << latlon0.lon << endl << flush; // } -// FAIL(); + + if(false) { + LatLonDegrees64 latlon(32.5735,-100.05); + STARE index2; + int resolution = 1; + STARE_ArrayIndexSpatialValue idx = index2.ValueFromLatLonDegrees(latlon.lat,latlon.lon,resolution); + + idx = shiftSpatialIdAtLevel( idx, resolution, 1 ); + + cout << endl << "base..." << endl << flush; + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx << dec + << scientific + << " latlon = " + << latlon.lat << "," << latlon.lon + << dec << " rLevel = " << resolution + << endl << flush; + + + cout << endl << "increasing resolution" << endl << flush; + for(int resolution_ = 0; resolution_ < 12; ++ resolution_) { + STARE_ArrayIndexSpatialValue idx1 = shiftSpatialIdAtLevel( idx, resolution_, 1 ); + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + int resolution1 = idx1 & spatialLevelMask(); + + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx1 << dec + << scientific + << " latlon = " + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + } + + cout << endl << "decreasing resolution" << endl << flush; + for(int resolution_ = 0; resolution_ < 12; ++ resolution_) { + STARE_ArrayIndexSpatialValue idx1 = shiftSpatialIdAtLevel( idx, resolution_, -1 ); + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + int resolution1 = idx1 & spatialLevelMask(); + + cout + << setprecision(18) + << setw(23) + << " idx = 0x" << hex << idx1 << dec + << scientific + << " latlon = " + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + } + + { + cout << endl << "incrementing at resolution" << endl << flush; + uint64 resolution_ = 0; + STARE_ArrayIndexSpatialValue idx1 = 0x0000000000000001; + stringstream ss; + + for( resolution_ = 0; resolution_ < 8; ++resolution_) { + for(uint64 itmp = 0; itmp < 7; ++itmp) { + + // cout << resolution_ << "," << itmp << " res,i" << endl << flush; + + try { + idx1 = shiftSpatialIdAtLevel( 0x0000000000000000+resolution_, resolution_, itmp ); + // idx1 = shiftSpatialIdAtLevel( 0ul+resolution_, resolution_, itmp ); + } catch (const SpatialException & e ) { + cout << "Exception: " << dec << e.what() << endl << flush; + FAIL(); + } + + LatLonDegrees64 latlon1 = index2.LatLonDegreesFromValue(idx1); + + int resolution1 = idx1 & spatialLevelMask(); + + ss.clear(); ss.str(string()); + ss << "Increment by " << itmp << " at resolution level " << resolution_; + + cout + << " idx = 0x" << setw(16) << hex << idx1 << dec + << " cmp = 0x" << setw(16) << hex << (itmp << 59) << dec + << scientific << setprecision(18) << setw(23) + << " latlon = " + << setfill(' ') + << latlon1.lat << "," << latlon1.lon + << dec << " rLevel = " << resolution1 + << " resLevel = " << resolution_ + << endl << flush; + + ASSERT_EQUALM(ss.str().c_str(),((itmp << (59 - 2*resolution_)) || resolution_),idx1); + + } + } + } + } + // FAIL(); } diff --git a/tests/CUTE/Test.cpp b/tests/CUTE/Test.cpp index 89371e7..804aa83 100644 --- a/tests/CUTE/Test.cpp +++ b/tests/CUTE/Test.cpp @@ -1589,22 +1589,59 @@ void htmRangeMultiLevel() { A = rangeFromSymbols("N0001", "N0030"); B = rangeFromSymbols("N0002", "N0030"); +#ifdef DIAG +#define hexOut(a,b) cout << " 0x" << hex << setfill('0') << setw(16) << a << " 0x" << hex << setfill('0') << setw(16) << b << dec << endl << flush; + cout << "A "; hexOut(A.lo,A.hi); + cout << "B "; hexOut(B.lo,B.hi); + cout << "A+"; hexOut(leftJustified.increment(A.lo,3),0); + cout << "B-"; hexOut(leftJustified.decrement(B.lo,3),0); +#undef hexOut +#endif + ASSERT_EQUAL(B.lo,leftJustified.increment(A.lo,3)); ASSERT_EQUAL(A.lo,leftJustified.decrement(B.lo,3)); ASSERT_EQUAL(A.lo,leftJustified.increment(leftJustified.decrement(A.lo,3),3)); A = rangeFromSymbols("S0000", "N0030"); - ASSERT_EQUAL(0,leftJustified.decrement(A.lo,3)); + string failureMessage; uint64 ljx = 0; + failureMessage = "'"; + try { + ljx = leftJustified.decrement(A.lo,3); + } catch ( SpatialFailure& failure ) { + failureMessage += failure.what(); + } + failureMessage += "'"; + ASSERT_EQUALM("Underflow from S0000","'EmbeddedLevelNameEncoding::error-decrement-wrap-around'",failureMessage); + // ASSERT_EQUAL(0,leftJustified.decrement(A.lo,3)); A = rangeFromSymbols("N3333333333333333333333333333","N3333333333333333333333333333"); // Level 27 + failureMessage = "'"; + try { + ljx = leftJustified.increment(A.lo,27); + } catch ( SpatialFailure& failure ) { + failureMessage += failure.what(); + } + failureMessage += "'"; + ASSERT_EQUALM("++A.lo,27","'EmbeddedLevelNameEncoding::error-increment-overflow'",failureMessage); + //#define hexOut(a,b) cout << a << "0x" << hex << setfill('0') << setw(16) << b << dec << endl << flush; // cout << "level(A.lo) " << leftJustified.levelById(A.lo) << endl << flush; // hexOut("A.lo ",A.lo); //#undef hexOut - ASSERT_EQUAL(0,leftJustified.increment(A.lo,27)); + // ASSERT_EQUAL(0,leftJustified.increment(A.lo,27)); A = rangeFromSymbols("N3333","N3333"); - ASSERT_EQUAL(0,leftJustified.increment(A.lo,3)); + + failureMessage = "'"; + try { + ljx = leftJustified.increment(A.lo,3); + } catch ( SpatialFailure& failure ) { + failureMessage += failure.what(); + } + failureMessage += "'"; + ASSERT_EQUALM("++A.lo,3","'EmbeddedLevelNameEncoding::error-increment-overflow'",failureMessage); + + // ASSERT_EQUAL(0,leftJustified.increment(A.lo,3)); A = rangeFromSymbols("N33000","N33330"); B = rangeFromSymbols("N3301","N3333"); @@ -3964,9 +4001,18 @@ void LeftJustifiedDecrementBug() { // cout << "ljx: " << hex << ljx << dec << endl << flush; ASSERT_EQUALM("Max left inc2 dec2",0xfff0000000000004,ljx); - ljx = lj.increment(ljx,4); + string failureMessage = "'"; + try { + ljx = lj.increment(ljx,4); + } catch ( SpatialFailure& failure ) { + failureMessage += failure.what(); + } + failureMessage += "'"; // cout << "ljx: " << hex << ljx << dec << endl << flush; - ASSERT_EQUALM("Max left inc3 dec2",0,ljx); + ASSERT_EQUALM("Max left inc3 dec2","'EmbeddedLevelNameEncoding::error-increment-overflow'",failureMessage); + + // 0 Used to be an error return code. + // ASSERT_EQUALM("Max left inc3 dec2",0,ljx); // cout << "--" << endl; @@ -3990,9 +4036,16 @@ void LeftJustifiedDecrementBug() { // cout << "ljx: " << hex << ljx << dec << endl << flush; ASSERT_EQUALM("Min left dec2 inc2",0x8000000000000004,ljx); - ljx = lj.decrement(ljx,4); + failureMessage = "'"; + try { + ljx = lj.decrement(ljx,4); + } catch ( SpatialFailure& failure ) { + failureMessage += failure.what(); + } + failureMessage += "'"; // cout << "ljx: " << hex << ljx << dec << endl << flush; - ASSERT_EQUALM("Min left dec3 inc2",0,ljx); + // ASSERT_EQUALM("Min left dec3 inc2",0,ljx); + ASSERT_EQUALM("Min left dec3 inc2","'EmbeddedLevelNameEncoding::error-decrement-wrap-around'",failureMessage); // FAIL(); } From d0481156caa4763e89bb1491a82c1a8845c3d393 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Mon, 9 Sep 2019 13:40:41 -0400 Subject: [PATCH 02/31] Initial commit of new SpatialRange. SpatialRange provides a way for STARE applications to reach into and use the multi-level HstmRange code. There's overhead in setting up the range and in translation back and forth between STARE format and the underlying native left-justified format. --- include/CMakeLists.txt | 3 +- include/EmbeddedLevelNameEncoding.h | 15 +-- include/HstmRange.h | 4 +- include/SpatialRange.h | 59 ++++++++++ src/CMakeLists.txt | 1 + src/EmbeddedLevelNameEncoding.C | 30 +++-- src/HstmRange.C | 41 +++++-- src/HtmRangeMultiLevel.cpp | 2 + src/STARE.C | 27 +---- src/SpatialRange.cpp | 109 ++++++++++++++++++ tests/CUTE/CMakeLists.txt | 2 + tests/CUTE/EmbeddedLevelNameEncoding_test.cpp | 58 ++++++++++ tests/CUTE/SpatialRange_test.cpp | 81 +++++++++++++ tests/CUTE/Test.cpp | 2 + tests/CUTE/Test.h | 2 +- 15 files changed, 387 insertions(+), 49 deletions(-) create mode 100644 include/SpatialRange.h create mode 100644 src/SpatialRange.cpp create mode 100644 tests/CUTE/EmbeddedLevelNameEncoding_test.cpp create mode 100644 tests/CUTE/SpatialRange_test.cpp diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index 352c1f1..0f16164 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -33,9 +33,10 @@ set ( ${CMAKE_CURRENT_SOURCE_DIR}/SpatialIndex.hxx ${CMAKE_CURRENT_SOURCE_DIR}/SpatialInterface.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialInterface.hxx - ${CMAKE_CURRENT_SOURCE_DIR}/SpatialSign.h + ${CMAKE_CURRENT_SOURCE_DIR}/SpatialSign.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialVector.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialVector.hxx + ${CMAKE_CURRENT_SOURCE_DIR}/SpatialRange.h ${CMAKE_CURRENT_SOURCE_DIR}/SpatialRotation.h ${CMAKE_CURRENT_SOURCE_DIR}/STARE.h ${CMAKE_CURRENT_SOURCE_DIR}/TemporalGeneral.h diff --git a/include/EmbeddedLevelNameEncoding.h b/include/EmbeddedLevelNameEncoding.h index d5adb91..1fc59ec 100644 --- a/include/EmbeddedLevelNameEncoding.h +++ b/include/EmbeddedLevelNameEncoding.h @@ -63,8 +63,12 @@ class EmbeddedLevelNameEncoding: public NameEncoding { uint32 aLevel ) const; ///< Somewhat dangerous to use. uint64 idFromTerminatorAndLevel_NoDepthBit(uint64 terminator, uint32 level); ///< Also a little dangerous. - bool terminatorp(); - bool terminatorp(uint64 terminatorCandidate); + + bool terminatorp() const; + bool terminatorp(uint64 terminatorCandidate) const; + + bool SciDBterminatorp() const; + bool SciDBterminatorp(uint64 terminatorCandidate) const; /// What triangle is just after the terminator? uint64 successorToTerminator_NoDepthBit(uint64 terminator, uint32 level) const; @@ -101,11 +105,6 @@ class EmbeddedLevelNameEncoding: public NameEncoding { uint64 increment(uint64 lowerBound, uint32 level, int steps = 1) const; uint64 decrement(uint64 lowerBound, uint32 level, int steps = 1) const; - bool terminatorp(uint64 terminator) const { - uint64 levelBits = terminator & levelMask; - return levelBits == levelMask; - } - EmbeddedLevelNameEncoding& operator=(EmbeddedLevelNameEncoding obj) { this->setId(obj.getId()); return *this; @@ -222,4 +221,6 @@ class EmbeddedLevelNameEncoding: public NameEncoding { }; +void EmbeddedLevelNameEncoding_test(); + #endif /* INCLUDE_EMBEDDEDLEVELNAMEENCODING_H_ */ diff --git a/include/HstmRange.h b/include/HstmRange.h index ec55f72..8d6db3e 100644 --- a/include/HstmRange.h +++ b/include/HstmRange.h @@ -36,6 +36,8 @@ namespace std { class HstmRange { public: HstmRange(); + HstmRange(HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel *range); + // TOO Deep copy? HstmRange(HstmRange range); virtual ~HstmRange(); // TODO Note the int in the following is a return code, not an index. @@ -53,8 +55,6 @@ class HstmRange { /** * For STARE we may need to coarsen position information to the resolution level since it points to a 7cm triangle inside a larger triangle. * - * - * */ void addRange(Key a, Key b); void addRange(HstmRange* r); diff --git a/include/SpatialRange.h b/include/SpatialRange.h new file mode 100644 index 0000000..5cb73f9 --- /dev/null +++ b/include/SpatialRange.h @@ -0,0 +1,59 @@ +/* + * SpatialRange.h + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#ifndef INCLUDE_SPATIALRANGE_H_ +#define INCLUDE_SPATIALRANGE_H_ + +#include "STARE.h" +#include "HstmRange.h" + +namespace std { + +/** + * A wrapper for an HstmRange that knows about STARE to provide htm-range-like functions. + * + */ +class SpatialRange { +public: + SpatialRange(); + SpatialRange(STARE_SpatialIntervals intervals); + SpatialRange(HstmRange *range) { this->range = range; } + virtual ~SpatialRange(); + + void addSpatialIntervals(STARE_SpatialIntervals intervals); + void addSpatialRange(SpatialRange r); + + STARE_SpatialIntervals toSpatialIntervals(); + int getNextSpatialInterval(STARE_SpatialIntervals &interval); + + HstmRange *range; + + /////////// private: Maybe? //////////// + + // Maybe inherit these? + // TODO Note the int in the following is a return code, not an index. + int getNext(KeyPair &kp) { + int istat = range->getNext(kp); +// cout << "" << flush; + return istat; + }; + void reset() { range->reset(); } // range not null? + void purge() { range->purge(); } // what if range null? + +}; + +inline SpatialRange operator& ( const SpatialRange& a, const SpatialRange& b) { + return SpatialRange(new HstmRange(a.range->range->HtmRangeMultiLevelAtLevelFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. +} + +void SpatialRange_test(); + +} /* namespace std */ + +#endif /* INCLUDE_SPATIALRANGE_H_ */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f636b08..f44aded 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,7 @@ set ( SpatialInterface.cpp SpatiallyAdaptiveDataCover.C SpatialVector.cpp + SpatialRange.cpp SpatialRotation.C STARE.C TemporalIndex.cpp diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index 120ebfe..2371b26 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -192,7 +192,7 @@ uint64 EmbeddedLevelNameEncoding::idFromTerminatorAndLevel_NoDepthBit(uint64 ter int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat(uint64 leftId) const { uint64 id_NoLevelBit = leftId & stripLevelBitMask; // Note this covers bits 0-5. - uint64 level = leftId & levelMask; + uint64 level = leftId & levelMask; // TODO Repent the sin of redundant code. int64 leftId_scidb = id_NoLevelBit; @@ -223,7 +223,7 @@ int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat() const { // *** TODO *** NEED ERROR CHECKING & SIGNALLING FOR INVALID IDs - int64 leftId_scidb = 0; // Erf... // What's the invalid id now? Anything negative. + int64 leftId_scidb = 0; // Erf... // What's the invalid id now? Anything negative? uint64 id_NoLevelBit = this->maskOffLevelBit(); int64 level = this->getLevel(); @@ -240,6 +240,9 @@ int64 EmbeddedLevelNameEncoding::getSciDBLeftJustifiedFormat() const { } int64 EmbeddedLevelNameEncoding::getSciDBTerminatorLeftJustifiedFormat() const { + if(terminatorp()) { + return getSciDBLeftJustifiedFormat(); // Just return what we have. + } return getSciDBLeftJustifiedFormat(getIdTerminator_NoDepthBit()); } @@ -252,6 +255,10 @@ void EmbeddedLevelNameEncoding::setIdFromSciDBLeftJustifiedFormat( int64 id_scid iTmp = iTmp | level; iTmp = iTmp | TopBit; // This says we have a valid ID. How does this relate to HTM-ID? // this->id = iTmp; // hacking... + // if(SciDBterminatorp(id_scidb)) { + if( level == levelMaskSciDB ) { + iTmp = iTmp | levelMask; // Fill in the extra bit. + } this->setId(iTmp); } @@ -355,15 +362,20 @@ uint64 EmbeddedLevelNameEncoding::predecessorToLowerBound_NoDepthBit(uint64 lowe return terminator; } -bool EmbeddedLevelNameEncoding::terminatorp() { +bool EmbeddedLevelNameEncoding::terminatorp() const { return terminatorp(this->id); } - -bool EmbeddedLevelNameEncoding::terminatorp(uint64 terminatorCandidate) { - uint64 level = terminatorCandidate & levelMask; - return level == 63; +bool EmbeddedLevelNameEncoding::terminatorp(uint64 terminator) const { + uint64 levelBits = terminator & levelMask; + return levelBits == levelMask; +} +bool EmbeddedLevelNameEncoding::SciDBterminatorp() const { + return SciDBterminatorp(this->id); +} +bool EmbeddedLevelNameEncoding::SciDBterminatorp(uint64 terminator) const { + uint64 levelBits = terminator & levelMaskSciDB; + return levelBits == levelMaskSciDB; } - uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int n) const { /// TODO Error checking of overflow not trustworthy here. using namespace std; @@ -510,3 +522,5 @@ uint64 EmbeddedLevelNameEncoding::decrement(uint64 lowerBound, uint32 level, int } + + diff --git a/src/HstmRange.C b/src/HstmRange.C index 19ff0ba..7c10654 100644 --- a/src/HstmRange.C +++ b/src/HstmRange.C @@ -6,16 +6,19 @@ */ #include "HstmRange.h" +#include namespace std { HstmRange::HstmRange() { - // TODO Auto-generated constructor stub range = new HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel(); } +HstmRange::HstmRange(HtmRangeMultiLevel_NameSpace::HtmRangeMultiLevel *range){ + this->range = range; +} + HstmRange::~HstmRange() { - // TODO Auto-generated destructor stub delete range; } @@ -35,13 +38,22 @@ void HstmRange::addRange(Key a_, Key b_) { // cout << dec << 902 << " : aLevel = " << aLevel << endl << flush; if( aLevel != bLevel ) { - cout << "HstmRange::addRange::ERROR::KeyLevelMismatch " - << aLevel << " != " << bLevel << " " - << "a = " << a_ << " " - << "b = " << b_ << " " - << endl << flush; - // TODO Throw? - throw SpatialException("HstmRange::addRange::ERROR::KeyLevelMismatch"); + if( !leftJustifiedEncoding.terminatorp(b_) ) { +#ifdef DIAG + cout << "HstmRange::addRange::ERROR::KeyLevelMismatch " + << aLevel << " != " << bLevel << " " + << "a = " << a_ << " " + << "b = " << b_ << " " + << endl << flush; +#endif + // TODO Throw? + throw SpatialException("HstmRange::addRange::ERROR::KeyLevelMismatch"); + } else { + b = b | leftJustifiedEncoding.levelMask; + } + } else { + leftJustifiedEncoding.setId(b_); + b = leftJustifiedEncoding.getIdTerminator_NoDepthBit(); } // Put a into HtmRange without change @@ -65,8 +77,19 @@ void HstmRange::addRange(Key a_, Key b_) { Key bSave = b; */ + /* 2019-0909 MLR THIS PIECE "WORKS". leftJustifiedEncoding.setId(b_); b = leftJustifiedEncoding.getIdTerminator_NoDepthBit(); + */ + +// #define DIAG +#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << " " << s << endl << flush; +#ifdef DIAG + IDOUT(cout,"hr::ar a_: ",a_) + IDOUT(cout,"hr::ar a : ",a) + IDOUT(cout,"hr::ar b_: ",b_) + IDOUT(cout,"hr::ar b : ",b) +#endif // cout << dec << 1000 << " : " << hex << "0x" << a_ << " 0x" << bSave << endl << flush; diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index b21cd40..933a3ab 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -1581,8 +1581,10 @@ int HtmRangeMultiLevel::getNext(Key &lo, Key &hi) // cout << " " << hi << " " << flush; // OLD if (hi <= (Key) 0){ if (hi < (Key) 0){ +#if DIAG cout << endl; cout << " getNext error!! " << endl << flush; +#endif // hi = lo = (Key) 0; hi = lo = (Key) -1; return 0; diff --git a/src/STARE.C b/src/STARE.C index d4691bb..dba1146 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -403,6 +403,12 @@ STARE_SpatialIntervals STARE::CoverCircleFromLatLonRadiusDegrees(float64 latDegr return intervals; } +/** + * Returns a list of STARE spatial IDs on the hull about @param points. + * + * Unfortunately, it erroneously returns a few spatial IDs from the interior as well. + * + */ STARE_SpatialIntervals STARE::ConvexHull(LatLonDegrees64ValueVector points,int force_resolution_level) { STARE_SpatialIntervals cover; @@ -601,28 +607,7 @@ bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId) { return leftJustifiedWithResolution.terminatorp(); } -/** - * Construct a range object (spatial region) from a vector of intervals. - * - * TODO: Have a SpatialRange class instead of an HstmRange? - */ -HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals) { - HstmRange range; - EmbeddedLevelNameEncoding leftJustified; - for(auto i0=intervals.begin(); i0 != intervals.end(); ++i0) { - leftJustified.setIdFromSciDBLeftJustifiedFormat(*i0); - uint64 a = leftJustified.getId(), b = a; - auto i1 = (i0+1); - if(terminatorp(*i1)) { - leftJustified.setIdFromSciDBLeftJustifiedFormat(*i1); - b = leftJustified.getId(); - ++i0; // Skip to next - } - range.addRange(a,b); - } - return range; -} STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( STARE_ArrayIndexSpatialValue spatialStareId, diff --git a/src/SpatialRange.cpp b/src/SpatialRange.cpp new file mode 100644 index 0000000..5545545 --- /dev/null +++ b/src/SpatialRange.cpp @@ -0,0 +1,109 @@ +/* + * SpatialRange.cpp + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "SpatialRange.h" + +namespace std { + +SpatialRange::SpatialRange() { + this->range = new HstmRange; +} + +SpatialRange::SpatialRange(STARE_SpatialIntervals intervals) { + this->range = new HstmRange; + this->addSpatialIntervals(intervals); +} + +SpatialRange::~SpatialRange() { + delete this->range; // TODO mlr. Um. Dangerous if range is from elsewhere? Maybe use some sort of smart pointer? +} + +void SpatialRange::addSpatialRange(SpatialRange range) { + this->range->addRange(range.range); +} + +/** + * Construct a range object (spatial region) from a vector of intervals. + * + * TODO: Have a SpatialRange class instead of an HstmRange? + */ +// #define DIAG +void SpatialRange::addSpatialIntervals(STARE_SpatialIntervals intervals) { + EmbeddedLevelNameEncoding leftJustified; + for(auto i0=intervals.begin(); i0 != intervals.end(); ++i0) { + leftJustified.setIdFromSciDBLeftJustifiedFormat(*i0); + uint64 a = leftJustified.getId(), b = a; + auto i1 = (i0+1); + if( i1 <= intervals.end() ) { + if(terminatorp(*i1)) { + leftJustified.setIdFromSciDBLeftJustifiedFormat(*i1); + b = leftJustified.getId(); + ++i0; // Skip to next + } + } + +#ifdef DIAG + cout << "sr::addsi " + << setw(16) << setfill('0') << hex << a << " " + << setw(20) << dec << a << " " + << setw(16) << setfill('0') << hex << b << " " + << setw(20) << dec << b << " " + << endl << flush; +#endif + this->range->addRange(a,b); +#ifdef DIAG + cout << "sr::addsi nr = " << this->range->range->nranges() << endl << flush; +#endif + } +} + + +int SpatialRange::getNextSpatialInterval(STARE_SpatialIntervals &interval) { + KeyPair kp(-1,-2); + int istat = this->getNext(kp); +#ifdef DIAG + cout << "sr::gnsi istat = " << istat << ", kp = " << kp.lo << ", " << kp.hi << endl << flush; +#endif + if( istat > 0 ) { + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setId(kp.lo); + interval.push_back(leftJustified.getSciDBLeftJustifiedFormat()); + if( kp.lo != kp.hi ) { + STARE_ArrayIndexSpatialValue term_lo = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + leftJustified.setId(kp.hi); // Note: hi should be a terminator already. + STARE_ArrayIndexSpatialValue term_hi = leftJustified.getSciDBLeftJustifiedFormat(); + +#ifdef DIAG + cout << "lo,term_lo,term_hi: " + << setw(16) << setfill('0') << hex + << kp.lo << "," + << setw(16) << setfill('0') << hex + << term_lo << "," + << setw(16) << setfill('0') << hex + << term_hi << endl << flush << dec; +#endif + + if( term_hi != term_lo ) { + interval.push_back( term_hi ); + } + } + } + return istat; +} +#undef DIAG + +STARE_SpatialIntervals SpatialRange::toSpatialIntervals() { + STARE_SpatialIntervals intervals; + this->range->reset(); + while( this->getNextSpatialInterval(intervals) > 0 ); + return intervals; +} + + +} /* namespace std */ diff --git a/tests/CUTE/CMakeLists.txt b/tests/CUTE/CMakeLists.txt index 4d4b2e5..202c1ef 100644 --- a/tests/CUTE/CMakeLists.txt +++ b/tests/CUTE/CMakeLists.txt @@ -1,10 +1,12 @@ set ( TestSrcFiles + EmbeddedLevelNameEncoding_test.cpp Test.cpp SpatialFailure_test.cpp SpatialIndex_test.cpp SpatialInterface_test.cpp + SpatialRange_test.cpp SpatialRotation_test.cpp SpatialVector_test.cpp STARE_test.cpp diff --git a/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp new file mode 100644 index 0000000..c4d7d22 --- /dev/null +++ b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp @@ -0,0 +1,58 @@ +/* + * EmbeddedLevelNameEncoding_test.cpp + * + * Created on: Sep 9, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "Test.h" + +void EmbeddedLevelNameEncoding_test() { + STARE index; + STARE_ArrayIndexSpatialValue sid = index.ValueFromLatLonDegrees( 30, 30, 8 ); + + EmbeddedLevelNameEncoding leftJustified; + + leftJustified.setIdFromSciDBLeftJustifiedFormat(sid); +#ifdef DIAG + cout << "sid tp? " << leftJustified.terminatorp() << endl << flush; +#endif + + STARE_ArrayIndexSpatialValue term_sid1 = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + leftJustified.setIdFromSciDBLeftJustifiedFormat(term_sid1); +#ifdef DIAG + cout << "term_sid1 tp? " << leftJustified.terminatorp() << endl << flush; +#endif + + STARE_ArrayIndexSpatialValue term_sid2 = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + +#ifdef DIAG + cout << "test: id,mask,level " + << setw(16) << setfill('0') << hex + << leftJustified.getId() + << " " << dec << leftJustified.levelMask + << " " << (leftJustified.getId() & leftJustified.levelMask) + << endl << flush; + cout << "testS: id,mask,level " + << setw(16) << setfill('0') << hex + << leftJustified.getId() + << " " << dec << leftJustified.levelMaskSciDB + << " " << (leftJustified.getId() & leftJustified.levelMaskSciDB) + << endl << flush; + cout << endl << flush; + +#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << endl << flush; + IDOUT(cout,"sid ", sid); + IDOUT(cout,"term_sid1 ", term_sid1); + IDOUT(cout,"term_sid1*", leftJustified.getId()); + IDOUT(cout,"term_sid1*", leftJustified.getSciDBLeftJustifiedFormat()); + IDOUT(cout,"term_sid1t", leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + IDOUT(cout,"term_sid2 ", term_sid2); +#undef IDOUT +#endif + + ASSERT(leftJustified.terminatorp(term_sid1)); + ASSERT_EQUALM("getSciDBTerminatorLeftJustifiedFormat idempotent?",term_sid1,term_sid2); +} diff --git a/tests/CUTE/SpatialRange_test.cpp b/tests/CUTE/SpatialRange_test.cpp new file mode 100644 index 0000000..99aaf51 --- /dev/null +++ b/tests/CUTE/SpatialRange_test.cpp @@ -0,0 +1,81 @@ +/* + * SpatialRange_test.cpp + * + * Created on: Sep 8, 2019 + * Author: mrilee + * + * Copyright (C) 2019 Rilee Systems Technologies LLC + */ + +#include "Test.h" + +namespace std { + +void SpatialRange_test () { + + STARE index; + STARE_ArrayIndexSpatialValue sid = index.ValueFromLatLonDegrees( 30, 30, 8 ); + STARE_SpatialIntervals sids; sids.push_back(sid); + + SpatialRange range(sids); + STARE_SpatialIntervals sids_out; + + sids_out = range.toSpatialIntervals(); + +// #define DIAG +#ifdef DIAG + for(int i=0; i < sids.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sids[i] << dec << endl << flush; + } + for(int i=0; i < sids_out.size(); ++i ) { + cout << i << " i,so: " << setw(16) << setfill('0') << hex << sids_out[i] << dec << endl << flush; + } +#endif + + ASSERT_EQUAL(1,sids_out.size()); + ASSERT_EQUAL(sids[0],sids_out[0]); + + // sids.push_back(index.ValueFromLatLonDegrees( 45, 45, 8)); + + sid = index.ValueFromLatLonDegrees( 15, 15, 8 ); + + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setIdFromSciDBLeftJustifiedFormat(sid); + sids.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + +#ifdef DIAG + for(int i=0; i < sids.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sids[i] << dec << endl << flush; + } +// for(int i=0; i < sids_out.size(); ++i ) { +// cout << i << " i,so: " << setw(16) << setfill('0') << hex << sids_out[i] << dec << endl << flush; +// } +#endif + + // cout << 50 << endl << flush; + + range.purge(); + range.addSpatialIntervals(sids); +// cout << 100 +// << " nr: " << range.range->range->nranges() +// << endl << flush; + sids_out = range.toSpatialIntervals(); + +#ifdef DIAG + cout << 110 << " sids_out.s: " << sids_out.size() << endl << flush; + for(int i=0; i Date: Tue, 10 Sep 2019 23:37:39 -0400 Subject: [PATCH 03/31] Adding a SpatialRange and expandIntervals. SpatialRange provides a STARE-based way to get at HstmRange. The expandIntervals routines expand intervals into individual trixels. This provides a way to force intervals to a single resolution. Acts as a conversion from STARE_SpatialIntervals to STARE_ArrayIndexSpatialValues. The latter has no intervals. --- include/CMakeLists.txt | 1 + include/EmbeddedLevelNameEncoding.h | 1 + include/STARE.h.in | 5 ++ include/SpatialRange.h | 3 - src/CMakeLists.txt | 2 +- src/EmbeddedLevelNameEncoding.C | 20 ++++++ src/STARE.C | 100 +++++++++++++++++++++++++++- src/SpatialRange.cpp | 2 - tests/CUTE/CMakeLists.txt | 4 +- tests/CUTE/STARE_test.cpp | 65 ++++++++++++++++++ tests/CUTE/SpatialRange_test.cpp | 3 - tests/CUTE/Test.h | 3 +- 12 files changed, 195 insertions(+), 14 deletions(-) diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index 0f16164..b1a5c45 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -86,6 +86,7 @@ install(FILES STARE.h SpatialIndex.hxx SpatialInterface.h SpatialInterface.hxx + SpatialRange.h SpatialSign.h SpatialVector.h SpatialVector.hxx diff --git a/include/EmbeddedLevelNameEncoding.h b/include/EmbeddedLevelNameEncoding.h index 1fc59ec..4853b1b 100644 --- a/include/EmbeddedLevelNameEncoding.h +++ b/include/EmbeddedLevelNameEncoding.h @@ -102,6 +102,7 @@ class EmbeddedLevelNameEncoding: public NameEncoding { EmbeddedLevelNameEncoding clearDeeperThanLevel(uint64 level); + void increment_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const; uint64 increment(uint64 lowerBound, uint32 level, int steps = 1) const; uint64 decrement(uint64 lowerBound, uint32 level, int steps = 1) const; diff --git a/include/STARE.h.in b/include/STARE.h.in index c37e0e3..1cf72f6 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -145,6 +145,11 @@ bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId); /// Check if the HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); +void STARE_ArrayIndexSpatialValues_insert(STARE_ArrayIndexSpatialValues& v, STARE_ArrayIndexSpatialValue siv); + +STARE_SpatialIntervals expandInterval(STARE_SpatialIntervals interval, int64 force_resolution = -1); +STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, int64 force_resolution = -1); + STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( STARE_ArrayIndexSpatialValue spatialStareId, int resolution, diff --git a/include/SpatialRange.h b/include/SpatialRange.h index 5cb73f9..f612813 100644 --- a/include/SpatialRange.h +++ b/include/SpatialRange.h @@ -13,8 +13,6 @@ #include "STARE.h" #include "HstmRange.h" -namespace std { - /** * A wrapper for an HstmRange that knows about STARE to provide htm-range-like functions. * @@ -54,6 +52,5 @@ inline SpatialRange operator& ( const SpatialRange& a, const SpatialRange& b) { void SpatialRange_test(); -} /* namespace std */ #endif /* INCLUDE_SPATIALRANGE_H_ */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f44aded..d16f5c8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,9 +24,9 @@ set ( SpatialInterface.cpp SpatiallyAdaptiveDataCover.C SpatialVector.cpp - SpatialRange.cpp SpatialRotation.C STARE.C + SpatialRange.cpp TemporalIndex.cpp TemporalWordFormat1.cpp VarStr.C diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index 2371b26..422c402 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -376,6 +376,20 @@ bool EmbeddedLevelNameEncoding::SciDBterminatorp(uint64 terminator) const { uint64 levelBits = terminator & levelMaskSciDB; return levelBits == levelMaskSciDB; } + +void EmbeddedLevelNameEncoding::increment_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const { + one_mask_to_level = 0; + one_at_level = one; + for(uint64 shift = 2; + shift <= ( (topBitPosition-3) - 2*(level) ); + shift +=2 ) { + one_mask_to_level = one_mask_to_level << 2; + one_mask_to_level += 3; + one_at_level = one_at_level << 2; + } +// one_at_level = one_at_level >> 2; +} + uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int n) const { /// TODO Error checking of overflow not trustworthy here. using namespace std; @@ -384,6 +398,8 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int uint64 successor = lowerBound; // Bump up one, but we still need the level. // Should clean up successor just in case terminator non-3 prefix is not consistent with level. + + /* uint64 one_mask_to_level = 0; uint64 one_at_level = one; for(uint64 shift = 2; @@ -394,6 +410,10 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int one_at_level = one_at_level << 2; } // one_at_level = one_at_level >> 2; +*/ + + uint64 one_mask_to_level, one_at_level; + increment_LevelToMaskDelta(level,one_mask_to_level,one_at_level); #ifdef DIAG cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; diff --git a/src/STARE.C b/src/STARE.C index dba1146..22f4d35 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -14,6 +14,7 @@ #include "SpatialDomain.h" #include "SpatialInterface.h" #include +#include /** * @brief Version function with C linkage to aid in finding the library with autoconf @@ -22,6 +23,21 @@ extern "C" const char *STARE_version() { return (const char *)STARE_VERSION; } + +void STARE_ArrayIndexSpatialValues_insert(STARE_ArrayIndexSpatialValues& v, STARE_ArrayIndexSpatialValue siv) { + // cout << dec << 1000 << endl; + STARE_ArrayIndexSpatialValues::iterator + it = std::lower_bound( v.begin(), v.end(), siv); + // cout << dec << 1010 << endl; + if(it == v.end()) { + v.insert( it, siv ); + } else if( (*it) != siv) { + // cout << dec << 1020 << endl; + v.insert( it, siv ); + // cout << dec << 1030 << endl; + } + // cout << dec << 1040 << endl; +} /** * @@ -250,6 +266,9 @@ STARE_ArrayIndexSpatialValue sTerminator(STARE_ArrayIndexSpatialValue spatialSta /** * Compare two spatial array index values a, b. * + * Suffers some overhead, but should handle the overlap. For other + * kinds of intersection we need to handle intervals or move to a range. + * * Returns * 1 if b is in a * -1 if a is in b @@ -607,8 +626,6 @@ bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId) { return leftJustifiedWithResolution.terminatorp(); } - - STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( STARE_ArrayIndexSpatialValue spatialStareId, int resolution, @@ -635,3 +652,82 @@ uint64 spatialLevelMask() { EmbeddedLevelNameEncoding leftJustified; return leftJustified.levelMaskSciDB; } + +STARE_ArrayIndexSpatialValues expandInterval(STARE_SpatialIntervals interval, int64 force_resolution) { + // STARE_SpatialIntervals interval should just be one interval, i.e. a value or value+terminator. + // cout << dec << 200 << endl << flush; + STARE_ArrayIndexSpatialValue siv0 = interval[0]; + EmbeddedLevelNameEncoding leftJustified; + // cout << dec << 220 << endl << flush; + uint64 return_resolution = siv0 & leftJustified.levelMaskSciDB; + // cout << dec << 225 << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush; + if( force_resolution > -1 ) { + siv0 = ( siv0 & ~leftJustified.levelMaskSciDB ) | force_resolution; + return_resolution = force_resolution; + } + // cout << dec << 230 << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); + // cout << dec << 235 << endl << flush; + STARE_ArrayIndexSpatialValue siv_term; + if( interval.size() > 1 ) { + siv_term = interval[1]; + } else { + siv_term = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + } + // cout << dec << 240 << endl << flush; + uint64 one_mask_to_resolution, one_at_resolution; + leftJustified.increment_LevelToMaskDelta(siv0 & leftJustified.levelMaskSciDB,one_mask_to_resolution,one_at_resolution); + // cout << dec << 245 << endl << flush; + // Give as much rope as needed. + siv0 = (siv0 & ~one_mask_to_resolution) | return_resolution; + STARE_ArrayIndexSpatialValues expanded_interval; + while( siv0 < siv_term ) { + expanded_interval.push_back(siv0); + siv0 += one_at_resolution; + } + // cout << dec << 250 << endl << flush; + return expanded_interval; +} + +/** + * + */ +STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, int64 force_resolution) { + STARE_ArrayIndexSpatialValues expanded_values; + EmbeddedLevelNameEncoding leftJustified; + + int i=0; + while( i < intervals.size() ) { + // cout << dec << 100 << endl << flush; + STARE_ArrayIndexSpatialValue siv0, siv1; + STARE_SpatialIntervals interval; + siv0 = intervals[i]; + interval.push_back(siv0); + // cout << dec << 110 << endl << flush; + ++i; + if( i < intervals.size() ) { + // peek + // cout << dec << 120 << endl << flush; + siv1 = intervals[i]; + // cout << dec << 120 << " " << i << " " << setw(16) << setfill('0') << hex << siv1 << dec << endl << flush; + if( (siv1 & leftJustified.levelMaskSciDB) == leftJustified.levelMaskSciDB ) { + // cout << dec << 121 << " " << i << endl << flush; + interval.push_back(siv1); + ++i; + } + } + // cout << dec << 130 << " " << i << endl << flush; + STARE_SpatialIntervals expandOne = expandInterval(interval,force_resolution); + // cout << dec << 140 << " " << i << endl << flush; + for(int j=0; j < expandOne.size(); ++j) { +// cout << dec << 142 << " " << i << " " << j +// << setw(16) << setfill('0') << hex << expandOne[j] << dec +// << endl << flush; + STARE_ArrayIndexSpatialValues_insert( expanded_values, expandOne[j] ); +// cout << dec << 143 << " " << i << " " << j << endl << flush; + } +// cout << dec << 150 << " " << i << endl << flush; + } +// cout << dec << 160 << " " << i << endl << flush; + return expanded_values; +} diff --git a/src/SpatialRange.cpp b/src/SpatialRange.cpp index 5545545..ebd5220 100644 --- a/src/SpatialRange.cpp +++ b/src/SpatialRange.cpp @@ -9,7 +9,6 @@ #include "SpatialRange.h" -namespace std { SpatialRange::SpatialRange() { this->range = new HstmRange; @@ -106,4 +105,3 @@ STARE_SpatialIntervals SpatialRange::toSpatialIntervals() { } -} /* namespace std */ diff --git a/tests/CUTE/CMakeLists.txt b/tests/CUTE/CMakeLists.txt index 202c1ef..f8b4089 100644 --- a/tests/CUTE/CMakeLists.txt +++ b/tests/CUTE/CMakeLists.txt @@ -1,8 +1,8 @@ set ( TestSrcFiles + Test.h EmbeddedLevelNameEncoding_test.cpp - Test.cpp SpatialFailure_test.cpp SpatialIndex_test.cpp SpatialInterface_test.cpp @@ -11,7 +11,7 @@ set ( SpatialVector_test.cpp STARE_test.cpp TemporalIndex_test.cpp - Test.h + Test.cpp ) # SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -I${PROJECT_SOURCE_DIR}/include" ) diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index 57db842..8a4f93b 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -1008,6 +1008,71 @@ void STARE_test() { } } } + +#define SIVOUT(m,siv) +// #define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; + if(true) { + // cout << "expandIntervals 10" << endl << flush; + EmbeddedLevelNameEncoding leftJustified; + uint64 one_mask_to_level, one_at_level; + leftJustified.increment_LevelToMaskDelta(8,one_mask_to_level, one_at_level); + + // cout << "expandIntervals 20" << endl << flush; + + // STARE Spatial index value and interval arrays. + // STARE_ArrayIndexSpatialValue siv0 = 0x3300000000000008; + STARE_ArrayIndexSpatialValue siv0 = 0x3213213213213208; + SIVOUT("siv0",siv0) + + STARE_SpatialIntervals intervals; + intervals.push_back(siv0); + + // cout << "expandIntervals 30" << endl << flush; + + STARE_ArrayIndexSpatialValues expanded_values = expandIntervals(intervals,8); + ASSERT_EQUAL((siv0 & ~one_mask_to_level) | 8,expanded_values[0]); + + // cout << "expandIntervals 40" << endl << flush; + + intervals.push_back(0x2000000000000008); + expanded_values = expandIntervals(intervals,8); + + ASSERT_EQUAL(intervals[1],expanded_values[0]); + + siv0 = 0x00000000000000008; + intervals.push_back(siv0); + siv0+=6*one_at_level; + SIVOUT("siv0+6@8",siv0); + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); + intervals.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + SIVOUT("siv0+6@8",leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + STARE_ArrayIndexSpatialValue expected6_term = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); + + expanded_values = expandIntervals(intervals,8); + + ASSERT_EQUAL(0x0000300000000008,expanded_values[3]); + + leftJustified.setIdFromSciDBLeftJustifiedFormat(expanded_values[6]); + SIVOUT("6's ",leftJustified.getSciDBLeftJustifiedFormat()); + SIVOUT("6's term",leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + ASSERT_EQUAL(expected6_term,leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + + expanded_values = expandIntervals(intervals,7); + +#ifdef DIAG + cout << "intervals "<< dec << endl < Date: Thu, 12 Sep 2019 11:50:48 -0400 Subject: [PATCH 04/31] Initial commit of SpatialRange-based intersection. Returning a SpatialRange didn't work, but returning SpatialRange* did. Could use some more testing. Note that the use of TopBit (compared to SciDB format) can make low-level htm-id output confusing or hard to compare with inputs to SpatialRange. --- include/HtmRangeMultiLevel.h | 184 +++++++++++++++-------------- include/SpatialRange.h | 7 +- src/HstmRange.C | 5 +- src/HtmRangeMultiLevel.cpp | 114 ++++++++++++++++-- src/SpatialIndex.cpp | 5 +- src/SpatialRange.cpp | 36 +++++- tests/CUTE/SpatialRange_test.cpp | 193 ++++++++++++++++++++++++++----- 7 files changed, 410 insertions(+), 134 deletions(-) diff --git a/include/HtmRangeMultiLevel.h b/include/HtmRangeMultiLevel.h index b3f0072..18d1cd4 100644 --- a/include/HtmRangeMultiLevel.h +++ b/include/HtmRangeMultiLevel.h @@ -22,11 +22,11 @@ namespace HtmRangeMultiLevel_NameSpace { enum InclusionType { - InclOutside = 0, /* */ - InclInside, /* */ - InclLo, /* number is on low end of an interval */ - InclHi, /* number is on high end of an interval */ - InclAdjacentXXX + InclOutside = 0, /* */ + InclInside, /* */ + InclLo, /* number is on low end of an interval */ + InclHi, /* number is on high end of an interval */ + InclAdjacentXXX }; struct TInsideResult { @@ -47,89 +47,97 @@ struct TInsideResult { */ class LINKAGE HtmRangeMultiLevel { - public: - static int HIGHS; - static int LOWS; - static int BOTH; - int getNext(Key &lo, Key &hi); - int getNext(Key *lo, Key *hi); - KeyPair getNext(); - int getNext(KeyPair &kp); - - void setSymbolic(bool flag); - - void addRange(const Key lohi); - void addRange(const Key lo, const Key hi); - void addRange(HtmRangeMultiLevel *range); - void mergeRange(const Key lo, const Key hi); - void defrag(); - void defrag(Key gap); - void CompressionPass(); - void purge(); - int isIn(Key key); - int isIn(Key lo, Key hi); - int isIn(HtmRangeMultiLevel & otherRange); - Key bestgap(Key desiredSize); - - HtmRangeMultiLevel getSpan(); - void parse(std::string rangeString); - - bool equalp(HtmRangeMultiLevel *other); - - int stats(int desiredSize); - int nranges(); - int nindexes_in_ranges(); - - int getLosLength() { - return my_los->getLength(); - } - int getHisLength() { - return my_his->getLength(); - } - void reset(); - - void print(std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... - void print(int what, std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... - - int compare(const HtmRangeMultiLevel & other) const; - - /** Hold the encoding scheme for the symbolic representation. - */ - EmbeddedLevelNameEncoding *encoding; - /** Change or set symbolic encoding. - */ - void setEncoding(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; } - /** Return the encoding for translations, etc. - */ - EmbeddedLevelNameEncoding *getEncoding() { return encoding; } - - // Moved here per ajmendez. - // TODO MLR Why two skip lists for ranges? Why not one skip list? Or an interval arithmetic package? - SkipList *my_los; - SkipList *my_his; - bool symbolicOutput; - - HtmRangeMultiLevel(); - HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding); - ~HtmRangeMultiLevel(){ - purge(); - delete encoding; - delete my_los; - delete my_his; - }; - - HtmRangeMultiLevel *HtmRangeMultiLevelAtLevelFromIntersection(HtmRangeMultiLevel *range2, int htmIdLevel=-1); - int contains(Key a, Key b); - - friend LINKAGE ostream& operator<<(ostream& os, const HtmRangeMultiLevel& range); - - protected: - TInsideResult tinside(const Key mid) const; - // const char buffer[256]; -// bool symbolicOutput; -// private: -// SkipList *my_los; -// SkipList *my_his; +public: + + // For the stored indices + bool embeddedLevel = true; // Set to false for deprecated legacy behavior. + + // For the range data type + static int HIGHS; + static int LOWS; + static int BOTH; + int getNext(Key &lo, Key &hi); + int getNext(Key *lo, Key *hi); + KeyPair getNext(); + int getNext(KeyPair &kp); + + void setSymbolic(bool flag); + + void addRange(const Key lohi); + void addRange(const Key lo, const Key hi); + void addRange(HtmRangeMultiLevel *range); + void mergeRange(const Key lo, const Key hi); + void defrag(); + void defrag(Key gap); + void CompressionPass(); + void purge(); + int isIn(Key key); + int isIn(Key lo, Key hi); + int isIn(HtmRangeMultiLevel & otherRange); + Key bestgap(Key desiredSize); + + HtmRangeMultiLevel getSpan(); + void parse(std::string rangeString); + + bool equalp(HtmRangeMultiLevel *other); + + int stats(int desiredSize); + int nranges(); + int nindexes_in_ranges(); + + int getLosLength() { + return my_los->getLength(); + } + int getHisLength() { + return my_his->getLength(); + } + void reset(); + + void print(std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... + void print(int what, std::ostream& os, bool symbolic = false); // FIX THIS, so caller does not set symbolic here.... + + int compare(const HtmRangeMultiLevel & other) const; + + /** Hold the encoding scheme for the symbolic representation. + */ + EmbeddedLevelNameEncoding *encoding; + /** Change or set symbolic encoding. + */ + void setEncoding(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; } + /** Return the encoding for translations, etc. + */ + EmbeddedLevelNameEncoding *getEncoding() { return encoding; } + + // Moved here per ajmendez. + // TODO MLR Why two skip lists for ranges? Why not one skip list? Or an interval arithmetic package? + SkipList *my_los; + SkipList *my_his; + bool symbolicOutput; + + HtmRangeMultiLevel(); + HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding); + ~HtmRangeMultiLevel(){ + purge(); + delete encoding; + delete my_los; + delete my_his; + }; + + HtmRangeMultiLevel *HtmRangeMultiLevelAtLevelFromIntersection(HtmRangeMultiLevel *range2, int htmIdLevel=-1); + HtmRangeMultiLevel *RangeFromIntersection(HtmRangeMultiLevel *range2, int force_htmIdLevel=-1); + + int contains(Key a, Key b); + // TODO KeyPair contains(Key a, Key b); + + friend LINKAGE ostream& operator<<(ostream& os, const HtmRangeMultiLevel& range); + +protected: + TInsideResult tinside(const Key mid) const; + // const char buffer[256]; + // bool symbolicOutput; + // private: + // SkipList *my_los; + // SkipList *my_his; }; #define SKIP_PROB (0.5) diff --git a/include/SpatialRange.h b/include/SpatialRange.h index f612813..d90f6c6 100644 --- a/include/SpatialRange.h +++ b/include/SpatialRange.h @@ -46,10 +46,11 @@ class SpatialRange { }; -inline SpatialRange operator& ( const SpatialRange& a, const SpatialRange& b) { - return SpatialRange(new HstmRange(a.range->range->HtmRangeMultiLevelAtLevelFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. -} +SpatialRange* sr_intersect(const SpatialRange&a, const SpatialRange& b); +inline SpatialRange* operator& ( const SpatialRange& a, const SpatialRange& b) { + return new SpatialRange(new HstmRange(a.range->range->RangeFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. +} void SpatialRange_test(); diff --git a/src/HstmRange.C b/src/HstmRange.C index 7c10654..4f82c9b 100644 --- a/src/HstmRange.C +++ b/src/HstmRange.C @@ -28,7 +28,7 @@ HstmRange::~HstmRange() { * @param b_ a Key (int64) */ void HstmRange::addRange(Key a_, Key b_) { - Key a = leftJustifiedEncoding.maskOffLevelBit(a_); + Key a = leftJustifiedEncoding.maskOffLevelBit(a_); // The level bit is a bit at the top used as a validity check. A deprecated feature and not needed for a left-justified index value. int aLevel = leftJustifiedEncoding.levelById(a_); Key b = leftJustifiedEncoding.maskOffLevelAndLevelBit(b_); int bLevel = leftJustifiedEncoding.levelById(b_); @@ -83,7 +83,8 @@ void HstmRange::addRange(Key a_, Key b_) { */ // #define DIAG -#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << " " << s << endl << flush; +#undef DIAG +#define IDOUT(p,m,s) p << m << " " << setw(16) << setfill('0') << hex << s << dec << " " << s << endl << flush; #ifdef DIAG IDOUT(cout,"hr::ar a_: ",a_) IDOUT(cout,"hr::ar a : ",a) diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index 933a3ab..b1f3220 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -12,7 +12,7 @@ #include // setw() #include // various *RepresentationString elements #include -#include // levelOfId +#include // levelOfId // Do wwe really need this? // TODO If levelOfId not needed, remove. #ifdef _WIN32 #include @@ -27,6 +27,15 @@ using namespace std; using namespace HtmRangeMultiLevel_NameSpace; +uint64 HRML_levelOfId(uint64 id,bool embeddedLevel,uint64 levelMask) { + if( embeddedLevel ) { + return id & levelMask; + } else { + // legacy behavior + return levelOfId(id); + } +} + /** * Translate an HtmRangeMultiLevel to one at a greater level. If the desired level is * less that the level implicit in the input range (lo & hi), then just return @@ -41,11 +50,11 @@ using namespace HtmRangeMultiLevel_NameSpace; */ KeyPair HtmRangeMultiLevelAtLevelFromHtmRangeMultiLevel(int htmIdLevel, Key lo, Key hi) { // htmIdLevel is used to set maxlevel in the index. aka olevel. - int levelLo = levelOfId(lo); + int levelLo = HRML_levelOfId(lo,false,63); if(levelLonranges()<=0)||(range2->nranges()<=0)) return 0; @@ -78,14 +87,14 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::HtmRangeMultiLevelAtLevelFromIntersectio if (!indexp1) return 0; if(htmIdLevel<0) { - htmIdLevel = levelOfId(lo1); + htmIdLevel = HRML_levelOfId(lo1,false,63); } // cout << "indexp1: " << indexp1 << endl << flush; // cout << "l,lo,hi1: " << htmIdLevel << " " << lo1 << " " << hi1 << endl << flush; - // cout << "a" << flush; + cout << "a" << flush; do { - // cout << "b" << endl << flush; + cout << "b" << endl << flush; KeyPair testRange1 = HtmRangeMultiLevelAtLevelFromHtmRangeMultiLevel(htmIdLevel,lo1,hi1); range2->reset(); uint64 indexp2 = range2->getNext(lo2,hi2); @@ -119,8 +128,90 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::HtmRangeMultiLevelAtLevelFromIntersectio return resultRange; } +KeyPair HRML_AtLevelFromMultiLevel(uint64 htmIdLevel, Key lo, Key hi, uint64 levelMask) { + uint64 levelLo = lo & levelMask; + const uint64 one = 1; + if( levelLo < htmIdLevel) { + lo = (lo & ~levelMask) | htmIdLevel; + // Ignore weird cases where level(lo) != level(hi) + // Make hi into a terminator for htmIdLevel. + hi = hi | ( (one << (6+54-2*htmIdLevel)) - 1 ); + } + KeyPair levelAdaptedRange; + levelAdaptedRange.lo = lo; + levelAdaptedRange.hi = hi; + return levelAdaptedRange; +} + +/** + * Find intersection of two ranges and return as a range. + * + * TODO Replace double-nested-loop with SkipLists search or find functionality. + * + */ +HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel *range2, int force_htmIdLevel) { + HtmRangeMultiLevel *range1 = this; // Just an alias + if((!range1)||(!range2)) return 0; + if((range1->nranges()<=0)||(range2->nranges()<=0)) return 0; + Key lo1,hi1,lo2,hi2; + range1->reset(); + uint64 indexp1 = range1->getNext(lo1,hi1); + if(!indexp1) return 0; + if(force_htmIdLevel<0) { + force_htmIdLevel = lo1 & this->encoding->levelMask; // TODO Establish 31 or 63? + } + HtmRangeMultiLevel *resultRange = new HtmRangeMultiLevel(); resultRange->purge(); + do { + KeyPair testRange1 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo1,hi1,this->encoding->levelMask); + range2->reset(); // Sigh. Reset and loop from the beginning. TODO Avoid restarting loop. There must be a faster way. + uint64 indexp2 = range2->getNext(lo2,hi2); // TODO Implement a find or search for inserting. + if(!indexp2) return 0; + bool intersects = false; +#define FMTX(x) setw(16) << setfill('0') << hex << x << dec + do { + KeyPair testRange2 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo2,hi2,this->encoding->levelMask); + intersects = testRange2.lo <= testRange1.hi + && testRange2.hi >= testRange1.lo; +#ifdef DIAG + cout << "lh1,lh2: " + << FMTX(lo1) << " " << FMTX(hi1) << ", " + << FMTX(lo2) << " " << FMTX(hi2) << ", " + << intersects << flush; +#endif + if(intersects){ + Key lo_ = max(testRange1.lo,testRange2.lo); + Key hi_ = min(testRange1.hi,testRange2.hi); + resultRange->addRange(lo_,hi_); +#ifdef DIAG + cout << ", added " + << FMTX(lo_) << " " + << FMTX(hi_) << flush; +#endif + } +#ifdef DIAG + cout << "." << endl << flush; +#endif + } while (range2->getNext(lo2,hi2)); + } while (range1->getNext(lo1,hi1)); // TODO Can we replace getNext with some sort of find or search. +#undef FMTX + // cout << "d" << flush; + // cout << "d nr " << resultRange->nranges() << endl << flush; + // cout << "d rr " << hex << resultRange << dec << endl << flush; + if(resultRange->nranges()>0)resultRange->defrag(); + // cout << "e" << flush; + return resultRange; +} + + // Note the default use of EmbeddedLevelNameEncoding. -HtmRangeMultiLevel::HtmRangeMultiLevel() : HtmRangeMultiLevel(new EmbeddedLevelNameEncoding()) {} +// HtmRangeMultiLevel::HtmRangeMultiLevel() : HtmRangeMultiLevel(new EmbeddedLevelNameEncoding()) {} +HtmRangeMultiLevel::HtmRangeMultiLevel() { + encoding = new EmbeddedLevelNameEncoding(); + my_los = new SkipList(SKIP_PROB); + // cout << "hrml my_los " << hex << my_los << dec << endl << flush; + my_his = new SkipList(SKIP_PROB); + symbolicOutput = false; +} HtmRangeMultiLevel::HtmRangeMultiLevel(EmbeddedLevelNameEncoding *encoding) { this->encoding = encoding; @@ -1076,14 +1167,17 @@ void HtmRangeMultiLevel::reset() /// The number of ranges. int HtmRangeMultiLevel::nranges() { -// cout << "z000" << endl << flush; + // cout << "z000" << endl << flush; Key lo; // Key hi; int n_ranges; n_ranges = 0; + // cout << "z001" << endl << flush; + // cout << "z001 my_los " << hex << my_los << dec << endl << flush; my_los->reset(); + // cout << "z002" << endl << flush; my_his->reset(); -// cout << "z010" << endl << flush; + // cout << "z010" << endl << flush; // This is a problem when lo can be zero. Is it? // getkey returns -1 if nothing is found, maybe fix the following using >= 0? Worry about id 0. Should be okay this low in the code. MLR 2019-0327 diff --git a/src/SpatialIndex.cpp b/src/SpatialIndex.cpp index 8d68780..f93db8f 100644 --- a/src/SpatialIndex.cpp +++ b/src/SpatialIndex.cpp @@ -1457,14 +1457,15 @@ int depthOfId(uint64 htmId) { int levelOfId(uint64 htmId) { int i; uint32 size; + // cout << "loid: " << setw(16) << setfill('0') << hex << htmId << dec << endl << flush; // determine index of first set bit for(i = 0; i < IDSIZE; i+=2) { if ( (htmId << i) & IDHIGHBIT ) break; if ( (htmId << i) & IDHIGHBIT2 ) // invalid id - throw SpatialFailure("SpatialIndex:nameById: invalid ID"); + throw SpatialFailure("SpatialIndex:nameById: invalid ID 1"); } if(htmId == 0) - throw SpatialFailure("SpatialIndex:nameById: invalid ID"); + throw SpatialFailure("SpatialIndex:nameById: invalid ID 2"); size=(IDSIZE-i) >> 1; return size-2; } diff --git a/src/SpatialRange.cpp b/src/SpatialRange.cpp index ebd5220..6d00725 100644 --- a/src/SpatialRange.cpp +++ b/src/SpatialRange.cpp @@ -62,7 +62,7 @@ void SpatialRange::addSpatialIntervals(STARE_SpatialIntervals intervals) { } } - +// #define DIAG int SpatialRange::getNextSpatialInterval(STARE_SpatialIntervals &interval) { KeyPair kp(-1,-2); int istat = this->getNext(kp); @@ -99,9 +99,39 @@ int SpatialRange::getNextSpatialInterval(STARE_SpatialIntervals &interval) { STARE_SpatialIntervals SpatialRange::toSpatialIntervals() { STARE_SpatialIntervals intervals; - this->range->reset(); - while( this->getNextSpatialInterval(intervals) > 0 ); + if(this->range) { + this->range->reset(); + while( this->getNextSpatialInterval(intervals) > 0 ); + } return intervals; } +/* + * Odd. The following does not seem to work if we just return the SpatialRange itself. Some of the pointers seem to be either corrupted or eliminated. + */ +SpatialRange *sr_intersect(const SpatialRange&a, const SpatialRange& b) { + HstmRange *range = new HstmRange(a.range->range->RangeFromIntersection(b.range->range)); // NOTE mlr Probably about the safest way to inst. SpatialRange. +// #define DIAG +#ifdef DIAG + KeyPair kp; range->reset(); range->getNext(kp); + cout << "sr_i range,r->r,nr " << range << " " << range->range << " " << range->range->nranges() << " : " + << setw(16) << setfill('0') << hex << kp.lo << " " + << setw(16) << setfill('0') << hex << kp.hi << " " + << dec + << endl << flush; + EmbeddedLevelNameEncoding leftJustified; + leftJustified.setId(kp.lo); cout << "kp.lo lj " << setw(16) << setfill('0') << hex << leftJustified.getSciDBLeftJustifiedFormat() << endl << flush; + leftJustified.setId(kp.hi); cout << "kp.hi lj " << setw(16) << setfill('0') << hex << leftJustified.getSciDBLeftJustifiedFormat() << endl << flush; + cout << " r-r-my_los " << hex << range->range->my_los << endl << flush; + cout << dec; +#endif + SpatialRange *sr = new SpatialRange(range); +#ifdef DIAG + cout << "sr-r " << hex << sr->range << endl << flush; + cout << "sr-r-r " << hex << sr->range->range << endl << flush; + cout << "sr-r-r-my_los " << hex << sr->range->range->my_los << endl << flush; +#endif + return sr; +} + diff --git a/tests/CUTE/SpatialRange_test.cpp b/tests/CUTE/SpatialRange_test.cpp index 80e1388..ae81a4b 100644 --- a/tests/CUTE/SpatialRange_test.cpp +++ b/tests/CUTE/SpatialRange_test.cpp @@ -12,67 +12,208 @@ void SpatialRange_test () { STARE index; - STARE_ArrayIndexSpatialValue sid = index.ValueFromLatLonDegrees( 30, 30, 8 ); - STARE_SpatialIntervals sids; sids.push_back(sid); + STARE_ArrayIndexSpatialValue siv = index.ValueFromLatLonDegrees( 30, 30, 8 ); + STARE_SpatialIntervals sis; sis.push_back(siv); - SpatialRange range(sids); - STARE_SpatialIntervals sids_out; + SpatialRange range(sis); + STARE_SpatialIntervals sis_out; - sids_out = range.toSpatialIntervals(); + sis_out = range.toSpatialIntervals(); // #define DIAG #ifdef DIAG - for(int i=0; i < sids.size(); ++i ) { - cout << i << " i,si: " << setw(16) << setfill('0') << hex << sids[i] << dec << endl << flush; + for(int i=0; i < sis.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; } - for(int i=0; i < sids_out.size(); ++i ) { - cout << i << " i,so: " << setw(16) << setfill('0') << hex << sids_out[i] << dec << endl << flush; + for(int i=0; i < sis_out.size(); ++i ) { + cout << i << " i,so: " << setw(16) << setfill('0') << hex << sis_out[i] << dec << endl << flush; } #endif - ASSERT_EQUAL(1,sids_out.size()); - ASSERT_EQUAL(sids[0],sids_out[0]); + ASSERT_EQUAL(1,sis_out.size()); + ASSERT_EQUAL(sis[0],sis_out[0]); // sids.push_back(index.ValueFromLatLonDegrees( 45, 45, 8)); - sid = index.ValueFromLatLonDegrees( 15, 15, 8 ); + siv = index.ValueFromLatLonDegrees( 15, 15, 8 ); EmbeddedLevelNameEncoding leftJustified; - leftJustified.setIdFromSciDBLeftJustifiedFormat(sid); - sids.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv); + sis.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); #ifdef DIAG - for(int i=0; i < sids.size(); ++i ) { - cout << i << " i,si: " << setw(16) << setfill('0') << hex << sids[i] << dec << endl << flush; + for(int i=0; i < sis.size(); ++i ) { + cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; } -// for(int i=0; i < sids_out.size(); ++i ) { -// cout << i << " i,so: " << setw(16) << setfill('0') << hex << sids_out[i] << dec << endl << flush; +// for(int i=0; i < sis_out.size(); ++i ) { +// cout << i << " i,so: " << setw(16) << setfill('0') << hex << sis_out[i] << dec << endl << flush; // } #endif // cout << 50 << endl << flush; range.purge(); - range.addSpatialIntervals(sids); + range.addSpatialIntervals(sis); // cout << 100 // << " nr: " << range.range->range->nranges() // << endl << flush; - sids_out = range.toSpatialIntervals(); + sis_out = range.toSpatialIntervals(); #ifdef DIAG - cout << 110 << " sids_out.s: " << sids_out.size() << endl << flush; - for(int i=0; irange) << dec << endl << flush; + cout << 1991 << " dR.r->r " << hex << (deltaRange->range->range) << dec << endl << flush; + cout << 1991 << " dR.r->r->my_los " << hex << (deltaRange->range->range->my_los) << dec << endl << flush; +#endif + if(!deltaRange->range->range) { cout << "Error deltaRange is null!" << endl << flush; } + // cout << 200 << " dR nR = " << deltaRange->range->range->nranges() << endl << flush; + STARE_SpatialIntervals deltaSis = deltaRange->toSpatialIntervals(); + // cout << 300 << endl << flush; + SISOUT("deltaSis",deltaSis); + // cout << 400 << endl << flush; + + ASSERT_EQUAL(sis1[0],deltaSis[0]); + ASSERT_EQUAL(sis0[1],deltaSis[1]); + +#undef SIVOUT +#undef SISOUT + + } + + + if(true) { + // TODO Fix the following to use the new 0/1 variables defined above. +// #define DIAG +#ifndef DIAG +#define SIVOUT(m,siv) +#define SISOUT(m,sis) +#else +#define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SISOUT(m,sis) { cout << endl << m << endl; for(int i=0; i < sis.size(); ++i ) { cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; }} +#endif + // #define DIAGOUT + // #define DIAGOUT(m) {cout << m << endl << flush;} + // Note level is 8. + uint64 one_mask_to_level, one_at_level; + leftJustified.increment_LevelToMaskDelta(8,one_mask_to_level, one_at_level); + + STARE_SpatialIntervals sis0; + SpatialRange range0; + STARE_ArrayIndexSpatialValue siv0 = 0x00000000000000008; + sis0.push_back(siv0); + siv0+=6*one_at_level; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); // TODO Rename this STARE, or provide a STARE interface for this. + sis0.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + range0.addSpatialIntervals(sis0); + SISOUT("sis0",sis0); + + STARE_SpatialIntervals sis1; + SpatialRange range1; + STARE_ArrayIndexSpatialValue siv1 = 0x00000000000000008; + siv1+=3*one_at_level; + siv1 = ( siv1 & ~leftJustified.levelMaskSciDB ) | 10; + sis1.push_back(siv1); + siv1+=6*one_at_level; + leftJustified.setIdFromSciDBLeftJustifiedFormat(siv1); + sis1.push_back(leftJustified.getSciDBTerminatorLeftJustifiedFormat()); + range1.addSpatialIntervals(sis1); + SISOUT("sis1",sis1); + + // cout << 100 << endl << flush; + SpatialRange *deltaRange; + try { + // cout << 110 << endl << flush; + deltaRange = range0 & range1; + // deltaRange = sr_intersect(range0, range1); + // cout << 120 << endl << flush; + } catch ( SpatialException e ) { + // cout << 130 << endl << flush; + cout << e.what() << endl << flush; + // cout << 140 << endl << flush; + } + +#ifdef DIAG + cout << 199 << endl << flush; + cout << 1991 << " dR.r " << hex << (deltaRange->range) << dec << endl << flush; + cout << 1991 << " dR.r->r " << hex << (deltaRange->range->range) << dec << endl << flush; + cout << 1991 << " dR.r->r->my_los " << hex << (deltaRange->range->range->my_los) << dec << endl << flush; +#endif + if(!deltaRange->range->range) { cout << "Error deltaRange is null!" << endl << flush; } + // cout << 200 << " dR nR = " << deltaRange->range->range->nranges() << endl << flush; + STARE_SpatialIntervals deltaSis = deltaRange->toSpatialIntervals(); + // cout << 300 << endl << flush; + SISOUT("deltaSis",deltaSis); + // cout << 400 << endl << flush; + + ASSERT_EQUAL(sis1[0],deltaSis[0]); + ASSERT_EQUAL(sis0[1],deltaSis[1]); + +#undef SIVOUT +#undef SISOUT +#undef DIAG + } +// TODO Add intersection tests. } From 71a156747cba549cb576db41e641d1a02fbcc2e2 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Thu, 12 Sep 2019 11:53:06 -0400 Subject: [PATCH 05/31] Updating notes. --- NOTES | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NOTES b/NOTES index 36aaabc..9fe37b8 100644 --- a/NOTES +++ b/NOTES @@ -1,4 +1,9 @@ +* 2019-09-X Version 0.8.0 + +Added SpatialRange to use SkipList to generate intersections. Note the intersect +doesn't use find or search, so it's not clear we're getting the benefit of the SkipLists. + * 2019-09-05 Version 0.7.0 Added Griessenbaum's new apps and python interface. From f7e9a964dbe57687d7d19d1ed2fd4ceac5475c6f Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Sun, 15 Sep 2019 01:08:11 -0400 Subject: [PATCH 06/31] More progress adding SpatialRange support. Improving understanding of RangeConvex etc. See VizHTM. --- include/STARE.h.in | 4 +++- include/SpatialRange.h | 2 +- src/HtmRangeMultiLevel.cpp | 2 ++ src/RangeConvex.cpp | 4 ++-- src/STARE.C | 8 ++++++++ 5 files changed, 16 insertions(+), 4 deletions(-) diff --git a/include/STARE.h.in b/include/STARE.h.in index 1cf72f6..1368dd0 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -143,7 +143,7 @@ STARE_ArrayIndexSpatialValue sTerminator(STARE_ArrayIndexSpatialValue spatialSta bool terminatorp(STARE_ArrayIndexSpatialValue spatialStareId); /// Check if the index value is a terminator. -HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); +// Deprecated... cf. SpatialRange - HstmRange SpatialRangeFromSpatialIntervals(STARE_SpatialIntervals intervals); void STARE_ArrayIndexSpatialValues_insert(STARE_ArrayIndexSpatialValues& v, STARE_ArrayIndexSpatialValue siv); @@ -155,6 +155,8 @@ STARE_ArrayIndexSpatialValue shiftSpatialIdAtLevel( int resolution, int shiftAmount ); + +STARE_SpatialIntervals spatialIntervalFromHtmIDKeyPair(KeyPair kp); uint64 spatialLevelMask(); diff --git a/include/SpatialRange.h b/include/SpatialRange.h index d90f6c6..d367113 100644 --- a/include/SpatialRange.h +++ b/include/SpatialRange.h @@ -46,7 +46,7 @@ class SpatialRange { }; -SpatialRange* sr_intersect(const SpatialRange&a, const SpatialRange& b); +SpatialRange* sr_intersect(const SpatialRange& a, const SpatialRange& b); inline SpatialRange* operator& ( const SpatialRange& a, const SpatialRange& b) { return new SpatialRange(new HstmRange(a.range->range->RangeFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index b1f3220..f01e92d 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -168,6 +168,8 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel if(!indexp2) return 0; bool intersects = false; #define FMTX(x) setw(16) << setfill('0') << hex << x << dec + // Search forward until we find an intersection. Once an intersection is found, + // figure out what the intersection is and add it to the result range. do { KeyPair testRange2 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo2,hi2,this->encoding->levelMask); intersects = testRange2.lo <= testRange1.hi diff --git a/src/RangeConvex.cpp b/src/RangeConvex.cpp index b9915e8..db5f8a7 100644 --- a/src/RangeConvex.cpp +++ b/src/RangeConvex.cpp @@ -920,7 +920,7 @@ SpatialMarkup RangeConvex::testTrixel(uint64 nodeIndex) { // cout << endl << flush; // cout << " check children & partials " << endl << flush; - // NEW NEW algorithm Disabled when enablenew is 0 + // NEW < algorithm Disabled when enablenew is 0 { childID = indexNode->childID_[0]; if (childID != 0) { @@ -951,7 +951,7 @@ SpatialMarkup RangeConvex::testTrixel(uint64 nodeIndex) { If partial, then we look ahead to see how many children are rejected. But ah, next iteration could benefit from having computed this already. - If two chidlren are rejected, then we stop + If two children are rejected, then we stop If one or 0 nodes are rejected, then we */ // cout << " mark: " << mark << endl << flush; diff --git a/src/STARE.C b/src/STARE.C index 22f4d35..5f12e06 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -731,3 +731,11 @@ STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, // cout << dec << 160 << " " << i << endl << flush; return expanded_values; } + +STARE_SpatialIntervals spatialIntervalFromHtmIDKeyPair(KeyPair kp) { + STARE_SpatialIntervals interval; + interval.push_back(EmbeddedLevelNameEncoding(BitShiftNameEncoding(kp.lo).leftJustifiedId()).getSciDBLeftJustifiedFormat()); + interval.push_back(EmbeddedLevelNameEncoding(BitShiftNameEncoding(kp.hi).leftJustifiedId()).getSciDBTerminatorLeftJustifiedFormat()); + return interval; +} + From 05e20714b9fcd8717c07e49dfcb6bf97a3e3cacc Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Mon, 16 Sep 2019 08:09:47 -0400 Subject: [PATCH 07/31] ConvexHull, intersection improvements. --- include/STARE.h.in | 1 + include/SpatialConstraint.h | 4 ++++ include/SpatialInterface.h | 6 +++++- src/STARE.C | 26 +++++++++++++++-------- src/SpatialConstraint.cpp | 8 ++++++++ src/SpatialInterface.cpp | 41 ++++++++++++++++++++++--------------- 6 files changed, 59 insertions(+), 27 deletions(-) diff --git a/include/STARE.h.in b/include/STARE.h.in index 1368dd0..e66fd20 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -60,6 +60,7 @@ public: // Spatial array index functions. [Maybe change the name StareId to spatialStareId in the following...] STARE_ArrayIndexSpatialValue ValueFromLatLonDegrees(float64 latDegrees, float64 lonDegrees, int resolutionLevel = STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); LatLonDegrees64 LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spatialStareId); + SpatialVector SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId); uint32 ResolutionLevelFromValue(STARE_ArrayIndexSpatialValue spatialStareId); Triangle TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel = -1); diff --git a/include/SpatialConstraint.h b/include/SpatialConstraint.h index de6e908..109fe9e 100644 --- a/include/SpatialConstraint.h +++ b/include/SpatialConstraint.h @@ -81,6 +81,10 @@ class LINKAGE SpatialConstraint : public SpatialSign { /// Initialization constructor SpatialConstraint(SpatialVector, float64); + /// Construct a zero constraint through two points + // Makes a zero-angle constraint directed towards v0^v1. + SpatialConstraint(SpatialVector v0, SpatialVector v1); + /// Copy constructor SpatialConstraint(const SpatialConstraint &); diff --git a/include/SpatialInterface.h b/include/SpatialInterface.h index c0bd648..7465dab 100644 --- a/include/SpatialInterface.h +++ b/include/SpatialInterface.h @@ -72,6 +72,8 @@ class LINKAGE htmInterface { saveDepth parameter can be specified to keep the given amount of levels in memory. This can also be altered by changeDepth. */ + htmInterface(const SpatialIndex *index); + htmInterface(size_t searchlevel = 5, size_t buildevel = 5, SpatialRotation rot = rot_identity ); // [ed:gyuri:saveDepth was 2] /** Destructor. */ @@ -187,7 +189,8 @@ class LINKAGE htmInterface { * @return */ const HTMRangeValueVector & convexHull( LatLonDegrees64ValueVector latlon, - size_t steps = -1); + size_t steps = -1, + bool interiorp = false); /** Request all triangles in the convex hull of a given set of points. @@ -302,6 +305,7 @@ class LINKAGE htmInterface { // void setPolyCorner(SpatialVector &v); // this routine does the work for all convexHull calls + bool hull_interiorp_ = false; const HTMRangeValueVector & doHull(); }; diff --git a/src/STARE.C b/src/STARE.C index 5f12e06..7a002c2 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -181,6 +181,14 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati } +SpatialVector STARE::SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId) { + uint64 htmID = htmIDFromValue(spatialStareId,STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); // Max resolution + SpatialVector v; + /// This returns the center of the triangle (at index.search_level). Need to extract the position information. + sIndex.pointByHtmId(v, htmID); + return v; +} + Triangle STARE::TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel) { // Users are going to expect the default resolution level to be that embedded in the sStareId. uint64 htmID = -1; @@ -433,30 +441,30 @@ STARE_SpatialIntervals STARE::ConvexHull(LatLonDegrees64ValueVector points,int f STARE_SpatialIntervals cover; int hullSteps = points.size(); htmInterface *htm; - cout << dec << 1000 << " hullSteps: " << hullSteps << endl << flush; + // cout << dec << 1000 << " hullSteps: " << hullSteps << endl << flush; if( force_resolution_level > -1 ) { - cout << dec << 1100 << endl << flush; + // cout << dec << 1100 << endl << flush; // htm = htmInterface(&index_); htm = new htmInterface( this->getIndex(force_resolution_level).getMaxlevel(), this->getIndex(force_resolution_level).getBuildLevel(), this->getIndex(force_resolution_level).getRotation()); - cout << dec << 1101 << endl << flush; + // cout << dec << 1101 << endl << flush; } else { - cout << dec << 1200 << endl << flush; + // cout << dec << 1200 << endl << flush; // htm = htmInterface(&index_); htm = new htmInterface( this->getIndex(8).getMaxlevel(), this->getIndex(8).getBuildLevel(), this->getIndex(8).getRotation()); - cout << dec << 1201 << endl << flush; + // cout << dec << 1201 << endl << flush; } - cout << dec << "a2000" << endl << flush; + // cout << dec << "a2000" << endl << flush; - HTMRangeValueVector htmRangeVector = htm->convexHull(points,hullSteps); + HTMRangeValueVector htmRangeVector = htm->convexHull(points,hullSteps,true); // TODO FIX interiorp = false is broken - cout << dec << "a3000 hrv.size: " << htmRangeVector.size() << endl << flush; + // cout << dec << "a3000 hrv.size: " << htmRangeVector.size() << endl << flush; for(int i=0; i < htmRangeVector.size(); ++i) { uint64 lo = ValueFromHtmID(htmRangeVector[i].lo); // TODO Should this be a function? @@ -468,7 +476,7 @@ STARE_SpatialIntervals STARE::ConvexHull(LatLonDegrees64ValueVector points,int f } } - cout << dec << "a4000" << endl << flush; + // cout << dec << "a4000" << endl << flush; delete htm; // TODO Hopefully this will not also delete the index we passed in. return cover; diff --git a/src/SpatialConstraint.cpp b/src/SpatialConstraint.cpp index 7adcfd1..e9da5fd 100644 --- a/src/SpatialConstraint.cpp +++ b/src/SpatialConstraint.cpp @@ -60,6 +60,14 @@ SpatialConstraint::SpatialConstraint(SpatialVector a, float64 d) : if(d_ >= gEpsilon) sign_ = pOS; } +/////////////CONSTRUCTOR Zero Constraint through two points////////////////////////////////// +// +SpatialConstraint::SpatialConstraint(SpatialVector v0, SpatialVector v1) +{ + a_ = v0 ^ v1; a_.normalize(); d_ = 0; s_ = acos(d_); +} + + /////////////COPY CONSTRUCTOR///////////////////////////// // SpatialConstraint::SpatialConstraint(const SpatialConstraint & old) : diff --git a/src/SpatialInterface.cpp b/src/SpatialInterface.cpp index b9f9ccb..53bf625 100644 --- a/src/SpatialInterface.cpp +++ b/src/SpatialInterface.cpp @@ -40,6 +40,9 @@ extern long long atoll (const char *str); //============================================================== ///////////CONSTRUCTOR/////////////////////// +htmInterface::htmInterface(const SpatialIndex *index) { + index_ = new SpatialIndex(index->maxlevel_, index->buildlevel_, index->rot_); // TODO delete bait? maxlevel is searchlevel, no? +} htmInterface::htmInterface(size_t searchlevel, size_t buildlevel, SpatialRotation rot) : t_(NULL) { index_ = new SpatialIndex(searchlevel, buildlevel, rot); } @@ -387,24 +390,25 @@ htmInterface::convexHull( } const HTMRangeValueVector & -htmInterface::convexHull( LatLonDegrees64ValueVector latlon, size_t steps ) { +htmInterface::convexHull( LatLonDegrees64ValueVector latlon, size_t steps, bool interiorp ) { + hull_interiorp_ = interiorp; polyCorners_.clear(); - cout << " ch " << 2000 << " latlon-size=" << latlon.size() << flush ; - cout << endl; + // cout << " ch " << 2000 << " latlon-size=" << latlon.size() << flush ; + // cout << endl; if (steps == (uint64) -1) { steps = latlon.size(); } else { steps = min(steps,latlon.size()); } for(size_t i = 0; i < steps; i++) { - cout << " " << i << flush; - cout << " ( " << latlon[i].lat << " " << latlon[i].lon << ")"; + // cout << " " << i << flush; + // cout << " ( " << latlon[i].lat << " " << latlon[i].lon << ")"; float64 *x = xyzFromLatLonDegrees(latlon[i].lat,latlon[i].lon); SpatialVector v(x[0],x[1],x[2]); setPolyCorner(v); - cout << endl << flush; + // cout << endl << flush; } - cout << endl << flush << 2100 << endl << flush; + // cout << endl << flush << 2100 << endl << flush; return doHull(); } @@ -478,7 +482,7 @@ htmInterface::doHull() { // dom.convexes_[0].boundingCircle_.write(cout); // dom.write(cout); //%%%%%%%%%%%%%%%%%%%%%%%%%%% - cout << 3999 << endl << flush; + // cout << 3999 << endl << flush; return domain(dom); } @@ -698,31 +702,34 @@ const HTMRangeValueVector & htmInterface::domain( SpatialDomain & domain ) { HtmRange htmRange; - cout << 4000 << endl << flush; +// cout << 4000 << endl << flush; +// cout << 4001 << " hull_interiorp_ " << hull_interiorp_ << endl << flush; Key gapsize; // SpatialIndex idx(20, 5); // domain.setOlevel(20); - domain.intersect(index_, &htmRange, false); + // domain.intersect(index_, &htmRange, false); + // domain.intersect(index_, &htmRange, true); + domain.intersect(index_, &htmRange, hull_interiorp_); - cout << 4100 << endl << flush; +// cout << 4100 << endl << flush; // gapsize = htmRange.bestgap(MAX_RANGES); // htmRange.defrag(gapsize); - cout << 4200 << endl << flush; -// htmRange.defrag(); +// cout << 4200 << endl << flush; + htmRange.defrag(); // DONT FORGET to: get best gap and defrag htmRange.defrag(bestgap); - cout << 4300 << endl << flush; +// cout << 4300 << endl << flush; // Construct the valuevector... fillValueVec( htmRange, range_); - cout << 4997 << endl << flush; +// cout << 4997 << endl << flush; htmRange.reset(); - cout << 4998 << endl << flush; +// cout << 4998 << endl << flush; htmRange.purge(); - cout << 4999 << endl << flush; +// cout << 4999 << endl << flush; return range_; } From 26cdc17bafacdec4ce7331e3ea89cdafc38f108d Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Mon, 16 Sep 2019 11:06:13 -0400 Subject: [PATCH 08/31] Improved intersection code to make use of SkipList. Avoids iteration over pairs. Also, assumes segments in a range are accessed in order. --- include/SkipListElement.h | 3 +- src/HtmRangeMultiLevel.cpp | 69 +++++++++++++++++++++++++++----------- 2 files changed, 51 insertions(+), 21 deletions(-) diff --git a/include/SkipListElement.h b/include/SkipListElement.h index 062bec2..a072100 100644 --- a/include/SkipListElement.h +++ b/include/SkipListElement.h @@ -36,7 +36,8 @@ typedef int64 Key; // key type TODO Why not uint64? // typedef uint64 Key; // key type TODO Why not uint64? -typedef int Value; // value type +// typedef int Value; // value type +typedef int64 Value; // mlr... // class ostream; class SkipListElement; diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index f01e92d..d975068 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -161,19 +161,29 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel force_htmIdLevel = lo1 & this->encoding->levelMask; // TODO Establish 31 or 63? } HtmRangeMultiLevel *resultRange = new HtmRangeMultiLevel(); resultRange->purge(); + do { KeyPair testRange1 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo1,hi1,this->encoding->levelMask); range2->reset(); // Sigh. Reset and loop from the beginning. TODO Avoid restarting loop. There must be a faster way. + /* Try to skip past by using SkipList functions. */ + // TODO do something like range2->findMAX... + Key loKey = range2->my_los->findMAX(testRange1.lo); + Value hiKey = range2->my_los->search(loKey,true); + Value vhi = range2->my_his->search(hiKey,true); + /**/ uint64 indexp2 = range2->getNext(lo2,hi2); // TODO Implement a find or search for inserting. - if(!indexp2) return 0; - bool intersects = false; + bool intersects = false, past_chance; #define FMTX(x) setw(16) << setfill('0') << hex << x << dec // Search forward until we find an intersection. Once an intersection is found, // figure out what the intersection is and add it to the result range. + // int kount=0; + if(indexp2) do { + // ++kount; KeyPair testRange2 = HRML_AtLevelFromMultiLevel(force_htmIdLevel,lo2,hi2,this->encoding->levelMask); intersects = testRange2.lo <= testRange1.hi && testRange2.hi >= testRange1.lo; +// #define DIAG #ifdef DIAG cout << "lh1,lh2: " << FMTX(lo1) << " " << FMTX(hi1) << ", " @@ -193,7 +203,10 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel #ifdef DIAG cout << "." << endl << flush; #endif - } while (range2->getNext(lo2,hi2)); +#undef DIAG + past_chance = (uint64) testRange2.lo > (uint64) testRange1.hi; + } while (range2->getNext(lo2,hi2) && !past_chance); + // cout << "kount = " << kount << endl << flush; } while (range1->getNext(lo1,hi1)); // TODO Can we replace getNext with some sort of find or search. #undef FMTX // cout << "d" << flush; @@ -542,7 +555,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // Add the first one. if( my_los->myHeader->getElement(0) == NIL ) { - my_los->insert(lo,100); + // my_los->insert(lo,100); + my_los->insert(lo,hi); my_his->insert(hi,100); // cout << 8002 << " First one inserted. " << endl << flush; return; @@ -601,7 +615,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // Don't know what's above h. Iterate. } else if( hi1 < l ) { // Case 1. A is below B. Just add - my_los->insert(lo1,10001); + // my_los->insert(lo1,10001); + my_los->insert(lo1,hi1); my_his->insert(hi1,10001); done = true; } else if( (lo1 < l) && ( (l <= hi1) && (hi1 <= h) ) ) { @@ -630,14 +645,16 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // At the same level, merge the two. my_los->freeRange(l_m,h_p); // Freeing up to h_p is okay because h_p==h is still part of current interval. my_his->freeRange(l_m,h_p); - my_los->insert(l_m,100021); + // my_los->insert(l_m,100021); + my_los->insert(l_m,h_p); my_his->insert(h_p,100021); } else { // The lower part overlaps an empty part. Just add. // ??? Don't need a freeRange ??? Okay... if(level > l_level) { // cout << "8000-1031" << endl << flush; // If the new interval's level is greater, just skip in the current, add before. - my_los->insert(l_m,100022); + // my_los->insert(l_m,100022); + my_los->insert(l_m,h_m); my_his->insert(h_m,100022); } else if(true) { // Case 2.3 level < l_level -- new interval wins // TODO WORRY -- What about collisions? If we have a collision, will we simply put in the value back in? @@ -654,10 +671,12 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) // } my_los->freeRange(l_m,h_p); my_his->freeRange(l_m,h_p); - my_los->insert(l_m,100023); // lo1 // this changes the level to level + // my_los->insert(l_m,100023); // lo1 // this changes the level to level + my_los->insert(l_m,h_0); // lo1 // this changes the level to level my_his->insert(h_0,100023); // hi1 if(h_0 < h_p) { - my_los->insert(l_p,100023); + // my_los->insert(l_p,100023); + my_los->insert(l_p,h_p); my_his->insert(h_p,100023); } my_los->reset(); my_his->reset(); @@ -690,7 +709,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) cout << "HtmRangeMultiLevel::mergeRange::ERROR!!! SUCC(H) > HI1, I.E. THEY'RE EQUIVALENT." << endl << flush; } if( l_m < h_m ) { // If lo1 and l are equivalent, current one wins, and we ignore the new one. - my_los->insert(l_m,1000421); + // my_los->insert(l_m,1000421); + my_los->insert(l_m,h_m); my_his->insert(h_m,1000421); // my_los->reset(); my_his->reset(); // } else { @@ -726,13 +746,16 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) my_los->freeRange(l_m,h_p); my_his->freeRange(l_m,h_p); if( l_m < h_m ) { // If l_m and l_0 are "equivalent", so the current interval wins and we ignore interval_m. - my_los->insert(l_m,10004); + // my_los->insert(l_m,10004); + my_los->insert(l_m,h_m); my_his->insert(h_m,10004); } - my_los->insert(l_0,10004); + // my_los->insert(l_0,10004); + my_los->insert(l_0,h_0); my_his->insert(h_0,10004); // TODO Subtle bug? Need to verify edge case. if( l_p < h_p ) { // Non equivalent h_0 and h_p. - my_los->insert(l_p,10004); + // my_los->insert(l_p,10004); + my_los->insert(l_p,h_p); my_his->insert(h_p,10004); } } @@ -787,10 +810,12 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) my_his->freeRange(l_m,h_0); // TODO NOTE: When predecessor is used, you have to check to see if pred(b) is less than inf(interval). if( l_m < h_m ) { // If l_m and l_0 are "equivalent", so the current interval wins and we ignore interval_m. - my_los->insert(l_m,100052); + // my_los->insert(l_m,100052); + my_los->insert(l_m,h_m); my_his->insert(h_m,100052); } - my_los->insert(l_0,100052); + // my_los->insert(l_0,100052); + my_los->insert(l_0,h_0); my_his->insert(h_0,100052); // TODO Subtle bug? Need to verify edge case. // update for iteration lo1 = l_p; @@ -808,7 +833,8 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) if(not done) { // The new interval goes at the end of the skiplists. // Case 6. We're at the top. Just add. - my_los->insert(lo1,10006); + // my_los->insert(lo1,10006); + my_los->insert(lo1,hi1); my_his->insert(hi1,10006); done = true; } @@ -834,7 +860,7 @@ void HtmRangeMultiLevel::mergeRange(const Key lo, const Key hi) */ void HtmRangeMultiLevel::addRange(const Key lo, const Key hi) { -// my_los->insert(lo, (Value) 0); // TODO Consider doing something useful with (Value)... +// my_los->insert(lo, (Value) 0); // TODO Consider doing something useful with (Value)... Like storing hi... // my_his->insert(hi, (Value) 0); // cout << "x200: " << hex << lo << " " << hi << endl; // cout << "x201: " << (lo == hi) << endl; @@ -1103,7 +1129,8 @@ void HtmRangeMultiLevel::CompressionPass() { Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); // oldLo..newLoPredecessor; newLo..hi0 Modify the skiplists. my_his->insert(newLoPredecessor,1024); - my_los->insert(newLo,1024); + // my_los->insert(newLo,1024); + my_los->insert(newLo,newLoPredecessor); // Set lists to the new lo my_los->reset(); my_his->reset(); // TODO Until we know better, start over. Bad, bad, bad. // TODO Perhaps instead try a find or a search that would set the iterators. @@ -1119,11 +1146,13 @@ void HtmRangeMultiLevel::CompressionPass() { } my_los->free(oldLo); --oldLo; // Reduce level - my_los->insert(oldLo,1025); + // my_los->insert(oldLo,1025); + my_los->insert(oldLo,my_his->getkey()); // What's the current key (for my_his)? // TODO Don't worry. Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); if(newLoPredecessor != hi0) { my_his->insert(newLoPredecessor,1025); - my_los->insert(newLo,1025); + // my_los->insert(newLo,1025); + my_los->insert(newLo,newLoPredecessor); } my_los->reset(); my_his->reset(); // TODO Reset is too drastic. Prefer to step back a little... Bad, bad, bad. } From 7ce3648c6d325656c97fd76bb620e1ac0ec4c4af Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Mon, 16 Sep 2019 12:33:32 -0400 Subject: [PATCH 09/31] Sketching for python api for intersection. --- include/PySTARE.h | 3 ++- src/PySTARE.cpp | 15 ++++++++++++++- src/PySTARE.i | 9 +++++++++ 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/include/PySTARE.h b/include/PySTARE.h index 133a3b4..0e79800 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -29,12 +29,13 @@ static STARE stare; // Spatial -void from_latlon(double* lat, int len_lat, double * lon, int len_lon, int64_t* indices, int level); +void from_latlon(double* lat, int len_lat, double * lon, int len_lon, int64_t* indices, int level); void to_latlon(int64_t* indices, int len, double* lat, double* lon); void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* levels); void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); void to_area(int64_t* indices, int len, double* areas); +// void intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 990fdc6..e6bd9a4 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -49,12 +49,25 @@ void to_triangle(int64_t* indices, int len) { } void to_area(int64_t* indices, int len, double* areas) { - for (int i=0; itoSpatialIntervals(); + leni = result->size(); + for(int i=0; i Date: Mon, 16 Sep 2019 13:26:55 -0400 Subject: [PATCH 10/31] SpatialRange intersect test improved for clarity. --- tests/CUTE/SpatialRange_test.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/CUTE/SpatialRange_test.cpp b/tests/CUTE/SpatialRange_test.cpp index ae81a4b..355c080 100644 --- a/tests/CUTE/SpatialRange_test.cpp +++ b/tests/CUTE/SpatialRange_test.cpp @@ -213,6 +213,19 @@ void SpatialRange_test () { #undef DIAG } + // TODO Write many more tests & consider edge cases. + if(true) { + STARE_ArrayIndexSpatialValue siv1[2] = { 0x0000000000000008, 0x000067ffffffffff }; + STARE_ArrayIndexSpatialValue siv2[2] = { 0x000030000000000a, 0x0000907fffffffff }; + STARE_SpatialIntervals sis1(siv1,siv1+2); + STARE_SpatialIntervals sis2(siv2,siv2+2); + SpatialRange r1(sis1), r2(sis2); + SpatialRange *ri = r1 & r2; + STARE_SpatialIntervals result = ri->toSpatialIntervals(); + ASSERT_EQUAL(0x000030000000000a,result[0]); + ASSERT_EQUAL(0x000067ffffffffff,result[1]); + } + // TODO Add intersection tests. } From daa661e1ca6c98ea5a0137602adbdfd7fc00ccc4 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Tue, 17 Sep 2019 17:43:02 -0400 Subject: [PATCH 11/31] Adding python functions. Fixed SetPolyCorner bug that broke hull. Also added switch for multi-level intersection. Somewhat odd behavior with single-level-ish intersection result. --- contrib/python/example-1.py | 9 +++++ include/HtmRangeMultiLevel.h | 2 +- include/PySTARE.h | 3 ++ include/SpatialRange.h | 6 ++-- src/HtmRangeMultiLevel.cpp | 10 ++++-- src/PySTARE.cpp | 34 ++++++++++++++++-- src/PySTARE.i | 37 +++++++++++++++++++ src/STARE.C | 3 +- src/SpatialInterface.cpp | 22 +++++++++++- src/SpatialRange.cpp | 4 +-- tests/CUTE/SpatialInterface_test.cpp | 54 ++++++++++++++++++++++++++++ tests/CUTE/SpatialRange_test.cpp | 41 +++++++++++++++++---- 12 files changed, 205 insertions(+), 20 deletions(-) create mode 100644 contrib/python/example-1.py diff --git a/contrib/python/example-1.py b/contrib/python/example-1.py new file mode 100644 index 0000000..986ea6a --- /dev/null +++ b/contrib/python/example-1.py @@ -0,0 +1,9 @@ + STARE_ArrayIndexSpatialValue siv1[2] = { 0x0000000000000008, 0x000067ffffffffff }; + STARE_ArrayIndexSpatialValue siv2[2] = { 0x000030000000000a, 0x0000907fffffffff }; + + +import numpy as np +import pystare as ps +a = np.array([0x0000000000000008, 0x000030000000000a, 0x000067ffffffffff, 0x000070000000000a, 0x0000907fffffffff],dtype=np.int64) +ps.from_intervals(a) + diff --git a/include/HtmRangeMultiLevel.h b/include/HtmRangeMultiLevel.h index 18d1cd4..ab0f787 100644 --- a/include/HtmRangeMultiLevel.h +++ b/include/HtmRangeMultiLevel.h @@ -124,7 +124,7 @@ class LINKAGE HtmRangeMultiLevel { }; HtmRangeMultiLevel *HtmRangeMultiLevelAtLevelFromIntersection(HtmRangeMultiLevel *range2, int htmIdLevel=-1); - HtmRangeMultiLevel *RangeFromIntersection(HtmRangeMultiLevel *range2, int force_htmIdLevel=-1); + HtmRangeMultiLevel *RangeFromIntersection(HtmRangeMultiLevel *range2, bool compress = true, int force_htmIdLevel=-1); int contains(Key a, Key b); // TODO KeyPair contains(Key a, Key b); diff --git a/include/PySTARE.h b/include/PySTARE.h index 0e79800..cbb13e6 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -35,6 +35,9 @@ void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* le void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); void to_area(int64_t* indices, int len, double* areas); + +void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); + // void intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); // Temporal diff --git a/include/SpatialRange.h b/include/SpatialRange.h index d367113..485ee0d 100644 --- a/include/SpatialRange.h +++ b/include/SpatialRange.h @@ -43,13 +43,13 @@ class SpatialRange { }; void reset() { range->reset(); } // range not null? void purge() { range->purge(); } // what if range null? - }; -SpatialRange* sr_intersect(const SpatialRange& a, const SpatialRange& b); +SpatialRange* sr_intersect(const SpatialRange& a, const SpatialRange& b, bool compress = false); inline SpatialRange* operator& ( const SpatialRange& a, const SpatialRange& b) { - return new SpatialRange(new HstmRange(a.range->range->RangeFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. + return sr_intersect(a,b); + // return new SpatialRange(new HstmRange(a.range->range->RangeFromIntersection(b.range->range))); // NOTE mlr Probably about the safest way to inst. SpatialRange. } void SpatialRange_test(); diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index d975068..49a6e4a 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -149,7 +149,8 @@ KeyPair HRML_AtLevelFromMultiLevel(uint64 htmIdLevel, Key lo, Key hi, uint64 lev * TODO Replace double-nested-loop with SkipLists search or find functionality. * */ -HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel *range2, int force_htmIdLevel) { +HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection( + HtmRangeMultiLevel *range2, bool compress, int force_htmIdLevel ) { HtmRangeMultiLevel *range1 = this; // Just an alias if((!range1)||(!range2)) return 0; if((range1->nranges()<=0)||(range2->nranges()<=0)) return 0; @@ -212,7 +213,12 @@ HtmRangeMultiLevel *HtmRangeMultiLevel::RangeFromIntersection(HtmRangeMultiLevel // cout << "d" << flush; // cout << "d nr " << resultRange->nranges() << endl << flush; // cout << "d rr " << hex << resultRange << dec << endl << flush; - if(resultRange->nranges()>0)resultRange->defrag(); + if(resultRange->nranges()>0) { + if(compress) { + resultRange->CompressionPass(); + } + resultRange->defrag(); + } // cout << "e" << flush; return resultRange; } diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index e6bd9a4..de55fd1 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -37,8 +37,7 @@ void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* le void to_level(int64_t* indices, int len, int* levels) { for (int i=0; i 0 ) { + // b is in the middle + polyCorners_[1].c_ = v; + } else { + // a is in the middle + polyCorners_[0].c_ = v; + } + } // dot_ab_av*dot_ab_vb == 0 is an error and > 0 means v is in the middle. } } else { // @@ -556,6 +575,7 @@ htmInterface::setPolyCorner(SpatialVector &v) { // if( (polyCorners_[i].c_ ^ polyCorners_[i+1==len ? 0 : i+1].c_)*v > tol2 ) { float64 delta = (polyCorners_[i].c_ ^ polyCorners_[i+1==len ? 0 : i+1].c_)*v; +// #define DIAG #ifdef DIAG DIAG_OUT << i << " i," diff --git a/src/SpatialRange.cpp b/src/SpatialRange.cpp index 6d00725..7c0b316 100644 --- a/src/SpatialRange.cpp +++ b/src/SpatialRange.cpp @@ -109,8 +109,8 @@ STARE_SpatialIntervals SpatialRange::toSpatialIntervals() { /* * Odd. The following does not seem to work if we just return the SpatialRange itself. Some of the pointers seem to be either corrupted or eliminated. */ -SpatialRange *sr_intersect(const SpatialRange&a, const SpatialRange& b) { - HstmRange *range = new HstmRange(a.range->range->RangeFromIntersection(b.range->range)); // NOTE mlr Probably about the safest way to inst. SpatialRange. +SpatialRange *sr_intersect(const SpatialRange&a, const SpatialRange& b, bool compress) { + HstmRange *range = new HstmRange(a.range->range->RangeFromIntersection(b.range->range,compress)); // NOTE mlr Probably about the safest way to inst. SpatialRange. // #define DIAG #ifdef DIAG KeyPair kp; range->reset(); range->getNext(kp); diff --git a/tests/CUTE/SpatialInterface_test.cpp b/tests/CUTE/SpatialInterface_test.cpp index 01386e5..de0e65e 100644 --- a/tests/CUTE/SpatialInterface_test.cpp +++ b/tests/CUTE/SpatialInterface_test.cpp @@ -114,5 +114,59 @@ void SpatialInterface_test() { v.setLatLonDegrees( 45, -140);ASSERT_EQUALDM("4polyC.3",htm->polyCorners_[3].c_,v,1.0e-12); delete htm; } + + if( true ) { + + cout << endl << "----" << endl << endl << endl << flush; + + // Try a grid... + int nLat = 11,nLon = 11, nLatLon=nLat*nLon; + float64 lat[nLat], lon[nLon]; + SpatialVector vs[nLatLon], cs[nLatLon], vs1[nLatLon]; + for( int i = 0; i < nLon; ++i ) { + lon[i] = i*5.0; + } + for( int j = 0; j < nLat; ++j ) { + lat[j] = j*5.0; + } + int k = 0; + for( int i = 0; i < nLon; ++i ) { + for( int j = 0; j < nLat; ++j ) { + vs[k++].setLatLonDegrees(lat[j], lon[i]); + } + } + + STARE index; + htmInterface *htm; + int resolution_level = 4; + htm = new htmInterface( + index.getIndex(resolution_level).getMaxlevel(), + index.getIndex(resolution_level).getBuildLevel(), + index.getIndex(resolution_level).getRotation()); + htm->polyCorners_.clear(); + SpatialVector v; + + k = 0; + for( int i=0; isetPolyCorner(vs[k]); + cout << k << " k, hpc size " << htm->polyCorners_.size() << endl << flush; + float64 lat_,lon_; vs[k].getLatLonDegrees(lat_,lon_); + cout << k << " k, try " + << lat_ << " " << lon_ << endl << flush; + for(int l=0; lpolyCorners_.size(); ++l) { + htm->polyCorners_[l].c_.getLatLonDegrees(lat_,lon_); + cout << l << " l, hpc " + << lat_ << " " << lon_ << endl << flush; + } + cout << "-----" << endl << flush; + ++k; + // if(k>5) exit(1); + } + } + + delete htm; + + } } } diff --git a/tests/CUTE/SpatialRange_test.cpp b/tests/CUTE/SpatialRange_test.cpp index 355c080..1ab94fb 100644 --- a/tests/CUTE/SpatialRange_test.cpp +++ b/tests/CUTE/SpatialRange_test.cpp @@ -9,6 +9,8 @@ #include "Test.h" +// #define DIAG + void SpatialRange_test () { STARE index; @@ -72,13 +74,17 @@ void SpatialRange_test () { ASSERT_EQUAL(sis[0],sis_out[0]); ASSERT_EQUAL(sis[1],sis_out[1]); - if(true) { // TODO Fix the following to use the new 0/1 variables defined above. + +#undef DIAG +#ifndef DIAG #define SIVOUT(m,siv) -// #define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; #define SISOUT(m,sis) -// #define SISOUT(m,sis) { cout << endl << m << endl; for(int i=0; i < sis.size(); ++i ) { cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; }} +#else +#define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SISOUT(m,sis) { cout << endl << m << endl; for(int i=0; i < sis.size(); ++i ) { cout << i << " i,si: " << setw(16) << setfill('0') << hex << sis[i] << dec << endl << flush; }} +#endif // #define DIAGOUT // #define DIAGOUT(m) {cout << m << endl << flush;} // Note level is 8. @@ -205,12 +211,9 @@ void SpatialRange_test () { SISOUT("deltaSis",deltaSis); // cout << 400 << endl << flush; + ASSERT_EQUAL(2,deltaSis.size()); ASSERT_EQUAL(sis1[0],deltaSis[0]); ASSERT_EQUAL(sis0[1],deltaSis[1]); - -#undef SIVOUT -#undef SISOUT -#undef DIAG } // TODO Write many more tests & consider edge cases. @@ -221,11 +224,35 @@ void SpatialRange_test () { STARE_SpatialIntervals sis2(siv2,siv2+2); SpatialRange r1(sis1), r2(sis2); SpatialRange *ri = r1 & r2; + // SpatialRange *ri = sr_intersect(r1, r2); + STARE_SpatialIntervals result = ri->toSpatialIntervals(); + ASSERT_EQUAL(2,result.size()); ASSERT_EQUAL(0x000030000000000a,result[0]); ASSERT_EQUAL(0x000067ffffffffff,result[1]); + + ri->range->range->CompressionPass(); + result = ri->toSpatialIntervals(); + ASSERT_EQUAL(4,result.size()); + ASSERT_EQUAL(0x0000300000000008,result[0]); + ASSERT_EQUAL(0x00003fffffffffff,result[1]); + ASSERT_EQUAL(0x0000400000000007,result[2]); + ASSERT_EQUAL(0x0000600000000008,result[3]); + + delete ri; + ri = sr_intersect(r1, r2, true); // Run a compression pass on the range. + result.clear(); + result = ri->toSpatialIntervals(); + ASSERT_EQUAL(4,result.size()); + ASSERT_EQUAL(0x0000300000000008,result[0]); + ASSERT_EQUAL(0x00003fffffffffff,result[1]); + ASSERT_EQUAL(0x0000400000000007,result[2]); + ASSERT_EQUAL(0x0000600000000008,result[3]); } +#undef SIVOUT +#undef SISOUT +#undef DIAG // TODO Add intersection tests. } From 0848af8c8d3adae8061aa8b60d25818a89c3e3c5 Mon Sep 17 00:00:00 2001 From: griessbaum Date: Wed, 18 Sep 2019 15:50:51 +0200 Subject: [PATCH 12/31] added intersect dummy to pystare --- include/PySTARE.h | 3 ++- src/PySTARE.cpp | 25 +++++++++++++----------- src/PySTARE.i | 50 ++++++++++++++++++++++++++++------------------- 3 files changed, 46 insertions(+), 32 deletions(-) diff --git a/include/PySTARE.h b/include/PySTARE.h index cbb13e6..74eb97d 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -25,6 +25,7 @@ #include #include "STARE.h" +#include "SpatialRange.h" static STARE stare; @@ -38,7 +39,7 @@ void to_area(int64_t* indices, int len, double* areas); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); -// void intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index de55fd1..78c2f69 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -56,6 +56,8 @@ void to_area(int64_t* indices, int len, double* areas) { /** * Go from an array of [id|id..term] to a pair of arrays [id][term] to aid tests for inclusion. */ + + void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ) { // cout << "len: " << len << endl << flush; // if(false) { @@ -83,18 +85,19 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ // } } -/* -void intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { - SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); - SpatialRange r1(si1), r2(si2); - SpatialRange *ri = r1 & r2; - SpatialIntervals result = ri->toSpatialIntervals(); - leni = result->size(); - for(int i=0; itoSpatialIntervals(); + //leni = result.size(); + //for(int i=0; i Date: Wed, 18 Sep 2019 10:12:29 -0400 Subject: [PATCH 13/31] Updated notes, examined pystare from_intervals. --- NOTES | 18 ++++++++++++++++++ contrib/python/example-1.py | 7 +++---- include/SpatialInterface.h | 3 ++- src/PySTARE.cpp | 4 ++++ src/STARE.C | 4 ++-- tests/CUTE/SpatialInterface_test.cpp | 14 +++++++------- 6 files changed, 36 insertions(+), 14 deletions(-) diff --git a/NOTES b/NOTES index 9fe37b8..71ba105 100644 --- a/NOTES +++ b/NOTES @@ -1,6 +1,24 @@ * 2019-09-X Version 0.8.0 +Added a number of tests using Constraint. + + + +Improved intersection behavior by injecting dependency on varlen, i.e. +interior vs. boundary calculation. Setting varlen to false tends to +force calculation all the way down to leaf, whereas varlen true +stops on reaching a fully enclosed trixel, leading to the multi-resolution +partitioning. + +Corrected bug in SetPolyCorner of ConvexHull calculation. 2-corner list would +throw out co-linear third point without checking if it extended the side +or not. + +Changed HtmRangeMultiLevel::RangeFromIntersection to eliminate iteration over +the range in favor of using the SkipList's findMAX. Also eliminated superfluous +iteration after intersection is found. + Added SpatialRange to use SkipList to generate intersections. Note the intersect doesn't use find or search, so it's not clear we're getting the benefit of the SkipLists. diff --git a/contrib/python/example-1.py b/contrib/python/example-1.py index 986ea6a..463d8c1 100644 --- a/contrib/python/example-1.py +++ b/contrib/python/example-1.py @@ -1,9 +1,8 @@ - STARE_ArrayIndexSpatialValue siv1[2] = { 0x0000000000000008, 0x000067ffffffffff }; - STARE_ArrayIndexSpatialValue siv2[2] = { 0x000030000000000a, 0x0000907fffffffff }; - import numpy as np import pystare as ps a = np.array([0x0000000000000008, 0x000030000000000a, 0x000067ffffffffff, 0x000070000000000a, 0x0000907fffffffff],dtype=np.int64) -ps.from_intervals(a) +starts,ends = ps.from_intervals(a) +for i in range(starts.size): + print(i,hex(starts[i]),hex(ends[i])) diff --git a/include/SpatialInterface.h b/include/SpatialInterface.h index 7465dab..047a16b 100644 --- a/include/SpatialInterface.h +++ b/include/SpatialInterface.h @@ -48,10 +48,11 @@ struct htmPolyCorner { /** htmInterface class. + The SpatialInterface class contains all methods to interface the HTM index with external applications. - // TODO Consider renaming htmInterface to SpatialInterface + // TODO Consider renaming htmInterface to SpatialInterface, but maybe "htm" is a good indicator that we're working at a low level with this interface. */ diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index de55fd1..2b49ca5 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -76,6 +76,10 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ } ++i; // Next } while(iconvexHull(points,hullSteps,true); // TODO FIX interiorp = false is broken + HTMRangeValueVector htmRangeVector = htm->convexHull(points,hullSteps,true); // Compress result // cout << dec << "a3000 hrv.size: " << htmRangeVector.size() << endl << flush; @@ -471,7 +471,7 @@ STARE_SpatialIntervals STARE::ConvexHull(LatLonDegrees64ValueVector points,int f uint64 lo = ValueFromHtmID(htmRangeVector[i].lo); // TODO Should this be a function? cover.push_back(lo); uint64 hi; - if( htmRangeVector[i].lo != htmRangeVector[i].lo ) { + if( htmRangeVector[i].lo != htmRangeVector[i].hi ) { hi = sTerminator(ValueFromHtmID(htmRangeVector[i].hi)); cover.push_back(hi); } diff --git a/tests/CUTE/SpatialInterface_test.cpp b/tests/CUTE/SpatialInterface_test.cpp index de0e65e..ed1bce0 100644 --- a/tests/CUTE/SpatialInterface_test.cpp +++ b/tests/CUTE/SpatialInterface_test.cpp @@ -117,7 +117,7 @@ void SpatialInterface_test() { if( true ) { - cout << endl << "----" << endl << endl << endl << flush; + // cout << endl << "----" << endl << endl << endl << flush; // Try a grid... int nLat = 11,nLon = 11, nLatLon=nLat*nLon; @@ -150,16 +150,16 @@ void SpatialInterface_test() { for( int i=0; isetPolyCorner(vs[k]); - cout << k << " k, hpc size " << htm->polyCorners_.size() << endl << flush; + // cout << k << " k, hpc size " << htm->polyCorners_.size() << endl << flush; float64 lat_,lon_; vs[k].getLatLonDegrees(lat_,lon_); - cout << k << " k, try " - << lat_ << " " << lon_ << endl << flush; + // cout << k << " k, try " + // << lat_ << " " << lon_ << endl << flush; for(int l=0; lpolyCorners_.size(); ++l) { htm->polyCorners_[l].c_.getLatLonDegrees(lat_,lon_); - cout << l << " l, hpc " - << lat_ << " " << lon_ << endl << flush; + // cout << l << " l, hpc " + // << lat_ << " " << lon_ << endl << flush; } - cout << "-----" << endl << flush; + // cout << "-----" << endl << flush; ++k; // if(k>5) exit(1); } From 2f6f45cb87ad1e61c76c4e2678a75e0a08cdb273 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Wed, 18 Sep 2019 11:25:44 -0400 Subject: [PATCH 14/31] Corrected linking order in pystare build and added sr_intersection. --- contrib/python/test.py | 42 ++++++++++++++++++++++++++++++++ src/CMakeLists.txt | 2 +- src/PySTARE.cpp | 19 +++++++-------- src/PySTARE.i | 7 ++++-- tests/CUTE/SpatialRange_test.cpp | 23 +++++++++++++++++ 5 files changed, 80 insertions(+), 13 deletions(-) create mode 100644 contrib/python/test.py diff --git a/contrib/python/test.py b/contrib/python/test.py new file mode 100644 index 0000000..0b9f81c --- /dev/null +++ b/contrib/python/test.py @@ -0,0 +1,42 @@ +#!/usr/bin/python3 + +import numpy +import pystare + + +lat = numpy.array([30,45,60], dtype=numpy.double) +lon = numpy.array([45,60,10], dtype=numpy.double) + + +indices = pystare.from_latlon(lat, lon, 12) +print('0 indices: ',[hex(i) for i in indices]) + +lat, lon = pystare.to_latlon(indices) +print(lat, lon) + +lat, lon, level = pystare.to_latlonlevel(indices) +print(lat, lon, level) + +level = pystare.to_level(indices) +print(level) + +area = pystare.to_area(indices) +print(area) + +datetime = numpy.array(['2002-02-03T13:56:03.172', '2016-01-05T17:26:00.172'], dtype=numpy.datetime64) +print(datetime) + +index = pystare.from_utc(datetime, 6) +print(index) + +intersected = numpy.zeros([3], dtype=numpy.int64) +#? leni = 0 +pystare._intersect(indices, indices, intersected) +print('1 indices: ',[hex(i) for i in indices]) +print('1 intersected: ',[hex(i) for i in intersected] ) + +intersected = pystare.intersect(indices, indices) +print('2 intersected: ',[hex(i) for i in intersected]) + + + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d16f5c8..d5fa51c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -80,9 +80,9 @@ if (SWIG_FOUND AND POLICY CMP0078 AND POLICY CMP0086) # swig_add_library( PySTARE TYPE SHARED LANGUAGE python SOURCES PySTARE.i ) swig_add_library( pystare TYPE SHARED LANGUAGE python SOURCES PySTARE.i ) - swig_link_libraries( pystare STARE ) swig_link_libraries( pystare PySTAREA ) swig_link_libraries( pystare ${PYTHON_LIBRARIES} ) + swig_link_libraries( pystare STARE ) endif() install(TARGETS STARE DESTINATION lib) diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 975f1e6..94d00ac 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -91,18 +91,17 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); - intersection[0] = 69; - //SpatialRange r1(si1), r2(si2); - //SpatialRange *ri = r1 & r2; - //STARE_SpatialIntervals result = ri->toSpatialIntervals(); - //leni = result.size(); - //for(int i=0; itoSpatialIntervals(); + leni = result.size(); + for(int i=0; itoSpatialIntervals(); + leni = result.size(); + for(int i=0; i Date: Wed, 18 Sep 2019 11:37:40 -0400 Subject: [PATCH 15/31] Added multilevel capability to intersection. --- include/PySTARE.h | 1 + src/PySTARE.cpp | 13 +++++++++++++ src/PySTARE.i | 7 +++++-- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/include/PySTARE.h b/include/PySTARE.h index 74eb97d..b2f8fa8 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -40,6 +40,7 @@ void to_area(int64_t* indices, int len, double* areas); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _intersect_multilevel(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 94d00ac..d3a760d 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -102,6 +102,19 @@ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_ } } +void _intersect_multilevel(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { + STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); + // intersection[0] = 69; + SpatialRange r1(si1), r2(si2); + // SpatialRange *ri = r1 & r2; + SpatialRange *ri = sr_intersect(r1,r2,true); + STARE_SpatialIntervals result = ri->toSpatialIntervals(); + leni = result.size(); + for(int i=0; i Date: Wed, 18 Sep 2019 11:38:51 -0400 Subject: [PATCH 16/31] Corrected multilevel intersect in pystare (boolean error). --- contrib/python/test.py | 2 ++ src/PySTARE.i | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index 0b9f81c..9aa0bda 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -38,5 +38,7 @@ intersected = pystare.intersect(indices, indices) print('2 intersected: ',[hex(i) for i in intersected]) +intersected = pystare.intersect(indices, indices, multilevel=True) +print('3 intersected: ',[hex(i) for i in intersected]) diff --git a/src/PySTARE.i b/src/PySTARE.i index e667563..4f875ff 100644 --- a/src/PySTARE.i +++ b/src/PySTARE.i @@ -280,9 +280,9 @@ def intersect(indices1, indices2, multilevel=True): intersected = numpy.full([out_length],-1,dtype=numpy.int64) leni = 0 if(multilevel): - _intersect(indices1, indices2, intersected) - else: _intersect_multilevel(indices1, indices2, intersected) + else: + _intersect(indices1, indices2, intersected) endarg = numpy.argmax(intersected < 0) # intersected = intersected.trim_zeros() intersected = intersected[:endarg] From 9e8329cb057477a4e203420ee656fa813ce348c5 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Wed, 18 Sep 2019 11:44:05 -0400 Subject: [PATCH 17/31] Changed the term level to resolution. --- contrib/python/test.py | 2 +- include/PySTARE.h | 2 +- src/PySTARE.cpp | 2 +- src/PySTARE.i | 8 +++----- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index 9aa0bda..34f3870 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -38,7 +38,7 @@ intersected = pystare.intersect(indices, indices) print('2 intersected: ',[hex(i) for i in intersected]) -intersected = pystare.intersect(indices, indices, multilevel=True) +intersected = pystare.intersect(indices, indices, multiresolution=True) print('3 intersected: ',[hex(i) for i in intersected]) diff --git a/include/PySTARE.h b/include/PySTARE.h index b2f8fa8..4e6e475 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -40,7 +40,7 @@ void to_area(int64_t* indices, int len, double* areas); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); -void _intersect_multilevel(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index d3a760d..7c997e2 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -102,7 +102,7 @@ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_ } } -void _intersect_multilevel(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { +void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); // intersection[0] = 69; SpatialRange r1(si1), r2(si2); diff --git a/src/PySTARE.i b/src/PySTARE.i index 4f875ff..940c5fc 100644 --- a/src/PySTARE.i +++ b/src/PySTARE.i @@ -265,22 +265,20 @@ (int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ) } - %pythonprepend from_utc(int64_t*, int, int64_t*, int) %{ import numpy datetime = datetime.astype(numpy.int64) %} - %pythoncode %{ import numpy -def intersect(indices1, indices2, multilevel=True): +def intersect(indices1, indices2, multiresolution=True): out_length = 2*max(len(indices1),len(indices2)) # intersected = numpy.zeros([out_length], dtype=numpy.int64) intersected = numpy.full([out_length],-1,dtype=numpy.int64) leni = 0 - if(multilevel): - _intersect_multilevel(indices1, indices2, intersected) + if(multiresolution): + _intersect_multiresolution(indices1, indices2, intersected) else: _intersect(indices1, indices2, intersected) endarg = numpy.argmax(intersected < 0) From fca51d9d76bb4fbe62e3e78f3e16273820c84a7a Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Wed, 18 Sep 2019 11:51:13 -0400 Subject: [PATCH 18/31] Added some intersection tests. --- contrib/python/test.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/contrib/python/test.py b/contrib/python/test.py index 34f3870..6c191ed 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -41,4 +41,10 @@ intersected = pystare.intersect(indices, indices, multiresolution=True) print('3 intersected: ',[hex(i) for i in intersected]) +indices1 = numpy.array([indices[1]], dtype=numpy.int64) +intersected = pystare.intersect(indices, indices1, multiresolution=True) +print('4 intersected: ',[hex(i) for i in intersected]) +indices1 = numpy.array([0x100000000000000c], dtype=numpy.int64) +intersected = pystare.intersect(indices, indices1, multiresolution=True) +print('5 intersected: ',[hex(i) for i in intersected]) From c67f92fdf036e4812b08ee36892fc668387b738a Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Wed, 18 Sep 2019 15:24:56 -0400 Subject: [PATCH 19/31] Added to_vertices to pystare. --- contrib/python/test.py | 26 +++++++++++++++++++++++++ include/PySTARE.h | 1 + src/PySTARE.cpp | 20 ++++++++++++++++++++ src/PySTARE.i | 40 +++++++++++++++++++++++++++++++++++++-- src/STARE.C | 24 ++++++++++++++++++++++- tests/CUTE/STARE_test.cpp | 16 ++++++++++++++++ 6 files changed, 124 insertions(+), 3 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index 6c191ed..220ad5a 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -48,3 +48,29 @@ indices1 = numpy.array([0x100000000000000c], dtype=numpy.int64) intersected = pystare.intersect(indices, indices1, multiresolution=True) print('5 intersected: ',[hex(i) for i in intersected]) +print("") + +vertices0 = numpy.zeros([3], dtype=numpy.int64) +vertices1 = numpy.zeros([3], dtype=numpy.int64) +vertices2 = numpy.zeros([3], dtype=numpy.int64) + +vertices0,vertices1,vertices2 = pystare.to_vertices(indices) +print('vertices0: ',[hex(i) for i in vertices0]) +print('vertices1: ',[hex(i) for i in vertices1]) +print('vertices2: ',[hex(i) for i in vertices2]) + +vertices0_lats,vertices0_lons = pystare.to_latlon(vertices0) +vertices1_lats,vertices1_lons = pystare.to_latlon(vertices1) +vertices2_lats,vertices2_lons = pystare.to_latlon(vertices2) + +for i in range(len(vertices0_lats)): + print(i," vert0 lat,lon: ",vertices0_lats[i],vertices0_lons[i]) + print(i," vert1 lat,lon: ",vertices1_lats[i],vertices1_lons[i]) + print(i," vert2 lat,lon: ",vertices2_lats[i],vertices2_lons[i]) + print("") + + + + + + diff --git a/include/PySTARE.h b/include/PySTARE.h index 4e6e475..6254645 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -35,6 +35,7 @@ void to_latlon(int64_t* indices, int len, double* lat, double* lon); void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* levels); void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); +void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2); void to_area(int64_t* indices, int len, double* areas); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 7c997e2..900f0a5 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -47,6 +47,26 @@ void to_triangle(int64_t* indices, int len) { } } +void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2) { + for (int i=0; i Date: Wed, 18 Sep 2019 17:19:07 -0400 Subject: [PATCH 20/31] Added to_compressed_range to make it easier to construct ranges later. --- contrib/python/test.py | 10 ++++++++++ include/PySTARE.h | 1 + src/PySTARE.cpp | 9 +++++++++ src/PySTARE.i | 26 +++++++++++++++++++------- src/STARE.C | 2 ++ 5 files changed, 41 insertions(+), 7 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index 220ad5a..8d9647a 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -53,11 +53,13 @@ vertices0 = numpy.zeros([3], dtype=numpy.int64) vertices1 = numpy.zeros([3], dtype=numpy.int64) vertices2 = numpy.zeros([3], dtype=numpy.int64) +print('') vertices0,vertices1,vertices2 = pystare.to_vertices(indices) print('vertices0: ',[hex(i) for i in vertices0]) print('vertices1: ',[hex(i) for i in vertices1]) print('vertices2: ',[hex(i) for i in vertices2]) +print('') vertices0_lats,vertices0_lons = pystare.to_latlon(vertices0) vertices1_lats,vertices1_lons = pystare.to_latlon(vertices1) @@ -69,6 +71,14 @@ print(i," vert2 lat,lon: ",vertices2_lats[i],vertices2_lons[i]) print("") +indices1 = numpy.array([0,0,0], dtype=numpy.int64) +pystare._to_compressed_range(indices,indices1) +print('_indices1: ',[hex(i) for i in indices1]) + +indices1 = numpy.array([0,0,0], dtype=numpy.int64) +indices1 = pystare.to_compressed_range(indices) +print('indices1: ',list(map(hex,indices1))) + diff --git a/include/PySTARE.h b/include/PySTARE.h index 6254645..b1e7093 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -38,6 +38,7 @@ void to_triangle(int64_t* indices, int len); void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2); void to_area(int64_t* indices, int len, double* areas); +void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 900f0a5..8a3fbe9 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -109,6 +109,15 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ // } } +void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri) { + STARE_SpatialIntervals si(indices, indices+len); + SpatialRange r(si); + STARE_SpatialIntervals result = r.toSpatialIntervals(); + for(int i=0; i 0): + endarg = endarg + 1 + # endarg = numpy.argmax(range_indices < 0) + range_indices = range_indices[:endarg] + return range_indices + def intersect(indices1, indices2, multiresolution=True): out_length = 2*max(len(indices1),len(indices2)) - # intersected = numpy.zeros([out_length], dtype=numpy.int64) intersected = numpy.full([out_length],-1,dtype=numpy.int64) leni = 0 if(multiresolution): @@ -318,10 +332,8 @@ def intersect(indices1, indices2, multiresolution=True): else: _intersect(indices1, indices2, intersected) endarg = numpy.argmax(intersected < 0) - # intersected = intersected.trim_zeros() intersected = intersected[:endarg] return intersected - %} %include "PySTARE.h" diff --git a/src/STARE.C b/src/STARE.C index 67e00c3..205f338 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -284,6 +284,8 @@ float64 STARE::AreaFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int re return sIndexes[resolutionLevel].areaByHtmId(htmID); } +/// TODO STARE::InfoFromValue // Get all the info from a volume, so we only translate between forms once. + STARE_ArrayIndexSpatialValue sTerminator(STARE_ArrayIndexSpatialValue spatialStareId) { EmbeddedLevelNameEncoding leftJustifiedWithResolution; leftJustifiedWithResolution.setIdFromSciDBLeftJustifiedFormat(spatialStareId); From e716a1df510dc5c1bb46ad4a4e0a65ade60ce196 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Wed, 18 Sep 2019 21:56:17 -0400 Subject: [PATCH 21/31] Many additions to pystare. Including hull, vertices, and centroids. --- NOTES | 9 +++++-- contrib/python/test.py | 24 ++++++++++++++---- include/PySTARE.h | 3 ++- src/PySTARE.cpp | 17 ++++++++++++- src/PySTARE.i | 57 ++++++++++++++++++++++++++++++++++-------- src/STARE.C | 8 ++++++ 6 files changed, 99 insertions(+), 19 deletions(-) diff --git a/NOTES b/NOTES index 71ba105..5ad3f8c 100644 --- a/NOTES +++ b/NOTES @@ -1,9 +1,14 @@ * 2019-09-X Version 0.8.0 -Added a number of tests using Constraint. - +Many improvements to pystare. +- Vertices & Centroids +- ConvexHull +- SpatialRange lists +Some use is made of the STARE spatial index as a way to convey +location information. +Added a number of tests using Constraint. Improved intersection behavior by injecting dependency on varlen, i.e. interior vs. boundary calculation. Setting varlen to false tends to diff --git a/contrib/python/test.py b/contrib/python/test.py index 8d9647a..ba48367 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -3,11 +3,9 @@ import numpy import pystare - lat = numpy.array([30,45,60], dtype=numpy.double) lon = numpy.array([45,60,10], dtype=numpy.double) - indices = pystare.from_latlon(lat, lon, 12) print('0 indices: ',[hex(i) for i in indices]) @@ -48,27 +46,30 @@ indices1 = numpy.array([0x100000000000000c], dtype=numpy.int64) intersected = pystare.intersect(indices, indices1, multiresolution=True) print('5 intersected: ',[hex(i) for i in intersected]) -print("") +print('') vertices0 = numpy.zeros([3], dtype=numpy.int64) vertices1 = numpy.zeros([3], dtype=numpy.int64) vertices2 = numpy.zeros([3], dtype=numpy.int64) -print('') +centroid = numpy.zeros([3], dtype=numpy.int64) -vertices0,vertices1,vertices2 = pystare.to_vertices(indices) +vertices0,vertices1,vertices2,centroid = pystare.to_vertices(indices) print('vertices0: ',[hex(i) for i in vertices0]) print('vertices1: ',[hex(i) for i in vertices1]) print('vertices2: ',[hex(i) for i in vertices2]) +print('centroid: ',[hex(i) for i in centroid]) print('') vertices0_lats,vertices0_lons = pystare.to_latlon(vertices0) vertices1_lats,vertices1_lons = pystare.to_latlon(vertices1) vertices2_lats,vertices2_lons = pystare.to_latlon(vertices2) +centroid_lats,centroid_lons = pystare.to_latlon(centroid) for i in range(len(vertices0_lats)): print(i," vert0 lat,lon: ",vertices0_lats[i],vertices0_lons[i]) print(i," vert1 lat,lon: ",vertices1_lats[i],vertices1_lons[i]) print(i," vert2 lat,lon: ",vertices2_lats[i],vertices2_lons[i]) + print(i," centr lat,lon: ",centroid_lats[i],centroid_lons[i]) print("") indices1 = numpy.array([0,0,0], dtype=numpy.int64) @@ -79,6 +80,19 @@ indices1 = pystare.to_compressed_range(indices) print('indices1: ',list(map(hex,indices1))) +indices1 = numpy.zeros([1000], dtype=numpy.int64) +result_size = numpy.zeros([1], dtype=numpy.int) +pystare._to_hull_range(indices,8,indices1,result_size) +# print('hull indices1: ',list(map(hex,indices1[0:100]))) +print('hull result_size: ',result_size) +# print('hull indices1: ',list(map(hex,indices1[899:910]))) +indices1=indices1[0:result_size[0]] +print('hull trimmed size ',indices1.size) +print('') + +hull_indices = pystare.to_hull_range(indices,8,2000) +print('hull_indices size: ',hull_indices.size) +print('') diff --git a/include/PySTARE.h b/include/PySTARE.h index b1e7093..818aabe 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -35,10 +35,11 @@ void to_latlon(int64_t* indices, int len, double* lat, double* lon); void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* levels); void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); -void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2); +void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid); void to_area(int64_t* indices, int len, double* areas); void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri); +void _to_hull_range(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 8a3fbe9..203e5c4 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -47,7 +47,7 @@ void to_triangle(int64_t* indices, int len) { } } -void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2) { +void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid) { for (int i=0; i Date: Thu, 19 Sep 2019 09:20:49 -0400 Subject: [PATCH 22/31] Added pystare.cmp_spatial. Need to watch out for intervals vs. values. --- contrib/python/test.py | 38 ++++++++++++++++++++++++++++++---- include/PySTARE.h | 1 + src/PySTARE.cpp | 11 ++++++++++ src/PySTARE.i | 47 ++++++++++++++++++++++++++++++++++++++++-- 4 files changed, 91 insertions(+), 6 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index ba48367..ca9cf97 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -94,7 +94,37 @@ print('hull_indices size: ',hull_indices.size) print('') - - - - +cmp = numpy.zeros([9],dtype=numpy.int64) +pystare._cmp_spatial(indices,indices,cmp) +print('cmp input: ',indices) +print('cmp: ',cmp,' i.e. diagonal') +print('') +cmp = numpy.zeros([3*1],dtype=numpy.int64) +indices1 = numpy.zeros([1],dtype=numpy.int64) +indices1[0] = indices[1] +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +pystare._cmp_spatial(indices,indices1,cmp) +print('cmp: ',cmp) +print('') +cmp = numpy.zeros([3*2],dtype=numpy.int64) +indices1 = numpy.zeros([2],dtype=numpy.int64) +indices1[0] = indices[1] +indices1[1] = indices[1]+3 +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +pystare._cmp_spatial(indices,indices1,cmp) +print('cmp: ',cmp) +pystare._cmp_spatial(indices1,indices,cmp) +print('cmp: ',cmp) +print('') +mp = numpy.zeros([3*2],dtype=numpy.int64) +indices1 = numpy.zeros([2],dtype=numpy.int64) +indices1[0] = indices[1]-2 +indices1[1] = indices[1] +print('cmp input1: ',indices) +print('cmp input2: ',indices1) +# pystare._cmp_spatial(indices,indices1,cmp) +cmp1 = pystare.cmp_spatial(indices,indices1) +print('cmp1: ',cmp1) +print('') diff --git a/include/PySTARE.h b/include/PySTARE.h index 818aabe..7129f78 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -44,6 +44,7 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); +void _cmp_spatial(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* cmp, int len12); // Temporal void from_utc(int64_t *datetime, int len, int64_t *indices, int resolution); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 203e5c4..ae157db 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -159,6 +159,17 @@ void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, } } +void _cmp_spatial(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* cmp, int len12) { + STARE_ArrayIndexSpatialValues sivs1(indices1, indices1+len1), sivs2(indices2, indices2+len2); + int k=0; + for(int i=0; i Date: Fri, 20 Sep 2019 00:00:20 -0400 Subject: [PATCH 23/31] Intermediate push. Lot's of debugging. idByPoint problem. --- contrib/python/.gitignore | 1 + contrib/python/example-viz-1.py | 151 +++++++++++++++++++++++++++++++ contrib/python/test.py | 12 +++ include/SpatialGeneral.h | 2 + src/PySTARE.cpp | 4 + src/PySTARE.i | 34 +------ src/STARE.C | 8 +- src/SpatialIndex.cpp | 132 +++++++++++++++++++++++++-- tests/CUTE/STARE_test.cpp | 155 ++++++++++++++++++++++++++++++++ 9 files changed, 459 insertions(+), 40 deletions(-) create mode 100644 contrib/python/.gitignore create mode 100644 contrib/python/example-viz-1.py diff --git a/contrib/python/.gitignore b/contrib/python/.gitignore new file mode 100644 index 0000000..facc968 --- /dev/null +++ b/contrib/python/.gitignore @@ -0,0 +1 @@ +example-track-1.py diff --git a/contrib/python/example-viz-1.py b/contrib/python/example-viz-1.py new file mode 100644 index 0000000..845fce5 --- /dev/null +++ b/contrib/python/example-viz-1.py @@ -0,0 +1,151 @@ + +from math import ceil +import pystare as ps + +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import cartopy.crs as ccrs + +import numpy as np + +def shiftarg_lon(lon): + if(lon>180): + return ((lon + 180.0) % 360.0)-180.0 + else: + return lon + +def triangulate(i0,i1,i2): + # i0,i1,i2,ic = ps.to_vertices(indices) + i0lat,i0lon = ps.to_latlon(i0) + i1lat,i1lon = ps.to_latlon(i1) + i2lat,i2lon = ps.to_latlon(i2) + lats = np.zeros([3*len(i0lat)],dtype=np.double) + lons = np.zeros([3*len(i0lat)],dtype=np.double) + intmat = [] + k=0 + for i in range(len(i0)): + lats[k] = i0lat[i] + lons[k] = i0lon[i] + lats[k+1] = i1lat[i] + lons[k+1] = i1lon[i] + lats[k+2] = i2lat[i] + lons[k+2] = i2lon[i] + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + return lons,lats,intmat + +# lat = np.array([0, 0,60], dtype=np.double) +# lon = np.array([0,60,30], dtype=np.double) + +lat = np.array([0, 0,60], dtype=np.double) +lon = np.array([0,60,30], dtype=np.double) + +s_resolution = 6 +indices = ps.from_latlon(lat, lon, s_resolution) + +def test1(indices): + "Plot the outline of a triangle using first three indices." + i0 = np.array([indices[0]], dtype=np.int64) + i1 = np.array([indices[1]], dtype=np.int64) + i2 = np.array([indices[2]], dtype=np.int64) + lons,lats,intmat = triangulate(i0,i1,i2) + # print('lats: ',lats) + # print('lons: ',lons) + triang = tri.Triangulation(lons,lats,intmat) + return triang + +def plot1(triang): + ax = plt.axes(projection=ccrs.PlateCarree()) + ax.set_xlim(-180,180) + ax.set_ylim(-90,90) + ax.coastlines() + # plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) + # plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) + # plt.triplot(triang,'ko-') + plt.triplot(triang,'r-') + # plt.show() + return + +plot1(test1(indices)) + +level = 3 +hull = ps.to_hull_range(indices,level,100) +h0,h1,h2,hc = ps.to_vertices(hull) +lons1,lats1,intmat1 = triangulate(h0,h1,h2) +print('0 hull len: ',len(hull)) +print('0 hull lats len: ',len(lats1)) +jtest=3 +j=jtest +print('0 hull lats1: ',j,lats1[j*3:(j+1)*3]) +print('0 hull lons1: ',j,lons1[j*3:(j+1)*3]) +j=0 +print('0 hull lats1: ',[i for i in lats1[j:j+12]]) +print('0 hull lons1: ',[i for i in lons1[j:j+12]]) +print('') + +tid = np.array([0x4c0000000000003],dtype=np.int64) +t0,t1,t2,tc = ps.to_vertices(tid) +print('t0: ',hex(t0[0])) +print('t1: ',hex(t1[0])) +print('t2: ',hex(t2[0])) +print('tc: ',hex(tc[0])) +print('t0 ll : ',ps.to_latlon(t0)) +print('t1 ll : ',ps.to_latlon(t1)) +print('t2 ll : ',ps.to_latlon(t2)) +print('tc ll : ',ps.to_latlon(tc)) + +print('') + +if True: + # i=9; ilen=2 + # i=10; ilen=1 + # i=jtest; ilen=4 + i=jtest; ilen=1 + id_test = np.array(hull[i:i+ilen],dtype=np.int64) + print('i,id : ',i,[hex(j) for j in id_test]) + i0=i*3; i1=(i+ilen)*3 + lats1 = lats1[i0:i1] + lons1 = lons1[i0:i1] + # intmat1 = [intmat1[i]] + intmat1 = [] + for j in range(ilen): + intmat1.append([3*j,3*j+1,3*j+2]) + # intmat1 = [[0,2,1]] + i0test,i1test,i2test,ictest = ps.to_vertices(id_test) + print('test id: ',[hex(i) for i in id_test]) + print('test ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) + print('test im : ',[i for i in intmat1]) + lonstest,latstest,intmattest = triangulate(i0test,i1test,i2test) + print('test lat: ',latstest) + print('test lon: ',lonstest) + triangtest = tri.Triangulation(lonstest,latstest,intmattest) + plt.triplot(triangtest,'g-') + plt.scatter(lonstest,latstest,s=50,c='g') + +print('level : ',level) +# print('hull : ',[hex(i) for i in hull]) +print('hull len: ',len(hull)) +print('hull ll len: ',len(lons1)) +print('hull im len: ',len(intmat1)) +print('hull ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) +print('hull im : ',[i for i in intmat1]) + +triang1 = tri.Triangulation(lons1,lats1,intmat1) +plt.triplot(triang1,'b-') + +# ax = plt.axes(projection=ccrs.PlateCarree()) +# ax.set_xlim(-180,180) +# ax.set_ylim(-90,90) +# ax.coastlines() +# # plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) +# # plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) +# # plt.triplot(triang,'ko-') +# plt.triplot(triang1,'r-') +plt.show() + + + + diff --git a/contrib/python/test.py b/contrib/python/test.py index ca9cf97..a95ea63 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -128,3 +128,15 @@ cmp1 = pystare.cmp_spatial(indices,indices1) print('cmp1: ',cmp1) print('') + +tid = numpy.array([0x4c0000000000003],dtype=numpy.int64) +t0,t1,t2,tc = pystare.to_vertices(tid) +print('t0: ',hex(t0[0])) +print('t1: ',hex(t1[0])) +print('t2: ',hex(t2[0])) +print('tc: ',hex(tc[0])) +print('t0 ll : ',pystare.to_latlon(t0)) +print('t1 ll : ',pystare.to_latlon(t1)) +print('t2 ll : ',pystare.to_latlon(t2)) +print('tc ll : ',pystare.to_latlon(tc)) +print('') diff --git a/include/SpatialGeneral.h b/include/SpatialGeneral.h index fe8644d..1f952ab 100644 --- a/include/SpatialGeneral.h +++ b/include/SpatialGeneral.h @@ -115,6 +115,8 @@ const float64 gPi = 3.1415926535897932385E0 ; const float64 gPio2 = 3.1415926535897932385E0/2.0 ; const float64 gPr = gPi/180.0; const float64 piDiv180 = gPr; + +// TODO Set these by environment variable or other property list. // const float64 gEpsilon = 1.0E-15; // const float64 gEpsilon = 1.0E-16; // const float64 gEpsilon = 1.0E-17; diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index ae157db..633cf5d 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -127,9 +127,13 @@ void _to_hull_range(int64_t* indices, int len, int resolution, int64_t* range_in cout << "_to_hull_range-warning: range_indices.size = " << len_ri << " too small." << endl << flush; cout << "_to_hull_range-warning: result size = " << result.size() << "." << endl << flush; } + int k=10; + cout << "thr "; for(int i=0; i < (len_ri < result.size() ? len_ri : result.size()); ++i) { + if(k-->0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } range_indices[i] = result[i]; } + cout << dec << endl << flush; result_size[0] = result.size(); } diff --git a/src/PySTARE.i b/src/PySTARE.i index 7a9eb21..1b6c308 100644 --- a/src/PySTARE.i +++ b/src/PySTARE.i @@ -185,8 +185,10 @@ $5 = (int*) array_data((PyArrayObject*)out3); } -/* maps ONE int64_t input array to THREE 1D int output array with the same length */ -/* We use this to convert STARE index to triangle vertices & centroid */ +/* + * maps ONE int64_t input array to THREE 1D int output array with the same length + * We use this to convert STARE index to triangle vertices & centroid + */ %typemap(in, numinputs=1) (int64_t* in_array, int length, int64_t* out_array1, int64_t* out_array2, int64_t* out_array3, int64_t* out_array4) (PyObject* out1=NULL, PyObject* out2=NULL, PyObject* out3=NULL, PyObject* out4=NULL) @@ -215,30 +217,6 @@ $6 = (int64_t*) array_data((PyArrayObject*)out4); } -/* maps TWO int64_t input arrays to ONE 1D int64_t output array with the cross-product length */ -/* We use this to compare two lists of spatial ids, i.e. a cross-join of sorts. */ -%typemap(in, numinputs=1) - (int64_t* in_array1, int length1, int64_t* in_array2, int length2, int64_t* out_array, int out_length) - (PyObject* out=NULL) -{ - int is_new_object=0; - npy_intp size[1] = { -1}; - PyArrayObject* array = obj_to_array_contiguous_allow_conversion($input, NPY_INT64, &is_new_object); - if (!array || !require_dimensions(array, 1)) SWIG_fail; - - size[0] = PyArray_DIM(array, 0); - - out1 = PyArray_SimpleNew(1, size, NPY_INT64); - if (!out1) SWIG_fail; - - $1 = (int64_t*) array_data(array); - $2 = (int) array_size(array,0); - $3 = (int64_t*) array_data(array); - $4 = (int) array_size(array,0); - $5 = (int64_t*) array_data((PyArrayObject*)out1); - $6 = (int) array_size(array,0); -} - /****************/ /* OUT typemaps */ /****************/ @@ -365,10 +343,6 @@ (int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ) } -# %apply (int64_t* in_array1, int length1, int64_t* in_array2, int length2, int64_t* out_array, int out_length) { -# (int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* cmp, int len12) -# } - %pythonprepend from_utc(int64_t*, int, int64_t*, int) %{ import numpy datetime = datetime.astype(numpy.int64) diff --git a/src/STARE.C b/src/STARE.C index a0979ea..a993604 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -164,6 +164,7 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati // cout << "sid: " << spatialStareId << endl << flush; uint64 htmID = htmIDFromValue(spatialStareId); + cout << "lldfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; SpatialVector v; /// This returns the center of the triangle (at index.search_level). Need to extract the position information. @@ -172,11 +173,13 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati float64 lat=-999, lon=-999; v.getLatLonDegrees(lat, lon); - // cout << "sid-latlon: " << lat << ", " << lon << endl << flush; + cout << "0 sid-latlon: " << lat << ", " << lon << endl << flush; // LatLonDegrees64 latlon = {.lat = lat, .lon = lon }; LatLonDegrees64 latlon(lat, lon); // = {.lat = lat, .lon = lon }; + cout << "1 sid-latlon: " << latlon.lat << ", " << latlon.lon << endl << flush; + // return latlon; return LatLonDegrees64(lat, lon); @@ -191,6 +194,7 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati SpatialVector STARE::SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId) { // uint64 htmID = htmIDFromValue(spatialStareId,STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); // Max resolution uint64 htmID = htmIDFromValue(spatialStareId); // Max resolution + cout << "svfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; SpatialVector v; /// This returns the center of the triangle (at index.search_level). Need to extract the position information. sIndex.pointByHtmId(v, htmID); @@ -506,7 +510,7 @@ STARE_SpatialIntervals STARE::ConvexHull(LatLonDegrees64ValueVector points,int f STARE_SpatialIntervals STARE::ConvexHull(STARE_ArrayIndexSpatialValues points,int force_resolution_level) { LatLonDegrees64ValueVector latlon; - for(auto i=points.begin(); i != points.end(); ++i) { + for( STARE_ArrayIndexSpatialValues::iterator i=points.begin(); i != points.end(); ++i) { latlon.push_back(LatLonDegreesFromValue(*i)); } return ConvexHull(latlon,force_resolution_level); diff --git a/src/SpatialIndex.cpp b/src/SpatialIndex.cpp index f93db8f..24e1c9c 100644 --- a/src/SpatialIndex.cpp +++ b/src/SpatialIndex.cpp @@ -1334,24 +1334,124 @@ SpatialIndex::idByPoint(SpatialVector & v) const { uint64 index; uint64 ID; + bool verbose = false; + bool nudge = false; + + /**/ // start with the 8 root triangles, find the one which v points to for(index=1; index <=8; index++) { - if( (V(0) ^ V(1)) * v < -gEpsilon) continue; - if( (V(1) ^ V(2)) * v < -gEpsilon) continue; - if( (V(2) ^ V(0)) * v < -gEpsilon) continue; - break; + if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; } + cout << "0 preamble i: " << index << endl << flush; + nudge = index == 9; // loop through matching child until leaves are reached while(ICHILD(0)!=0) { uint64 oldindex = index; for(size_t i = 0; i < 4; i++) { index = nodes_[oldindex].childID_[i]; - if( (V(0) ^ V(1)) * v < -gEpsilon) continue; - if( (V(1) ^ V(2)) * v < -gEpsilon) continue; - if( (V(2) ^ V(0)) * v < -gEpsilon) continue; - break; + if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; } } + cout << "1 preamble i: " << index << endl << flush; + /**/ + + // if(index == 9) { + + /* THE OLD WAY */ +// Diagnostics + SpatialVector dc; + + float64 dcs[8]; + for(index=1; index <=8; index++) { + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + dcs[index-1] = dc.length(); + cout << dec << index << " index,dc " << setprecision(16) << dcs[index-1] << endl << flush; + } + + int index_dcs_sort[8]; + index_dcs_sort[0]=0; + + for(int i = 1; i < 8; ++i) { + int icmp = 0; + while( icmp < i && (dcs[i] >= dcs[index_dcs_sort[icmp]]) ) {++icmp; } + if( icmp <= i) { + for( int j=i+1; j>=icmp; --j ) { + index_dcs_sort[j]=index_dcs_sort[j-1]; + } + } + index_dcs_sort[icmp] = i; + } + for(int i=0;i<8;++i) { + cout << i << " i,sort,dcs " << index_dcs_sort[i] << " " << dcs[index_dcs_sort[i]] << endl << flush; + } + + // start with the 8 root triangles, find the one which v points to + for(index=1; index <=8; index++) { + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } + + if(nudge) { + float64 epsilon = 1.0e-13; + SpatialVector ctr = V(0)+V(1)+V(2); ctr.normalize(); + v = (1-epsilon)*v + epsilon*ctr; v.normalize(); + for(index=1; index <=8; index++) { + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } + } + + cout << dec << "index " << index << endl << flush; + float64 dc_start, dc_end; + float64 dc_improvement = 2.0; + int attempt = 0; + int index_tried = index; + int itry = 0; + do { + ++attempt; + if(attempt == 9) { + throw SpatialFailure("SpatialIndex::idByPoint(sv): Lost Point Failure 1."); + } + if(attempt>1) { + index = index_dcs_sort[itry]+1; + if( index == index_tried ) { + ++itry; + if(itry == 8) { + throw SpatialFailure("SpatialIndex::idByPoint(sv): Lost Point Failure 2."); + } + index = index_dcs_sort[itry]+1; + } + ++itry; + index_tried = index; + } + dc_start = dcs[index-1]; + cout << attempt << " trying " << index << endl << flush; + // loop through matching child until leaves are reached + int k = 1; + while(ICHILD(0)!=0) { + uint64 oldindex = index; + for(size_t i = 0; i < 4; i++) { + index = nodes_[oldindex].childID_[i]; + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + cout << dec << k++ << " k,dc " << setprecision(16) << dc.length() << " i: " << index << endl << flush; + } + dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; + dc_end = dc.length(); + dc_improvement = dc_end/dc_start; + cout << "dc start, end, ratio: " << dc_start << " " << dc_end << " " << dc_improvement << endl << flush; + } while ( dc_improvement > 0.25 && dc_end > 0.25 ); + // } + /**/ + // what if we haven't gotten to build level? // cout << "idbp: 100 " << endl << flush; // return if we have reached maxlevel @@ -1386,9 +1486,18 @@ SpatialIndex::idByPoint(SpatialVector & v) const { // cout << "idbp: 300 maxlevel_-buildlevel_ = " << level << endl << flush; // TODO make this whole routine less ad-hoc + SpatialVector delta; + cout << endl << flush; if(level>0) { int level0 = buildlevel_; while(level--) { + cout << setprecision(16) << dec; + // cout << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush; + cout << level << " level,d0,d1,d2 "; + delta = v-v0; cout << delta.length() << " "; + delta = v-v1; cout << delta.length() << " "; + delta = v-v2; cout << delta.length() << " "; + cout << endl << flush; int subTriangleIndex = subTriangleIndexByPoint(v,v0,v1,v2); // cout << "idbp: 350 level = " << level << ", level0 = " << level0++ << ", subTriangleIndex = " << subTriangleIndex << ", name = " << name; name[len++] = '0'+char(subTriangleIndex); @@ -1400,6 +1509,13 @@ SpatialIndex::idByPoint(SpatialVector & v) const { throw SpatialFailure("SpatialIndex::idByPoint::LevelTooLow len < 2"); } } + cout << level << " level,d0,d1,d2 "; + delta = v-v0; cout << delta.length() << " "; + delta = v-v1; cout << delta.length() << " "; + delta = v-v2; cout << delta.length() << " "; + cout << endl << flush; + cout << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush; + name[len] = '\0'; // cout << "idbp: 400 name = '" << name << "', len = " << len << endl << flush; diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index be3cef6..6e00401 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -1089,6 +1089,161 @@ void STARE_test() { ASSERT_EQUALM("iv-to-v-roundtrip",siv,(siv1 & ~lj.levelMaskSciDB) | 8); } + if(true) { + STARE index; + SpatialVector v(1,1,1); v.normalize(); + STARE_ArrayIndexSpatialValue sid = index.ValueFromSpatialVector(v); + SpatialVector vr = index.SpatialVectorFromValue(sid); + STARE_ArrayIndexSpatialValue sidr = index.ValueFromSpatialVector(vr); + cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec << endl << flush; + cout << "sidr 0x" << setw(16) << setfill('0') << hex << sidr << dec << endl << flush; + cout << "v " << v << endl << flush; + cout << "vr " << vr << endl << flush; + cout << endl << flush; + } + + if(true) { + /* + i,id : 3 ['0x4c0000000000003'] + test id: ['0x4c0000000000003'] + test ll : [(14.255438319990454, 64.21872912284311), (-22.58410062366997, 70.08873514750778), (0.8929126045078797, 67.65325299601919)] + test im : [[0, 1, 2]] + test lat: [ 14.25543832 -22.58410062 0.8929126 ] + test lon: [64.21872912 70.08873515 67.653253 ] + */ + STARE index; + // STARE_ArrayIndexSpatialValue sid = 0x4c0000000000003; + STARE_ArrayIndexSpatialValue sid = 0x4c000000000001b; + STARE_ArrayIndexSpatialValue sidtmp; + SpatialVector v_sid = index.SpatialVectorFromValue(sid); + LatLonDegrees64 v_ll = index.LatLonDegreesFromValue(sid); + float64 v_lat,v_lon; + v_sid.getLatLonDegrees(v_lat, v_lon); + + cout << "v_sid " << v_sid << endl << flush; + cout << "v_sid ll " << v_lat << " " << v_lon << endl << flush; + cout << "v idx ll " << v_ll.lat << " " << v_ll.lon << endl << flush; + cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec << endl << flush; + cout << endl << flush; + + cout << "Triangle " << endl << flush; + Triangle tr = index.TriangleFromValue(sid); + for( int i=0; i<3; ++i) { + float64 lt,ln; + tr.vertices[i].getLatLonDegrees(lt,ln); + cout << setprecision(16); + cout << i << " i,v " << tr.vertices[i] << " -- " << lt << "," << ln << endl << flush; + } + cout << endl << flush; + + cout << "index.ValueFromSpatialVector triangle" << endl << flush; + for( int i=0; i<3; ++i) { + sidtmp = index.ValueFromSpatialVector(tr.vertices[i],27); + SpatialVector vtmp = index.SpatialVectorFromValue(sidtmp); + float64 lt,ln; + vtmp.getLatLonDegrees(lt,ln); + cout << i << " i,vt " << vtmp << " -- " << lt << "," << ln << endl << flush; + } + cout << endl << flush; + + SpatialVector vecs[3] = { + SpatialVector( 0.4619397639595428, 0.8446232009168338, 0.2705980468259219 ), + SpatialVector( 0.4619397725326853, 0.8446231986774957, 0.2705980391803436 ), + SpatialVector( 0.4619397591445266, 0.8446232068078616, 0.2705980366578833 ) + }; + + cout << "index... from Explicit vecs" << endl << flush; + for( int i=0; i<3; ++i) { + sidtmp = index.ValueFromSpatialVector(vecs[i]); + SpatialVector vtmp = index.SpatialVectorFromValue(sidtmp); + float64 lt,ln; + vtmp.getLatLonDegrees(lt,ln); + cout << i << " i,vt " << vtmp << " -- " << lt << "," << ln << endl << flush; + } + cout << endl << flush; + + + for( int i=0; i<3; ++i) { + float64 lat, lon; + SpatialVector vtmp = tr.vertices[i]; + vtmp.getLatLonDegrees(lat,lon); + // tr.vertices[i].getLatLonDegrees(lat,lon); + sidtmp = index.ValueFromSpatialVector(vtmp,27); + // sidtmp = index.ValueFromSpatialVector(vtmp); + // sidtmp = index.ValueFromSpatialVector(tr.vertices[i]); + LatLonDegrees64 latlontmp = index.LatLonDegreesFromValue(sidtmp); + cout << i << " i,tr latlon " << lat << " " << lon << " -- " + << latlontmp.lat << " " << latlontmp.lon + << " -- " + << "0x" << setw(16) << setfill('0') << hex << sidtmp << dec + << endl << flush; + } + cout << endl << flush; + + cout << setprecision(16); cout << "tr.v[1] " << tr.vertices[1] << endl << flush; + SpatialVector v = tr.vertices[1]; + sid = index.ValueFromSpatialVector(v); + SpatialVector vr = index.SpatialVectorFromValue(sid); + STARE_ArrayIndexSpatialValue sidr = index.ValueFromSpatialVector(vr); + float64 lat,lon; + LatLonDegrees64 ll = index.LatLonDegreesFromValue(sid); + LatLonDegrees64 llr = index.LatLonDegreesFromValue(sidr); + cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec + << " - " << ll.lat << " " << ll.lon + << endl << flush; + cout << "sidr 0x" << setw(16) << setfill('0') << hex << sidr << dec + << " - " << llr.lat << " " << llr.lon + << endl << flush; + v.getLatLonDegrees(lat, lon); + cout << setprecision(16); + cout << "tr.v[1] " << tr.vertices[1] << endl << flush; + cout << "v " << v + << " - " << lat << " " << lon + << endl << flush; + vr.getLatLonDegrees(lat, lon); + cout << "vr " << vr + << " - " << lat << " " << lon + << endl << flush; + cout << endl << flush; + + } + + if(true) { + cout << endl << "------" << endl << endl << flush; + STARE index; + int resolution = 27; + SpatialIndex sIndex = index.getIndex(resolution); + for( int i = 0; i < 11; ++i ) { + SpatialVector v0 ( 0.4619397639595428, 0.8446232009168338, 0.2705980468259219 ); + // SpatialVector delta = SpatialVector(0.0,0.0,(1-i*0.1)*1.0e-16); + SpatialVector delta = SpatialVector(0.0,0.0,(1-i*0.1)*1.0e-15); + cout << i << " i,delta " << delta << endl << flush; + // v0 = v0 + SpatialVector(0.0,0.0,0.9e-8); + v0 = v0 - delta; + v0.normalize(); + // SpatialVector v0 ( 1, 0, 0); // Works for this one. + STARE_ArrayIndexSpatialValue siv0 = index.ValueFromSpatialVector(v0,resolution); + BitShiftNameEncoding rightJustified(sIndex.idByPoint(v0)); + cout << "sI idbp 0x" << setw(16) << setfill('0') << hex << sIndex.idByPoint(v0) << dec << endl << flush; + cout << "rj name "<< rightJustified.getName() << endl << flush; + + EmbeddedLevelNameEncoding leftJustified(rightJustified.leftJustifiedId()); + EmbeddedLevelNameEncoding leftJustifiedWithResolution = leftJustified.atLevel(resolution, true); // True means keep all bits + STARE_ArrayIndexSpatialValue siv = leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); + cout << "siv0 " << setw(16) << setfill('0') << hex << siv0 << endl << flush; + cout << "siv " << setw(16) << setfill('0') << hex << siv << endl << flush; + SpatialVector v_siv0 = index.SpatialVectorFromValue(siv0); + SpatialVector v_siv = index.SpatialVectorFromValue(siv); + cout << setprecision(16); + cout << "v0 " << v0 << endl << flush; + cout << "v_siv0 " << v_siv0 << endl << flush; + cout << "v_siv " << v_siv << endl << flush; + cout << endl << "--" << endl << endl << flush; + } + + + } + // FAIL(); } From 5d8c9b66c05b4bf0dc590a46d1e3b35e4c76fb5b Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Fri, 20 Sep 2019 09:52:57 -0400 Subject: [PATCH 24/31] Made idByPoint more expensive, but more robust. --- src/STARE.C | 8 +-- src/SpatialIndex.cpp | 86 +++++++++++++++++++------------- tests/CUTE/STARE_test.cpp | 49 ++++++++++++++---- tests/CUTE/SpatialIndex_test.cpp | 7 ++- 4 files changed, 102 insertions(+), 48 deletions(-) diff --git a/src/STARE.C b/src/STARE.C index a993604..842cd10 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -164,7 +164,7 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati // cout << "sid: " << spatialStareId << endl << flush; uint64 htmID = htmIDFromValue(spatialStareId); - cout << "lldfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; + // cout << "lldfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; SpatialVector v; /// This returns the center of the triangle (at index.search_level). Need to extract the position information. @@ -173,12 +173,12 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati float64 lat=-999, lon=-999; v.getLatLonDegrees(lat, lon); - cout << "0 sid-latlon: " << lat << ", " << lon << endl << flush; + // cout << "0 sid-latlon: " << lat << ", " << lon << endl << flush; // LatLonDegrees64 latlon = {.lat = lat, .lon = lon }; LatLonDegrees64 latlon(lat, lon); // = {.lat = lat, .lon = lon }; - cout << "1 sid-latlon: " << latlon.lat << ", " << latlon.lon << endl << flush; + // cout << "1 sid-latlon: " << latlon.lat << ", " << latlon.lon << endl << flush; // return latlon; return LatLonDegrees64(lat, lon); @@ -194,7 +194,7 @@ LatLonDegrees64 STARE::LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spati SpatialVector STARE::SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId) { // uint64 htmID = htmIDFromValue(spatialStareId,STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); // Max resolution uint64 htmID = htmIDFromValue(spatialStareId); // Max resolution - cout << "svfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; + // cout << "svfv htmID " << setw(16) << setfill('0') << hex << htmID << dec << endl << flush; SpatialVector v; /// This returns the center of the triangle (at index.search_level). Need to extract the position information. sIndex.pointByHtmId(v, htmID); diff --git a/src/SpatialIndex.cpp b/src/SpatialIndex.cpp index 24e1c9c..def40b0 100644 --- a/src/SpatialIndex.cpp +++ b/src/SpatialIndex.cpp @@ -1314,6 +1314,15 @@ SpatialIndex::NeighborsAcrossVerticesFromEdges( workspace[jw++] = q8; } +#ifndef DIAG +#define DIAG1(expr) +#define DIAGOUT2(out,expr) +#define DIAGOUTDELTA(out,a,b) +#else +#define DIAG1(expr) expr; +#define DIAGOUT2(out,expr) out << expr; +#define DIAGOUTDELTA(out,a,b) {SpatialVector delta_ = a-b; cout << delta_.length() << " ";} +#endif //////////////////IDBYPOINT//////////////////////////////////////////////// /** Find a leaf node where a vector points to @@ -1342,8 +1351,9 @@ SpatialIndex::idByPoint(SpatialVector & v) const { for(index=1; index <=8; index++) { if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; } - cout << "0 preamble i: " << index << endl << flush; + DIAGOUT2(cout,"0 preamble i: " << index << endl << flush;); nudge = index == 9; + // loop through matching child until leaves are reached while(ICHILD(0)!=0) { uint64 oldindex = index; @@ -1352,20 +1362,23 @@ SpatialIndex::idByPoint(SpatialVector & v) const { if(isInsideBarycentric(v,V(0),V(1),V(2),verbose)) break; } } - cout << "1 preamble i: " << index << endl << flush; + DIAGOUT2(cout,"1 preamble i: " << index << endl << flush;); + // int index_barycentric = index; /**/ // if(index == 9) { - /* THE OLD WAY */ -// Diagnostics + /* THE OLD WAY + * TODO Why keep the legacy point-finding behavior? Right thing to do would be to examine where barycentric and old way differ. + */ + SpatialVector dc; float64 dcs[8]; for(index=1; index <=8; index++) { dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; dcs[index-1] = dc.length(); - cout << dec << index << " index,dc " << setprecision(16) << dcs[index-1] << endl << flush; + DIAGOUT2(cout,dec << index << " index,dc " << setprecision(16) << dcs[index-1] << endl << flush;); } int index_dcs_sort[8]; @@ -1381,9 +1394,7 @@ SpatialIndex::idByPoint(SpatialVector & v) const { } index_dcs_sort[icmp] = i; } - for(int i=0;i<8;++i) { - cout << i << " i,sort,dcs " << index_dcs_sort[i] << " " << dcs[index_dcs_sort[i]] << endl << flush; - } + DIAG1(for(int i=0;i<8;++i) {DIAGOUT2(cout,i << " i,sort,dcs " << index_dcs_sort[i] << " " << dcs[index_dcs_sort[i]] << endl << flush;);}); // start with the 8 root triangles, find the one which v points to for(index=1; index <=8; index++) { @@ -1394,18 +1405,20 @@ SpatialIndex::idByPoint(SpatialVector & v) const { } if(nudge) { - float64 epsilon = 1.0e-13; - SpatialVector ctr = V(0)+V(1)+V(2); ctr.normalize(); - v = (1-epsilon)*v + epsilon*ctr; v.normalize(); - for(index=1; index <=8; index++) { - if( (V(0) ^ V(1)) * v < -gEpsilon) continue; - if( (V(1) ^ V(2)) * v < -gEpsilon) continue; - if( (V(2) ^ V(0)) * v < -gEpsilon) continue; - break; - } + // if(false) { + float64 epsilon = 1.0e-13; // TODO inject this dependency, expose as environment variable + // float64 epsilon = 1.0e-17; // TODO inject this dependency, expose as environment variable + SpatialVector ctr = V(0)+V(1)+V(2); ctr.normalize(); + v = (1-epsilon)*v + epsilon*ctr; v.normalize(); + for(index=1; index <=8; index++) { + if( (V(0) ^ V(1)) * v < -gEpsilon) continue; + if( (V(1) ^ V(2)) * v < -gEpsilon) continue; + if( (V(2) ^ V(0)) * v < -gEpsilon) continue; + break; + } } - cout << dec << "index " << index << endl << flush; + DIAGOUT2(cout,dec << "index " << index << endl << flush;); float64 dc_start, dc_end; float64 dc_improvement = 2.0; int attempt = 0; @@ -1429,7 +1442,7 @@ SpatialIndex::idByPoint(SpatialVector & v) const { index_tried = index; } dc_start = dcs[index-1]; - cout << attempt << " trying " << index << endl << flush; + DIAG1(cout << attempt << " trying " << index << endl << flush;); // loop through matching child until leaves are reached int k = 1; while(ICHILD(0)!=0) { @@ -1442,13 +1455,14 @@ SpatialIndex::idByPoint(SpatialVector & v) const { break; } dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; - cout << dec << k++ << " k,dc " << setprecision(16) << dc.length() << " i: " << index << endl << flush; + DIAG1(cout << dec << k++ << " k,dc " << setprecision(16) << dc.length() << " i: " << index << endl << flush;); } dc = V(0)+V(1)+V(2); dc.normalize(); dc = dc - v; dc_end = dc.length(); dc_improvement = dc_end/dc_start; - cout << "dc start, end, ratio: " << dc_start << " " << dc_end << " " << dc_improvement << endl << flush; + DIAG1(cout << "dc start, end, ratio: " << dc_start << " " << dc_end << " " << dc_improvement << endl << flush;); } while ( dc_improvement > 0.25 && dc_end > 0.25 ); + // } /**/ @@ -1487,17 +1501,17 @@ SpatialIndex::idByPoint(SpatialVector & v) const { // TODO make this whole routine less ad-hoc SpatialVector delta; - cout << endl << flush; + DIAGOUT2(cout,endl << flush); if(level>0) { int level0 = buildlevel_; while(level--) { - cout << setprecision(16) << dec; + DIAGOUT2(cout,setprecision(16) << dec;); // cout << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush; - cout << level << " level,d0,d1,d2 "; - delta = v-v0; cout << delta.length() << " "; - delta = v-v1; cout << delta.length() << " "; - delta = v-v2; cout << delta.length() << " "; - cout << endl << flush; + DIAGOUT2(cout,level << " level,d0,d1,d2 ";); + DIAGOUTDELTA(cout,v,v0); + DIAGOUTDELTA(cout,v,v1); + DIAGOUTDELTA(cout,v,v2); + DIAGOUT2(cout,endl << flush;); int subTriangleIndex = subTriangleIndexByPoint(v,v0,v1,v2); // cout << "idbp: 350 level = " << level << ", level0 = " << level0++ << ", subTriangleIndex = " << subTriangleIndex << ", name = " << name; name[len++] = '0'+char(subTriangleIndex); @@ -1509,12 +1523,12 @@ SpatialIndex::idByPoint(SpatialVector & v) const { throw SpatialFailure("SpatialIndex::idByPoint::LevelTooLow len < 2"); } } - cout << level << " level,d0,d1,d2 "; - delta = v-v0; cout << delta.length() << " "; - delta = v-v1; cout << delta.length() << " "; - delta = v-v2; cout << delta.length() << " "; - cout << endl << flush; - cout << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush; + + DIAGOUT2(cout,level << " level,d0,d1,d2 ";); + DIAGOUTDELTA(cout,v,v0); + DIAGOUTDELTA(cout,v,v1); + DIAGOUTDELTA(cout,v,v2); + DIAGOUT2(cout,endl << flush << level << " level,v... " << v << " " << v0 << " " << v1 << " " << v2 << endl << flush;); name[len] = '\0'; @@ -1538,6 +1552,10 @@ SpatialIndex::idByPoint(SpatialVector & v) const { return ID; } +#undef DIAG +#undef DIAG1 +#undef DIAGOUT2 +#undef DIAGOUTDELTA uint64 SpatialIndex::indexAtNodeIndex(uint64 nodeIndex) { return nodes_[nodeIndex].index_; diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index 6e00401..9f8437f 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -12,6 +12,16 @@ #include "SpatialInterface.h" +#ifndef DIAG +#define DIAG1(expr) +#define DIAGOUT2(out,expr) +#define DIAGOUTDELTA(out,a,b) +#else +#define DIAG1(expr) expr; +#define DIAGOUT2(out,expr) out << expr; +#define DIAGOUTDELTA(out,a,b) {SpatialVector delta_ = a-b; cout << delta_.length() << " ";} +#endif + void STARE_test() { SpatialVector @@ -1095,11 +1105,14 @@ void STARE_test() { STARE_ArrayIndexSpatialValue sid = index.ValueFromSpatialVector(v); SpatialVector vr = index.SpatialVectorFromValue(sid); STARE_ArrayIndexSpatialValue sidr = index.ValueFromSpatialVector(vr); +#ifdef DIAG cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec << endl << flush; cout << "sidr 0x" << setw(16) << setfill('0') << hex << sidr << dec << endl << flush; cout << "v " << v << endl << flush; cout << "vr " << vr << endl << flush; cout << endl << flush; +#endif + ASSERT_EQUAL(sid,sidr); } if(true) { @@ -1120,14 +1133,21 @@ void STARE_test() { float64 v_lat,v_lon; v_sid.getLatLonDegrees(v_lat, v_lon); +#ifdef DIAG cout << "v_sid " << v_sid << endl << flush; cout << "v_sid ll " << v_lat << " " << v_lon << endl << flush; cout << "v idx ll " << v_ll.lat << " " << v_ll.lon << endl << flush; cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec << endl << flush; cout << endl << flush; +#endif + + ASSERT_EQUAL(v_lat,v_ll.lat); + ASSERT_EQUAL(v_lon,v_ll.lon); - cout << "Triangle " << endl << flush; Triangle tr = index.TriangleFromValue(sid); + +#ifdef DIAG + cout << "Triangle " << endl << flush; for( int i=0; i<3; ++i) { float64 lt,ln; tr.vertices[i].getLatLonDegrees(lt,ln); @@ -1145,6 +1165,7 @@ void STARE_test() { cout << i << " i,vt " << vtmp << " -- " << lt << "," << ln << endl << flush; } cout << endl << flush; +#endif SpatialVector vecs[3] = { SpatialVector( 0.4619397639595428, 0.8446232009168338, 0.2705980468259219 ), @@ -1152,6 +1173,7 @@ void STARE_test() { SpatialVector( 0.4619397591445266, 0.8446232068078616, 0.2705980366578833 ) }; +#ifdef DIAG cout << "index... from Explicit vecs" << endl << flush; for( int i=0; i<3; ++i) { sidtmp = index.ValueFromSpatialVector(vecs[i]); @@ -1162,7 +1184,6 @@ void STARE_test() { } cout << endl << flush; - for( int i=0; i<3; ++i) { float64 lat, lon; SpatialVector vtmp = tr.vertices[i]; @@ -1179,8 +1200,10 @@ void STARE_test() { << endl << flush; } cout << endl << flush; +#endif - cout << setprecision(16); cout << "tr.v[1] " << tr.vertices[1] << endl << flush; + + DIAGOUT2(cout,setprecision(16); cout << "tr.v[1] " << tr.vertices[1] << endl << flush;); SpatialVector v = tr.vertices[1]; sid = index.ValueFromSpatialVector(v); SpatialVector vr = index.SpatialVectorFromValue(sid); @@ -1188,6 +1211,7 @@ void STARE_test() { float64 lat,lon; LatLonDegrees64 ll = index.LatLonDegreesFromValue(sid); LatLonDegrees64 llr = index.LatLonDegreesFromValue(sidr); +#ifdef DIAG cout << "sid 0x" << setw(16) << setfill('0') << hex << sid << dec << " - " << ll.lat << " " << ll.lon << endl << flush; @@ -1205,11 +1229,14 @@ void STARE_test() { << " - " << lat << " " << lon << endl << flush; cout << endl << flush; +#endif + ASSERT_EQUAL(sid,sidr); + ASSERT_EQUALDM("v-to-vr",v,vr,1.0e-8); } if(true) { - cout << endl << "------" << endl << endl << flush; + DIAGOUT2(cout,endl << "------" << endl << endl << flush;); STARE index; int resolution = 27; SpatialIndex sIndex = index.getIndex(resolution); @@ -1217,31 +1244,35 @@ void STARE_test() { SpatialVector v0 ( 0.4619397639595428, 0.8446232009168338, 0.2705980468259219 ); // SpatialVector delta = SpatialVector(0.0,0.0,(1-i*0.1)*1.0e-16); SpatialVector delta = SpatialVector(0.0,0.0,(1-i*0.1)*1.0e-15); - cout << i << " i,delta " << delta << endl << flush; // v0 = v0 + SpatialVector(0.0,0.0,0.9e-8); v0 = v0 - delta; v0.normalize(); // SpatialVector v0 ( 1, 0, 0); // Works for this one. STARE_ArrayIndexSpatialValue siv0 = index.ValueFromSpatialVector(v0,resolution); BitShiftNameEncoding rightJustified(sIndex.idByPoint(v0)); +#ifdef DIAG + cout << i << " i,delta " << delta << endl << flush; cout << "sI idbp 0x" << setw(16) << setfill('0') << hex << sIndex.idByPoint(v0) << dec << endl << flush; cout << "rj name "<< rightJustified.getName() << endl << flush; +#endif EmbeddedLevelNameEncoding leftJustified(rightJustified.leftJustifiedId()); EmbeddedLevelNameEncoding leftJustifiedWithResolution = leftJustified.atLevel(resolution, true); // True means keep all bits STARE_ArrayIndexSpatialValue siv = leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); - cout << "siv0 " << setw(16) << setfill('0') << hex << siv0 << endl << flush; - cout << "siv " << setw(16) << setfill('0') << hex << siv << endl << flush; SpatialVector v_siv0 = index.SpatialVectorFromValue(siv0); SpatialVector v_siv = index.SpatialVectorFromValue(siv); +#ifdef DIAG + cout << "siv0 " << setw(16) << setfill('0') << hex << siv0 << endl << flush; + cout << "siv " << setw(16) << setfill('0') << hex << siv << endl << flush; cout << setprecision(16); cout << "v0 " << v0 << endl << flush; cout << "v_siv0 " << v_siv0 << endl << flush; cout << "v_siv " << v_siv << endl << flush; cout << endl << "--" << endl << endl << flush; +#endif + ASSERT_EQUALDM("v0-to-siv0",v0,v_siv0,1.0e-8); + ASSERT_EQUALDM("v0-to-siv",v0,v_siv,1.0e-8); } - - } // FAIL(); diff --git a/tests/CUTE/SpatialIndex_test.cpp b/tests/CUTE/SpatialIndex_test.cpp index a127d6b..33d8c7a 100644 --- a/tests/CUTE/SpatialIndex_test.cpp +++ b/tests/CUTE/SpatialIndex_test.cpp @@ -24,10 +24,15 @@ void SpatialIndex_test() { // cout << 100 << endl << flush; - ASSERT_EQUAL("S010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,gPi)).nameByPoint(zhat)).c_str()); // SpatialVector zhat_r = -1.0*zhat; SpatialVector zhat_r = zhat.reverse(); + + ASSERT_EQUAL("N010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,0)).nameByPoint(zhat)).c_str()); + ASSERT_EQUAL("S010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,0)).nameByPoint(zhat_r)).c_str()); + + ASSERT_EQUAL("S010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,gPi)).nameByPoint(zhat)).c_str()); ASSERT_EQUAL("N010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,gPi)).nameByPoint(zhat_r)).c_str()); + // TODO Figure out how to have the following work. Probably by littering the call sequence with const. // ASSERT_EQUAL("N010000",string(SpatialIndex(maxlevel,buildlevel,SpatialRotation(xhat,gPi)).nameByPoint(zhat.reverse())).c_str()); From 5858c0a7431cc6917d9c611b79942a3d5469d697 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Fri, 20 Sep 2019 16:56:45 -0400 Subject: [PATCH 25/31] Abandoning for now use of spatial ids due to instability in idByPoint. --- contrib/python/example-viz-1.py | 125 +++++++++++++++++++++++--------- include/PySTARE.h | 4 +- include/STARE.h.in | 9 ++- src/PySTARE.cpp | 56 ++++++++++++++ src/PySTARE.i | 90 ++++++++++++++++++++++- src/STARE.C | 23 +++++- src/SpatialIndex.cpp | 40 ++++++++-- 7 files changed, 301 insertions(+), 46 deletions(-) diff --git a/contrib/python/example-viz-1.py b/contrib/python/example-viz-1.py index 845fce5..7914545 100644 --- a/contrib/python/example-viz-1.py +++ b/contrib/python/example-viz-1.py @@ -16,6 +16,7 @@ def shiftarg_lon(lon): return lon def triangulate(i0,i1,i2): + print('triangulating...') # i0,i1,i2,ic = ps.to_vertices(indices) i0lat,i0lon = ps.to_latlon(i0) i1lat,i1lon = ps.to_latlon(i1) @@ -35,6 +36,20 @@ def triangulate(i0,i1,i2): k=k+3 for i in range(len(lons)): lons[i] = shiftarg_lon(lons[i]) + print('triangulating done.') + return lons,lats,intmat + +def triangulate1(lats,lons): + print('triangulating1...') + intmat=[] + npts=int(len(lats)/3) + k=0 + for i in range(npts): + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating1 done.') return lons,lats,intmat # lat = np.array([0, 0,60], dtype=np.double) @@ -57,27 +72,57 @@ def test1(indices): triang = tri.Triangulation(lons,lats,intmat) return triang +# proj = ccrs.Mollweide() +proj = ccrs.PlateCarree() +# proj = ccrs.Robinson() +# proj = ccrs.Geodetic() +plt.figure() +plt.subplot(projection=proj) +# ax = plt.axes(projection=ccrs.PlateCarree()) +# ax = plt.axes(projection=ccrs.Mollweide()) +# proj=ccrs.Mollweide() +ax = plt.axes(projection=proj) +ax.set_global() +# ax.set_xlim(-180,180) +# ax.set_ylim(-90,90) +# ax.set_xlim(-1,1) +# ax.set_ylim(-1,1) +ax.coastlines() +# plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) +# plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) +# plt.triplot(triang,'ko-') +# plt.show() + def plot1(triang): - ax = plt.axes(projection=ccrs.PlateCarree()) - ax.set_xlim(-180,180) - ax.set_ylim(-90,90) - ax.coastlines() - # plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) - # plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) - # plt.triplot(triang,'ko-') - plt.triplot(triang,'r-') + # plt.triplot(triang,'ro-',transform=ccrs.Geodetic()) + plt.triplot(triang,'ro-',transform=proj) # plt.show() return plot1(test1(indices)) -level = 3 -hull = ps.to_hull_range(indices,level,100) -h0,h1,h2,hc = ps.to_vertices(hull) -lons1,lats1,intmat1 = triangulate(h0,h1,h2) +level = 4 +# hull = ps.to_hull_range(indices,level,100) +hull = ps.to_hull_range_from_latlon(lat,lon,level,100) +print('0 hull len: ',len(hull)) + +# print(90) +# lats0 = np.zeros(len(hull)*4,dtype=np.int64) +# lons0 = np.zeros(len(hull)*4,dtype=np.int64) +# lats0,lons0 = ps._to_vertices_latlon(hull) +# print('lats0,lons0: ',len(lats0),len(lons0)) +print(100) +lath,lonh,lathc,lonhc = ps.to_vertices_latlon(hull) +print('lath,lathc: ',len(lath),len(lathc)) +print(110) +lons1,lats1,intmat1 = triangulate1(lath,lonh) + +## h0,h1,h2,hc = ps.to_vertices(hull) +## lons1,lats1,intmat1 = triangulate(h0,h1,h2) + print('0 hull len: ',len(hull)) print('0 hull lats len: ',len(lats1)) -jtest=3 +jtest=18 j=jtest print('0 hull lats1: ',j,lats1[j*3:(j+1)*3]) print('0 hull lons1: ',j,lons1[j*3:(j+1)*3]) @@ -86,24 +131,24 @@ def plot1(triang): print('0 hull lons1: ',[i for i in lons1[j:j+12]]) print('') -tid = np.array([0x4c0000000000003],dtype=np.int64) -t0,t1,t2,tc = ps.to_vertices(tid) -print('t0: ',hex(t0[0])) -print('t1: ',hex(t1[0])) -print('t2: ',hex(t2[0])) -print('tc: ',hex(tc[0])) -print('t0 ll : ',ps.to_latlon(t0)) -print('t1 ll : ',ps.to_latlon(t1)) -print('t2 ll : ',ps.to_latlon(t2)) -print('tc ll : ',ps.to_latlon(tc)) +# tid = np.array([0x4c0000000000003],dtype=np.int64) +# t0,t1,t2,tc = ps.to_vertices(tid) +# print('t0: ',hex(t0[0])) +# print('t1: ',hex(t1[0])) +# print('t2: ',hex(t2[0])) +# print('tc: ',hex(tc[0])) +# print('t0 ll : ',ps.to_latlon(t0)) +# print('t1 ll : ',ps.to_latlon(t1)) +# print('t2 ll : ',ps.to_latlon(t2)) +# print('tc ll : ',ps.to_latlon(tc)) print('') -if True: +if False: # i=9; ilen=2 # i=10; ilen=1 # i=jtest; ilen=4 - i=jtest; ilen=1 + i=jtest; ilen=10 id_test = np.array(hull[i:i+ilen],dtype=np.int64) print('i,id : ',i,[hex(j) for j in id_test]) i0=i*3; i1=(i+ilen)*3 @@ -114,16 +159,23 @@ def plot1(triang): for j in range(ilen): intmat1.append([3*j,3*j+1,3*j+2]) # intmat1 = [[0,2,1]] - i0test,i1test,i2test,ictest = ps.to_vertices(id_test) - print('test id: ',[hex(i) for i in id_test]) - print('test ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) - print('test im : ',[i for i in intmat1]) - lonstest,latstest,intmattest = triangulate(i0test,i1test,i2test) - print('test lat: ',latstest) - print('test lon: ',lonstest) + id_test_lats,id_test_lons,id_test_latsc,id_test_lonsc = ps.to_vertices_latlon(id_test) + lonstest,latstest,intmattest = triangulate1(id_test_lats,id_test_lons) + # i0test,i1test,i2test,ictest = ps.to_vertices(id_test) + # print('test id: ',[hex(i) for i in id_test]) + # print('test ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) + # print('test im : ',[i for i in intmat1]) + # lonstest,latstest,intmattest = triangulate(i0test,i1test,i2test) + print('test lat: ',len(latstest)) + print('test lon: ',len(lonstest)) + print('test im: ',len(intmattest)) + print('test im[0]: ',intmattest[0]) + print('test im: ',intmattest) triangtest = tri.Triangulation(lonstest,latstest,intmattest) - plt.triplot(triangtest,'g-') - plt.scatter(lonstest,latstest,s=50,c='g') + plt.triplot(triangtest,'go-',transform=proj) + plt.scatter(lonstest,latstest,s=50,c='g',transform=proj) + # plt.triplot(triangtest,'go-',transform=ccrs.Geodetic()) + # plt.scatter(lonstest,latstest,s=50,c='g',transform=ccrs.Geodetic()) print('level : ',level) # print('hull : ',[hex(i) for i in hull]) @@ -134,7 +186,10 @@ def plot1(triang): print('hull im : ',[i for i in intmat1]) triang1 = tri.Triangulation(lons1,lats1,intmat1) -plt.triplot(triang1,'b-') +# plt.triplot(triang1,'bo-') +plt.triplot(triang1,'bo-',transform=proj) +# plt.triplot(triang1,'bo-',transform=ccrs.Geodetic()) +# plt.triplot(triang1,'bo-',transform=proj) # ax = plt.axes(projection=ccrs.PlateCarree()) # ax.set_xlim(-180,180) diff --git a/include/PySTARE.h b/include/PySTARE.h index 7129f78..96f51a8 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -35,11 +35,13 @@ void to_latlon(int64_t* indices, int len, double* lat, double* lon); void to_latlonlevel(int64_t* indices, int len, double* lat, double* lon, int* levels); void to_level(int64_t* indices, int len, int* levels); void to_triangle(int64_t* indices, int len); -void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid); +// Broken or dangerous. void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid); +void _to_vertices_latlon(int64_t* indices, int len, double* triangle_info_lats, int dmy1, double* triangle_info_lons, int dmy2 ); void to_area(int64_t* indices, int len, double* areas); void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri); void _to_hull_range(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); +void _to_hull_range_from_latlon(double* lat, int len_lat, double* lon, int len_lon, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni); diff --git a/include/STARE.h.in b/include/STARE.h.in index e66fd20..a39e252 100644 --- a/include/STARE.h.in +++ b/include/STARE.h.in @@ -43,8 +43,11 @@ typedef std::vector STARE_ArrayIndexSpatialValues; /// TODO The following is a stopgap on the way to completely removing this hardcoded bit. #define STARE_HARDWIRED_RESOLUTION_LEVEL_MAX 27 +// < 2019-0919 #define STARE_HARDWIRED_BUILD_LEVEL_MAX 5 #define STARE_HARDWIRED_BUILD_LEVEL_MAX 5 - +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 7 +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 8 +// #define STARE_HARDWIRED_BUILD_LEVEL_MAX 10 bad alloc class STARE { public: @@ -61,12 +64,16 @@ public: STARE_ArrayIndexSpatialValue ValueFromLatLonDegrees(float64 latDegrees, float64 lonDegrees, int resolutionLevel = STARE_HARDWIRED_RESOLUTION_LEVEL_MAX); LatLonDegrees64 LatLonDegreesFromValue(STARE_ArrayIndexSpatialValue spatialStareId); SpatialVector SpatialVectorFromValue(STARE_ArrayIndexSpatialValue spatialStareId); + STARE_ArrayIndexSpatialValue ValueFromSpatialVector(SpatialVector v, int resolution); uint32 ResolutionLevelFromValue(STARE_ArrayIndexSpatialValue spatialStareId); Triangle TriangleFromValue(STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel = -1); float64 AreaFromValue (STARE_ArrayIndexSpatialValue spatialStareId, int resolutionLevel = -1); + + STARE_ArrayIndexSpatialValues toVertices(STARE_ArrayIndexSpatialValues spatialStareIds); STARE_SpatialIntervals ConvexHull(LatLonDegrees64ValueVector points,int force_resolution_level); + STARE_SpatialIntervals ConvexHull(STARE_ArrayIndexSpatialValues points,int force_resolution_level); SpatialIndex getIndex() { return sIndex; } SpatialIndex getIndex(int resolutionLevel); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 633cf5d..7c58cbd 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -47,6 +47,9 @@ void to_triangle(int64_t* indices, int len) { } } +/* + * Broken or dangerous + * void to_vertices(int64_t* indices, int len, int64_t* vertices0, int64_t* vertices1, int64_t* vertices2, int64_t* centroid) { for (int i=0; i0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } + range_indices[i] = result[i]; + } + cout << dec << endl << flush; + result_size[0] = result.size(); +} + void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); // intersection[0] = 69; diff --git a/src/PySTARE.i b/src/PySTARE.i index 1b6c308..a7f24ad 100644 --- a/src/PySTARE.i +++ b/src/PySTARE.i @@ -133,6 +133,32 @@ $4 = (double*) array_data((PyArrayObject*)out2); } +/* maps ONE int64_t input vector to TWO double output vectors of 4 times the length */ +/* We use this to convert STARE index to lat+lon */ +%typemap(in, numinputs=1) + (int64_t* in_array, int length, double* out_array1, int dmy1, double* out_array2, int dmy2) + (PyObject* out1=NULL, PyObject* out2=NULL) +{ + int is_new_object=0; + npy_intp size[1] = { -1}; + PyArrayObject* array = obj_to_array_contiguous_allow_conversion($input, NPY_INT64, &is_new_object); + if (!array || !require_dimensions(array, 1)) SWIG_fail; + + size[0] = 4*PyArray_DIM(array, 0); + + out1 = PyArray_SimpleNew(1, size, NPY_DOUBLE); + if (!out1) SWIG_fail; + out2 = PyArray_SimpleNew(1, size, NPY_DOUBLE); + if (!out2) SWIG_fail; + + $1 = (int64_t*) array_data(array); + $2 = (int) array_size(array,0); + $3 = (double*) array_data((PyArrayObject*)out1); + $4 = (int) array_size(array,0); + $5 = (double*) array_data((PyArrayObject*)out2); + $6 = (int) array_size(array,0); +} + /* maps ONE int64_t input vector to TWO int64_t output vectors of the same length */ /* We use this to convert STARE intervals to start/terminator arrays to aid comparison */ %typemap(in, numinputs=1) @@ -198,7 +224,7 @@ PyArrayObject* array = obj_to_array_contiguous_allow_conversion($input, NPY_INT64, &is_new_object); if (!array || !require_dimensions(array, 1)) SWIG_fail; - size[0] = PyArray_DIM(array, 0); + size[0] = PyArray_DIM(array, 0); out1 = PyArray_SimpleNew(1, size, NPY_INT64); if (!out1) SWIG_fail; @@ -252,7 +278,13 @@ PyTuple_SetItem($result, 0, (PyObject*)out1$argnum); PyTuple_SetItem($result, 1, (PyObject*)out2$argnum); } - +%typemap(argout) + (int64_t* in_array, int length, double* out_array1, int dmy1, double* out_array2, int dmy2) +{ + $result = PyTuple_New(2); + PyTuple_SetItem($result, 0, (PyObject*)out1$argnum); + PyTuple_SetItem($result, 1, (PyObject*)out2$argnum); +} %typemap(argout) (int64_t* in_array, int length, int64_t* out_array1, int64_t* out_array2) { @@ -286,9 +318,17 @@ $result = (PyObject*)out$argnum; } +# %typemap(argout) +# (double* in_array1, int length1, double* in_array2, int length2, int resolution, int64_t* out_array1, int out_length1, int64_t* out_array2, int out_length2) +# { +# $result = (PyObject*)out$argnum; +# } + + /* Applying the typemaps */ %apply (double * IN_ARRAY1, int DIM1) { - (double* lat, int len_lat) + (double* lat, int len_lat), + (double* lon, int len_lon) } %apply (int64_t * IN_ARRAY1, int DIM1) { @@ -306,6 +346,11 @@ (int64_t* cmp, int len12) } +%apply (double * INPLACE_ARRAY1, int DIM1) { + (double* triangle_info_lats, int dmy1), + (double* triangle_info_lons, int dmy2) +} + # %apply (int64_t * ARGOUT_ARRAY1, int DIM1 ) { # (int64_t* out_array, int out_length) # (int64_t* range_indices, int len_ri) @@ -343,6 +388,10 @@ (int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ) } +%apply (int64_t* in_array, int length, double* out_array1, int dmy1, double* out_array2, int dmy2) { + (int64_t* indices, int len, double* triangle_info_lats, int dmy1, double* triangle_info_lons, int dmy2) +} + %pythonprepend from_utc(int64_t*, int, int64_t*, int) %{ import numpy datetime = datetime.astype(numpy.int64) @@ -371,6 +420,41 @@ def to_hull_range(indices,resolution,range_size_limit=1000): range_indices = range_indices[:result_size[0]] return range_indices +def to_hull_range_from_latlon(lat,lon,resolution,range_size_limit=1000): + out_length = range_size_limit + range_indices = numpy.full([out_length],-1,dtype=numpy.int64) + result_size = numpy.full([1],-1,dtype=numpy.int64) + _to_hull_range_from_latlon(lat,lon,resolution,range_indices,result_size) + range_indices = range_indices[:result_size[0]] + return range_indices + +def to_vertices_latlon(indices): + out_length = len(indices) + lats = numpy.zeros([4*out_length],dtype=numpy.double) + lons = numpy.zeros([4*out_length],dtype=numpy.double) + # _to_vertices_latlon(indices,lats,lons,0) + lats,lons = _to_vertices_latlon(indices) + latsv = numpy.zeros([3*out_length],dtype=numpy.double) + lonsv = numpy.zeros([3*out_length],dtype=numpy.double) + latc = numpy.zeros([out_length],dtype=numpy.double) + lonc = numpy.zeros([out_length],dtype=numpy.double) + + k=0; l=0 + for i in range(out_length): + latsv[l] = lats[ k ] + lonsv[l] = lons[ k ] + + latsv[l+1] = lats[ k+1 ] + lonsv[l+1] = lons[ k+1 ] + + latsv[l+2] = lats[ k+2 ] + lonsv[l+2] = lons[ k+2 ] + + latc [i] = lats[ k+3 ] + lonc [i] = lons[ k+3 ] + k = k + 4; l = l + 3 + return latsv,lonsv,latc,lonc + def cmp_spatial(indices1, indices2): out_length = len(indices1)*len(indices2) cmp = numpy.zeros([out_length],dtype=numpy.int64) diff --git a/src/STARE.C b/src/STARE.C index 842cd10..5949baf 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -131,7 +131,28 @@ STARE_ArrayIndexSpatialValue STARE::ValueFromLatLonDegrees( } STARE_ArrayIndexSpatialValue STARE::ValueFromSpatialVector(SpatialVector v, int resolution) { - BitShiftNameEncoding rightJustified(sIndex.idByPoint(v)); + uint64 htmID; + SpatialVector vtry(v); + int k = 3; + while(k>0) { + try { + --k; + htmID = sIndex.idByPoint(vtry); + } catch( SpatialException e ) { + cerr << e.what(); + if( k > 0 ) { + cerr << " " << k << " Trying again... " << endl << flush; + vtry = vtry + 1.0e-10*SpatialVector(1,0,0); vtry.normalize(); + } else { + cerr << endl << flush; + stringstream ss; ss << setprecision(16); + ss << "STARE::ValueFromSpatialVector can't find vector v= " << v << endl; + throw SpatialFailure(ss.str().c_str()); + } + + } + } + BitShiftNameEncoding rightJustified(htmID); EmbeddedLevelNameEncoding leftJustified(rightJustified.leftJustifiedId()); EmbeddedLevelNameEncoding leftJustifiedWithResolution = leftJustified.atLevel(resolution, true); // True means keep all bits return leftJustifiedWithResolution.getSciDBLeftJustifiedFormat(); diff --git a/src/SpatialIndex.cpp b/src/SpatialIndex.cpp index def40b0..c2ab035 100644 --- a/src/SpatialIndex.cpp +++ b/src/SpatialIndex.cpp @@ -15,6 +15,7 @@ //# Jul 25, 2002 : Gyorgy Fekete -- Added pointById() //# +#include #include #include @@ -1363,7 +1364,7 @@ SpatialIndex::idByPoint(SpatialVector & v) const { } } DIAGOUT2(cout,"1 preamble i: " << index << endl << flush;); - // int index_barycentric = index; + int index_barycentric = index; /**/ // if(index == 9) { @@ -1423,18 +1424,47 @@ SpatialIndex::idByPoint(SpatialVector & v) const { float64 dc_improvement = 2.0; int attempt = 0; int index_tried = index; + int index_last_found = index; int itry = 0; do { ++attempt; if(attempt == 9) { - throw SpatialFailure("SpatialIndex::idByPoint(sv): Lost Point Failure 1."); + stringstream ss; + float64 lat,lon; + v.getLatLonDegrees(lat, lon); + ss << setprecision(16); + ss << "SpatialIndex::idByPoint(sv): Lost Point Failure 1. No convergence." + << " point = " << v + << " index_barycentric = " << index_barycentric + << " index_last_found = " << index_last_found + << " latlon = " << lat << "," << lon + << " nudge " << ( nudge ? "true" : "false" ) + << " dc_improvement = " << dc_improvement + << " dc_start = " << dc_start + << " dc_end = " << dc_end + << endl << flush; + throw SpatialFailure(ss.str().c_str()); } if(attempt>1) { index = index_dcs_sort[itry]+1; if( index == index_tried ) { ++itry; if(itry == 8) { - throw SpatialFailure("SpatialIndex::idByPoint(sv): Lost Point Failure 2."); + stringstream ss; + float64 lat,lon; + v.getLatLonDegrees(lat, lon); + ss << setprecision(16); + ss << "SpatialIndex::idByPoint(sv): Lost Point Failure 2. No convergence." + << " point = " << v + << " index_barycentric = " << index_barycentric + << " index_last_found = " << index_last_found + << " latlon = " << lat << "," << lon + << " nudge " << ( nudge ? "true" : "false" ) + << " dc_improvement = " << dc_improvement + << " dc_start = " << dc_start + << " dc_end = " << dc_end + << endl << flush; + throw SpatialFailure(ss.str().c_str()); } index = index_dcs_sort[itry]+1; } @@ -1461,8 +1491,8 @@ SpatialIndex::idByPoint(SpatialVector & v) const { dc_end = dc.length(); dc_improvement = dc_end/dc_start; DIAG1(cout << "dc start, end, ratio: " << dc_start << " " << dc_end << " " << dc_improvement << endl << flush;); - } while ( dc_improvement > 0.25 && dc_end > 0.25 ); - + index_last_found = index; + } while ( dc_improvement > 0.125 && dc_end > 0.25 && !nudge ); // } /**/ From 194b54b0bde094ac8786657efb9d7787f7c359af Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Fri, 20 Sep 2019 20:05:56 -0400 Subject: [PATCH 26/31] Improved python example plotting. --- contrib/python/example-viz-1.py | 159 +++++++++++++++++--------------- src/PySTARE.cpp | 14 +-- 2 files changed, 92 insertions(+), 81 deletions(-) diff --git a/contrib/python/example-viz-1.py b/contrib/python/example-viz-1.py index 7914545..cebe7eb 100644 --- a/contrib/python/example-viz-1.py +++ b/contrib/python/example-viz-1.py @@ -10,78 +10,95 @@ import numpy as np def shiftarg_lon(lon): + "If lon is outside +/-180, then correct back." if(lon>180): return ((lon + 180.0) % 360.0)-180.0 else: return lon def triangulate(i0,i1,i2): - print('triangulating...') - # i0,i1,i2,ic = ps.to_vertices(indices) - i0lat,i0lon = ps.to_latlon(i0) - i1lat,i1lon = ps.to_latlon(i1) - i2lat,i2lon = ps.to_latlon(i2) - lats = np.zeros([3*len(i0lat)],dtype=np.double) - lons = np.zeros([3*len(i0lat)],dtype=np.double) - intmat = [] - k=0 - for i in range(len(i0)): - lats[k] = i0lat[i] - lons[k] = i0lon[i] - lats[k+1] = i1lat[i] - lons[k+1] = i1lon[i] - lats[k+2] = i2lat[i] - lons[k+2] = i2lon[i] - intmat.append([k,k+1,k+2]) - k=k+3 - for i in range(len(lons)): - lons[i] = shiftarg_lon(lons[i]) - print('triangulating done.') - return lons,lats,intmat + "Prepare data structures for tri.Triangulate." + print('triangulating...') + # i0,i1,i2,ic = ps.to_vertices(indices) + i0lat,i0lon = ps.to_latlon(i0) + i1lat,i1lon = ps.to_latlon(i1) + i2lat,i2lon = ps.to_latlon(i2) + lats = np.zeros([3*len(i0lat)],dtype=np.double) + lons = np.zeros([3*len(i0lat)],dtype=np.double) + intmat = [] + k=0 + for i in range(len(i0)): + lats[k] = i0lat[i] + lons[k] = i0lon[i] + lats[k+1] = i1lat[i] + lons[k+1] = i1lon[i] + lats[k+2] = i2lat[i] + lons[k+2] = i2lon[i] + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating done.') + return lons,lats,intmat def triangulate1(lats,lons): - print('triangulating1...') - intmat=[] - npts=int(len(lats)/3) - k=0 - for i in range(npts): - intmat.append([k,k+1,k+2]) - k=k+3 - for i in range(len(lons)): - lons[i] = shiftarg_lon(lons[i]) - print('triangulating1 done.') - return lons,lats,intmat - + "Prepare data for tri.Triangulate." + print('triangulating1...') + intmat=[] + npts=int(len(lats)/3) + k=0 + for i in range(npts): + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating1 done.') + return lons,lats,intmat + +# Make a triangle's lat & lon corners # lat = np.array([0, 0,60], dtype=np.double) # lon = np.array([0,60,30], dtype=np.double) - lat = np.array([0, 0,60], dtype=np.double) lon = np.array([0,60,30], dtype=np.double) -s_resolution = 6 +# Make the indices out of the corners +s_resolution = 27 indices = ps.from_latlon(lat, lon, s_resolution) +print('lat len: ',len(lat)) +print('indices len: ',len(indices)) +print('indices: ',[hex(i) for i in indices]) +# Prepare to triangulate the indices (as locations) def test1(indices): "Plot the outline of a triangle using first three indices." i0 = np.array([indices[0]], dtype=np.int64) i1 = np.array([indices[1]], dtype=np.int64) i2 = np.array([indices[2]], dtype=np.int64) lons,lats,intmat = triangulate(i0,i1,i2) - # print('lats: ',lats) - # print('lons: ',lons) - triang = tri.Triangulation(lons,lats,intmat) - return triang - -# proj = ccrs.Mollweide() -proj = ccrs.PlateCarree() + print('test1 intmat len: ',len(intmat)) + print('test1 lats: ',lats) + print('test1 lons: ',lons) + print('test1 intm: ',intmat) + return lons,lats,intmat + # triang = tri.Triangulation(lons,lats,intmat) + # print('test1 triang : ',triang) + # return triang + +# Set up the projection and transformation +# proj = ccrs.PlateCarree() # proj = ccrs.Robinson() # proj = ccrs.Geodetic() +# proj = ccrs.Geodesic() +proj = ccrs.Mollweide() +transf = ccrs.Geodetic() +# transf = ccrs.PlateCarree() + plt.figure() -plt.subplot(projection=proj) +plt.subplot(projection=proj,transform=transf) # ax = plt.axes(projection=ccrs.PlateCarree()) # ax = plt.axes(projection=ccrs.Mollweide()) # proj=ccrs.Mollweide() -ax = plt.axes(projection=proj) +ax = plt.axes(projection=proj,transform=transf) ax.set_global() # ax.set_xlim(-180,180) # ax.set_ylim(-90,90) @@ -95,15 +112,19 @@ def test1(indices): def plot1(triang): # plt.triplot(triang,'ro-',transform=ccrs.Geodetic()) - plt.triplot(triang,'ro-',transform=proj) + plt.triplot(triang,'r-',transform=transf) # plt.show() return -plot1(test1(indices)) +test1_lons,test1_lats,test1_intmat = test1(indices) +plot1(tri.Triangulation(test1_lons,test1_lats,test1_intmat)) +plt.scatter(test1_lons,test1_lats,s=5,c='r',transform=transf) -level = 4 -# hull = ps.to_hull_range(indices,level,100) -hull = ps.to_hull_range_from_latlon(lat,lon,level,100) +# resolution = 2; ntri=20 +resolution = 4; ntri=100 +# resolution = 7; ntri=1000 +# hull = ps.to_hull_range(indices,resolution,100) +hull = ps.to_hull_range_from_latlon(lat,lon,resolution,ntri) print('0 hull len: ',len(hull)) # print(90) @@ -111,10 +132,10 @@ def plot1(triang): # lons0 = np.zeros(len(hull)*4,dtype=np.int64) # lats0,lons0 = ps._to_vertices_latlon(hull) # print('lats0,lons0: ',len(lats0),len(lons0)) -print(100) +# print(100) lath,lonh,lathc,lonhc = ps.to_vertices_latlon(hull) print('lath,lathc: ',len(lath),len(lathc)) -print(110) +# print(110) lons1,lats1,intmat1 = triangulate1(lath,lonh) ## h0,h1,h2,hc = ps.to_vertices(hull) @@ -127,8 +148,8 @@ def plot1(triang): print('0 hull lats1: ',j,lats1[j*3:(j+1)*3]) print('0 hull lons1: ',j,lons1[j*3:(j+1)*3]) j=0 -print('0 hull lats1: ',[i for i in lats1[j:j+12]]) -print('0 hull lons1: ',[i for i in lons1[j:j+12]]) +# print('0 hull lats1: ',[i for i in lats1[j:j+12]]) +# print('0 hull lons1: ',[i for i in lons1[j:j+12]]) print('') # tid = np.array([0x4c0000000000003],dtype=np.int64) @@ -141,8 +162,7 @@ def plot1(triang): # print('t1 ll : ',ps.to_latlon(t1)) # print('t2 ll : ',ps.to_latlon(t2)) # print('tc ll : ',ps.to_latlon(tc)) - -print('') +# print('') if False: # i=9; ilen=2 @@ -172,33 +192,24 @@ def plot1(triang): print('test im[0]: ',intmattest[0]) print('test im: ',intmattest) triangtest = tri.Triangulation(lonstest,latstest,intmattest) - plt.triplot(triangtest,'go-',transform=proj) - plt.scatter(lonstest,latstest,s=50,c='g',transform=proj) + plt.triplot(triangtest,'g-',transform=transf) + plt.scatter(lonstest,latstest,s=5,c='g',transform=transf) # plt.triplot(triangtest,'go-',transform=ccrs.Geodetic()) - # plt.scatter(lonstest,latstest,s=50,c='g',transform=ccrs.Geodetic()) + # plt.scatter(lonstest,latstest,s=10,c='g',transform=ccrs.Geodetic()) -print('level : ',level) +print('resolution: ',resolution) # print('hull : ',[hex(i) for i in hull]) print('hull len: ',len(hull)) print('hull ll len: ',len(lons1)) print('hull im len: ',len(intmat1)) -print('hull ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) -print('hull im : ',[i for i in intmat1]) +# print('hull ll : ',[(lats1[i],lons1[i]) for i in range(len(lats1))]) +# print('hull im : ',[i for i in intmat1]) triang1 = tri.Triangulation(lons1,lats1,intmat1) -# plt.triplot(triang1,'bo-') -plt.triplot(triang1,'bo-',transform=proj) -# plt.triplot(triang1,'bo-',transform=ccrs.Geodetic()) -# plt.triplot(triang1,'bo-',transform=proj) - -# ax = plt.axes(projection=ccrs.PlateCarree()) -# ax.set_xlim(-180,180) -# ax.set_ylim(-90,90) -# ax.coastlines() +plt.triplot(triang1,'b-',transform=transf,lw=1,markersize=3) +plt.scatter(lons1,lats1,s=10,c='b',transform=transf) # # plt.contourf(xg,yg,v0g,60,transform=ccrs.PlateCarree()) -# # plt.scatter(xg_flat,yg_flat,s=300,c=v_flat) -# # plt.triplot(triang,'ko-') -# plt.triplot(triang1,'r-') + plt.show() diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 7c58cbd..1d5728d 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -161,13 +161,13 @@ void _to_hull_range(int64_t* indices, int len, int resolution, int64_t* range_in cout << "_to_hull_range-warning: range_indices.size = " << len_ri << " too small." << endl << flush; cout << "_to_hull_range-warning: result size = " << result.size() << "." << endl << flush; } - int k=10; - cout << "thr "; + // int k=10; + // cout << "thr "; for(int i=0; i < (len_ri < result.size() ? len_ri : result.size()); ++i) { - if(k-->0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } + // if(k-->0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } range_indices[i] = result[i]; } - cout << dec << endl << flush; + // cout << dec << endl << flush; result_size[0] = result.size(); } @@ -184,12 +184,12 @@ void _to_hull_range_from_latlon(double* lat, int len_lat, double* lon, int len_l cout << "_to_hull_range-warning: result size = " << result.size() << "." << endl << flush; } int k=10; - cout << "thr "; + // cout << "thr "; for(int i=0; i < (len_ri < result.size() ? len_ri : result.size()); ++i) { - if(k-->0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } + // if(k-->0) { cout << "0x" << setw(16) << setfill('0') << hex << result[i] << " "; } range_indices[i] = result[i]; } - cout << dec << endl << flush; + // cout << dec << endl << flush; result_size[0] = result.size(); } From 152c210e3814088e7e92e9fa6c3c5d9d60b7addc Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Fri, 20 Sep 2019 23:46:36 -0400 Subject: [PATCH 27/31] Intermediate commit. Debugging pystare interface. Esp. intersect. --- contrib/python/example-viz-2.py | 96 +++++++++++++++++++++++++++++++++ src/EmbeddedLevelNameEncoding.C | 2 + src/HtmRangeMultiLevel.cpp | 46 ++++++++++++---- src/NOTES | 11 ++++ src/PySTARE.cpp | 12 ++++- 5 files changed, 156 insertions(+), 11 deletions(-) create mode 100644 contrib/python/example-viz-2.py diff --git a/contrib/python/example-viz-2.py b/contrib/python/example-viz-2.py new file mode 100644 index 0000000..94b0675 --- /dev/null +++ b/contrib/python/example-viz-2.py @@ -0,0 +1,96 @@ + +from math import ceil +import pystare as ps + +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import cartopy.crs as ccrs + +import numpy as np + +def shiftarg_lon(lon): + "If lon is outside +/-180, then correct back." + if(lon>180): + return ((lon + 180.0) % 360.0)-180.0 + else: + return lon + +def triangulate1(lats,lons): + "Prepare data for tri.Triangulate." + print('triangulating1...') + intmat=[] + npts=int(len(lats)/3) + k=0 + for i in range(npts): + intmat.append([k,k+1,k+2]) + k=k+3 + for i in range(len(lons)): + lons[i] = shiftarg_lon(lons[i]) + print('triangulating1 done.') + return lons,lats,intmat + +def plot1(lon,lat,lons,lats,triang,c0='r',c1='b',transf=None,lw=1): + if(lon is not None): + x=np.zeros([lon.size+1],dtype=np.double);x[:-1]=lon[:];x[-1]=lon[0] + y=np.zeros([lat.size+1],dtype=np.double); y[:-1]=lat[:]; y[-1]=lat[0] + ax.plot(x,y,True,transform=transf,c=c0) + plt.triplot(triang,c1+'-',transform=transf,lw=lw,markersize=3) + plt.scatter(lons,lats,s=10,c=c1,transform=transf) + return + +def make_hull(lat0,lon0,resolution0,ntri0): + hull0 = ps.to_hull_range_from_latlon(lat0,lon0,resolution0,ntri0) + lath0,lonh0,lathc0,lonhc0 = ps.to_vertices_latlon(hull0) + lons0,lats0,intmat0 = triangulate1(lath0,lonh0) + triang0 = tri.Triangulation(lons0,lats0,intmat0) + return lats0,lons0,triang0,hull0 + +resolution = 5 + +resolution0 = resolution; ntri0 = 1000 +lat0 = np.array([ 10, 5, 60,70], dtype=np.double) +lon0 = np.array([-30,-20,60,10], dtype=np.double) +lats0,lons0,triang0,hull0 = make_hull(lat0,lon0,resolution0,ntri0) +print('hull0: ',len(hull0)) + +resolution1 = resolution; ntri1 = 1000 +lat1 = np.array([10, 20, 30, 20 ], dtype=np.double) +lon1 = np.array([-60, 60, 60, -60], dtype=np.double) +lats1,lons1,triang1,hull1 = make_hull(lat1,lon1,resolution1,ntri1) +print('hull1: ',len(hull1)) + +if True: + intersected = np.full([1000],-1,dtype=np.int64) + # intersected = ps.intersect(hull0,hull1,multiresolution=True) + print('hull0: ',[hex(i) for i in hull0]) + print('hull1: ',[hex(i) for i in hull1]) + ps._intersect_multiresolution(hull0,hull1,intersected) + # print('intersected: ',len(intersected)) + # print('np.min: ',np.amin(intersected)) + # print('intersected: ',[hex(i) for i in intersected]) + endarg = np.argmax(intersected < 0) + intersected = intersected[:endarg] + # intersected = ps.intersect(hull0,hull1) + print('intersected: ',len(intersected)) + lati,loni,latci,lonci = ps.to_vertices_latlon(intersected) + lonsi,latsi,intmati = triangulate1(lati,loni) + triangi = tri.Triangulation(lonsi,latsi,intmati) + +# Set up the projection and transformation +# proj = ccrs.PlateCarree() +# proj = ccrs.Robinson() +# proj = ccrs.Geodesic() +proj = ccrs.Mollweide() +transf = ccrs.Geodetic() +# transf = ccrs.PlateCarree() +plt.figure() +plt.subplot(projection=proj,transform=transf) +ax = plt.axes(projection=proj,transform=transf) +ax.set_global() +ax.coastlines() + +plot1(lon0,lat0,lons0,lats0,triang0,c0='r',c1='b',transf=transf) +plot1(lon1,lat1,lons1,lats1,triang1,c0='g',c1='c',transf=transf) +plot1(None,None,lonsi,latsi,triangi,c0='y',c1='r',transf=transf,lw=4) +plt.show() diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index 422c402..dd2063a 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -9,6 +9,7 @@ #include #include #include +#include EmbeddedLevelNameEncoding::EmbeddedLevelNameEncoding() {} @@ -444,6 +445,7 @@ uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int // Check for overflow. Not a completely accurate check. if( (successor & TopBit) == TopBit ) { + #ifdef DIAG cout << "n: " << setw(23) << dec << n << endl << flush; cout << "lowerBound: " << setw(23) << hex << lowerBound << endl << flush; diff --git a/src/HtmRangeMultiLevel.cpp b/src/HtmRangeMultiLevel.cpp index 49a6e4a..a67ddfe 100644 --- a/src/HtmRangeMultiLevel.cpp +++ b/src/HtmRangeMultiLevel.cpp @@ -1142,24 +1142,50 @@ void HtmRangeMultiLevel::CompressionPass() { // TODO Perhaps instead try a find or a search that would set the iterators. } else { // cout << "400: " << endl << flush; +// cout << "400: blo = 0x" << setw(16) << setfill('0') << hex << bareLo << endl << flush; +// cout << "400: bhi = 0x" << setw(16) << setfill('0') << hex << bareHi << endl << flush; +// cout << "400: dlt = 0x" << setw(16) << setfill('0') << hex << delta << endl << flush; +// cout << "400: dlt = " << setw(16) << dec << delta << endl << flush; // Snip off as much as possible int numberToCoalesce = (delta+1) / 4; + // cout << "410: ntc " << dec << numberToCoalesce << endl << flush; Key oldLo = lo0; Key newLo = oldLo; - for(int i=0; iincrement(newLo,level0); - } - } + // cout << "420: lo0 = 0x" << setw(16) << setfill('0') << hex << lo0 << endl << flush; + +// for(int i=0; iincrement(newLo,level0); +// } +// } + my_los->free(oldLo); --oldLo; // Reduce level // my_los->insert(oldLo,1025); my_los->insert(oldLo,my_his->getkey()); // What's the current key (for my_his)? // TODO Don't worry. - Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); - if(newLoPredecessor != hi0) { - my_his->insert(newLoPredecessor,1025); - // my_los->insert(newLo,1025); - my_los->insert(newLo,newLoPredecessor); + + try { + // Scoop up a bunch + for(int i=0; iincrement(newLo,level0); + } + } + // Then break it in half. + Key newLoPredecessor = encoding->predecessorToLowerBound_NoDepthBit(newLo,level0); + if(newLoPredecessor != hi0) { + my_his->insert(newLoPredecessor,1025); + // my_los->insert(newLo,1025); + my_los->insert(newLo,newLoPredecessor); + } + } catch ( SpatialException &e ) { + // What if we're at the top index already? Ooops, there's no new low at the new break. + // cout << "400: " << e.what() << endl << flush; + if( string(e.what()) != string("EmbeddedLevelNameEncoding::error-increment-overflow") ) { + throw SpatialFailure("HtmRangeMultiLevel::Compress::unknown error while incrementing to newLo."); + } } + + my_los->reset(); my_his->reset(); // TODO Reset is too drastic. Prefer to step back a little... Bad, bad, bad. } } diff --git a/src/NOTES b/src/NOTES index 23b228e..23d2238 100644 --- a/src/NOTES +++ b/src/NOTES @@ -1,3 +1,14 @@ + + + +2019-0920 HtmRangeMultiLevel Compress bug + +TODO add a unit test for edge cases for this file +The compress bug resulted from the logic not accounting for what happens when you're at the last index value at a level. I.e. you're at the end of the road. + + +2019 Early + test if a triangle (id) is rejected: Implicitly, id is used in macro NV diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 1d5728d..07a608f 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -207,16 +207,26 @@ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_ } void _intersect_multiresolution(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_t* intersection, int leni) { + cout << "lens: " << len1 << " " << len2 << " " << leni << endl << flush; STARE_SpatialIntervals si1(indices1, indices1+len1), si2(indices2, indices2+len2); + cout << 100 << endl << flush; + // intersection[0] = 69; SpatialRange r1(si1), r2(si2); + cout << 200 << endl << flush; // SpatialRange *ri = r1 & r2; SpatialRange *ri = sr_intersect(r1,r2,true); - STARE_SpatialIntervals result = ri->toSpatialIntervals(); + // SpatialRange *ri = sr_intersect(r1,r2,false); + cout << 300 << endl << flush; + STARE_SpatialIntervals result_intervals = ri->toSpatialIntervals(); + STARE_ArrayIndexSpatialValues result = expandIntervals(result_intervals); + cout << 400 << endl << flush; leni = result.size(); + cout << 500 << " result size " << result.size() << endl << flush; for(int i=0; i Date: Sat, 21 Sep 2019 16:08:24 -0400 Subject: [PATCH 28/31] Found EmbeddedLevelNameEnc. was returning incorrect one_at_level. ExpandInterval attempted to avoid a function call and retrieved a mask and an interval increment (i.e. one_at_level), but returned it's own internal "one" and not the SciDB "one", which is shifted over one position. Compare increment_LevelToMaskDelta and SciDBincrement_LevelToMaskDelta for details. Note STARE uses SciDB formatted values. This error was causing SpatialRange's intersect to appear to drop triangles. --- include/EmbeddedLevelNameEncoding.h | 5 +++ src/EmbeddedLevelNameEncoding.C | 8 ++++ src/STARE.C | 41 +++++++++++++---- tests/CUTE/EmbeddedLevelNameEncoding_test.cpp | 42 ++++++++++++++++++ tests/CUTE/STARE_test.cpp | 44 +++++++++++++++---- tests/CUTE/SpatialRange_test.cpp | 4 +- tests/CUTE/Test.cpp | 5 ++- 7 files changed, 128 insertions(+), 21 deletions(-) diff --git a/include/EmbeddedLevelNameEncoding.h b/include/EmbeddedLevelNameEncoding.h index 4853b1b..15c9e1c 100644 --- a/include/EmbeddedLevelNameEncoding.h +++ b/include/EmbeddedLevelNameEncoding.h @@ -102,7 +102,12 @@ class EmbeddedLevelNameEncoding: public NameEncoding { EmbeddedLevelNameEncoding clearDeeperThanLevel(uint64 level); + void SciDBincrement_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const; void increment_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const; + + /* + * increment and decrement work on the internal bit format, not the SciDB format. + */ uint64 increment(uint64 lowerBound, uint32 level, int steps = 1) const; uint64 decrement(uint64 lowerBound, uint32 level, int steps = 1) const; diff --git a/src/EmbeddedLevelNameEncoding.C b/src/EmbeddedLevelNameEncoding.C index dd2063a..0fdf1f4 100644 --- a/src/EmbeddedLevelNameEncoding.C +++ b/src/EmbeddedLevelNameEncoding.C @@ -391,6 +391,14 @@ void EmbeddedLevelNameEncoding::increment_LevelToMaskDelta(uint32 level,uint64 & // one_at_level = one_at_level >> 2; } +void EmbeddedLevelNameEncoding::SciDBincrement_LevelToMaskDelta(uint32 level,uint64 &one_mask_to_level,uint64 &one_at_level) const { + increment_LevelToMaskDelta(level,one_mask_to_level,one_at_level); + one_mask_to_level = one_mask_to_level >> 1; + one_at_level = one_at_level >> 1; + +// one_at_level = one_at_level >> 2; +} + uint64 EmbeddedLevelNameEncoding::increment(uint64 lowerBound, uint32 level, int n) const { /// TODO Error checking of overflow not trustworthy here. using namespace std; diff --git a/src/STARE.C b/src/STARE.C index 5949baf..2310424 100644 --- a/src/STARE.C +++ b/src/STARE.C @@ -16,6 +16,12 @@ #include #include +#ifndef DIAG +#define DIAGOUT1(a) +#else +#define DIAGOUT1(a) a +#endif + /** * @brief Version function with C linkage to aid in finding the library with autoconf * @return The library's version. Points to static storage. @@ -721,17 +727,18 @@ uint64 spatialLevelMask() { STARE_ArrayIndexSpatialValues expandInterval(STARE_SpatialIntervals interval, int64 force_resolution) { // STARE_SpatialIntervals interval should just be one interval, i.e. a value or value+terminator. - // cout << dec << 200 << endl << flush; + DIAGOUT1(cout << endl << dec << 200 << endl << flush;) STARE_ArrayIndexSpatialValue siv0 = interval[0]; EmbeddedLevelNameEncoding leftJustified; - // cout << dec << 220 << endl << flush; + DIAGOUT1(cout << dec << 220 << endl << flush;) uint64 return_resolution = siv0 & leftJustified.levelMaskSciDB; - // cout << dec << 225 << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush; + DIAGOUT1(cout << dec << 225 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) if( force_resolution > -1 ) { siv0 = ( siv0 & ~leftJustified.levelMaskSciDB ) | force_resolution; return_resolution = force_resolution; } - // cout << dec << 230 << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush; + DIAGOUT1(cout << dec << 229 << " f & resolution: " << dec << force_resolution << " " << return_resolution << endl << flush;) + DIAGOUT1(cout << dec << 230 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) leftJustified.setIdFromSciDBLeftJustifiedFormat(siv0); // cout << dec << 235 << endl << flush; STARE_ArrayIndexSpatialValue siv_term; @@ -740,23 +747,39 @@ STARE_ArrayIndexSpatialValues expandInterval(STARE_SpatialIntervals interval, in } else { siv_term = leftJustified.getSciDBTerminatorLeftJustifiedFormat(); } - // cout << dec << 240 << endl << flush; + DIAGOUT1(cout << dec << 239 << " " << setw(16) << setfill('0') << hex << siv_term << dec << endl << flush;) + DIAGOUT1(cout << endl << dec << 240 << endl << flush;) uint64 one_mask_to_resolution, one_at_resolution; - leftJustified.increment_LevelToMaskDelta(siv0 & leftJustified.levelMaskSciDB,one_mask_to_resolution,one_at_resolution); - // cout << dec << 245 << endl << flush; + leftJustified.SciDBincrement_LevelToMaskDelta(siv0 & leftJustified.levelMaskSciDB,one_mask_to_resolution,one_at_resolution); + // cout << dec << 242 << endl << flush; + + uint64 delta = ((siv_term+1)-(siv0 & ~leftJustified.levelMaskSciDB)); + + DIAGOUT1(cout << endl;) + uint64 one = 1; + DIAGOUT1(cout << dec << 243 << " " << setw(16) << setfill('0') << hex << (one << (63-3-2*return_resolution)) << dec << endl << flush;) + DIAGOUT1(cout << dec << 244 << " " << setw(16) << setfill('0') << hex << leftJustified.getSciDBTerminatorLeftJustifiedFormat() << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << siv_term << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) + DIAGOUT1(cout << dec << 245 << " " << setw(16) << setfill('0') << hex << delta << " " << dec << (delta/one_at_resolution) << endl << flush;) + // Give as much rope as needed. + DIAGOUT1(cout << dec << 246 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) siv0 = (siv0 & ~one_mask_to_resolution) | return_resolution; + DIAGOUT1(cout << dec << 247 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << endl << flush;) + DIAGOUT1(cout << dec << 247 << " " << setw(16) << setfill('0') << hex << one_at_resolution << dec << endl << endl << flush;) STARE_ArrayIndexSpatialValues expanded_interval; while( siv0 < siv_term ) { + DIAGOUT1(cout << dec << 249 << " " << setw(16) << setfill('0') << hex << siv0 << dec << endl << flush;) expanded_interval.push_back(siv0); siv0 += one_at_resolution; } - // cout << dec << 250 << endl << flush; + DIAGOUT1(cout << dec << 250 << endl << endl << flush;) return expanded_interval; } /** - * + * TODO Fix expandIntervals... */ STARE_ArrayIndexSpatialValues expandIntervals(STARE_SpatialIntervals intervals, int64 force_resolution) { STARE_ArrayIndexSpatialValues expanded_values; diff --git a/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp index c4d7d22..f4f2e87 100644 --- a/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp +++ b/tests/CUTE/EmbeddedLevelNameEncoding_test.cpp @@ -55,4 +55,46 @@ void EmbeddedLevelNameEncoding_test() { ASSERT(leftJustified.terminatorp(term_sid1)); ASSERT_EQUALM("getSciDBTerminatorLeftJustifiedFormat idempotent?",term_sid1,term_sid2); + + if(true) { + STARE_ArrayIndexSpatialValue siv0 = 0x2324000000000005; + EmbeddedLevelNameEncoding lj; lj.setIdFromSciDBLeftJustifiedFormat(siv0); + STARE_ArrayIndexSpatialValue siv0_term = lj.getSciDBTerminatorLeftJustifiedFormat(); + + STARE_ArrayIndexSpatialValue not_siv1 = lj.increment(lj.getId(), lj.getLevel(), 1); + EmbeddedLevelNameEncoding lj1(not_siv1); + STARE_ArrayIndexSpatialValue siv1 = lj1.getSciDBLeftJustifiedFormat(); + STARE_ArrayIndexSpatialValue sivt_term = 0x2327ffffffffffff; + + int level = lj.getLevel(); + uint64 one_mask_to_level, one_at_level; + lj.SciDBincrement_LevelToMaskDelta(level, one_mask_to_level, one_at_level); + STARE_ArrayIndexSpatialValue siv2 = siv0 + one_at_level; + STARE_ArrayIndexSpatialValue siv2_term = siv2 | one_mask_to_level; + + STARE_ArrayIndexSpatialValue siv0_tchk = (siv0 | one_mask_to_level); + +#define IDOUT(p,m,s) +// #define IDOUT(p,m,s) p << m << " " << setw(16) << setfill(' ') << hex << s << dec << endl << flush; + IDOUT(cout,"siv0 ",siv0); + IDOUT(cout,"siv0_term ",siv0_term); + IDOUT(cout,"siv0_tchk ",siv0_tchk); + IDOUT(cout,"not_siv1 ",not_siv1); + IDOUT(cout,"siv1 ",siv1); + IDOUT(cout,"sivt_term ",sivt_term); + IDOUT(cout,"siv2 ",siv2); + IDOUT(cout,"siv2_term ",siv2_term); + // IDOUT(cout,"lj id ",lj.getId()); +#undef IDOUT + ASSERT_EQUAL( siv0 , 0x2324000000000005 ); + ASSERT_EQUAL( siv0_term , 0x2325ffffffffffff ); + ASSERT_EQUAL( siv0_tchk , 0x2325ffffffffffff ); + ASSERT_EQUAL( not_siv1 , 0xc64c000000000005 ); + ASSERT_EQUAL( siv1 , 0x2326000000000005 ); + ASSERT_EQUAL( sivt_term , 0x2327ffffffffffff ); + ASSERT_EQUAL( siv2 , 0x2326000000000005 ); + ASSERT_EQUAL( siv2_term , 0x2327ffffffffffff ); + } + + } diff --git a/tests/CUTE/STARE_test.cpp b/tests/CUTE/STARE_test.cpp index 9f8437f..3b93c2c 100644 --- a/tests/CUTE/STARE_test.cpp +++ b/tests/CUTE/STARE_test.cpp @@ -1019,13 +1019,21 @@ void STARE_test() { } } +// #define DIAG +#ifndef DIAG +#define DIAGOUT2(p,m) #define SIVOUT(m,siv) -// #define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SIVSOUT(p,m,v) +#else +#define DIAGOUT2(p,m) p << m; +#define SIVOUT(m,siv) cout << m << " " << setw(16) << setfill('0') << hex << siv << dec << endl << flush; +#define SIVSOUT(p,m,v) { p << m << " "; for(int l=0; l Date: Sat, 21 Sep 2019 16:22:41 -0400 Subject: [PATCH 29/31] Added expandIntervals to _intersect. --- contrib/python/example-viz-2.py | 12 +++++++----- src/PySTARE.cpp | 15 ++++++++------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/contrib/python/example-viz-2.py b/contrib/python/example-viz-2.py index 94b0675..53e3264 100644 --- a/contrib/python/example-viz-2.py +++ b/contrib/python/example-viz-2.py @@ -62,15 +62,17 @@ def make_hull(lat0,lon0,resolution0,ntri0): if True: intersected = np.full([1000],-1,dtype=np.int64) + intersected = ps.intersect(hull0,hull1,multiresolution=False) # intersected = ps.intersect(hull0,hull1,multiresolution=True) - print('hull0: ',[hex(i) for i in hull0]) - print('hull1: ',[hex(i) for i in hull1]) - ps._intersect_multiresolution(hull0,hull1,intersected) + # print('hull0: ',[hex(i) for i in hull0]) + # print('hull1: ',[hex(i) for i in hull1]) + # ps._intersect_multiresolution(hull0,hull1,intersected) # print('intersected: ',len(intersected)) # print('np.min: ',np.amin(intersected)) # print('intersected: ',[hex(i) for i in intersected]) - endarg = np.argmax(intersected < 0) - intersected = intersected[:endarg] + # The following are for _intersect_multiresolution's results + # endarg = np.argmax(intersected < 0) + # intersected = intersected[:endarg] # intersected = ps.intersect(hull0,hull1) print('intersected: ',len(intersected)) lati,loni,latci,lonci = ps.to_vertices_latlon(intersected) diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 07a608f..3f63f07 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -199,7 +199,8 @@ void _intersect(int64_t* indices1, int len1, int64_t* indices2, int len2, int64_ SpatialRange r1(si1), r2(si2); // SpatialRange *ri = r1 & r2; SpatialRange *ri = sr_intersect(r1,r2,false); - STARE_SpatialIntervals result = ri->toSpatialIntervals(); + STARE_SpatialIntervals result_intervals = ri->toSpatialIntervals(); + STARE_ArrayIndexSpatialValues result = expandIntervals(result_intervals); leni = result.size(); for(int i=0; itoSpatialIntervals(); STARE_ArrayIndexSpatialValues result = expandIntervals(result_intervals); - cout << 400 << endl << flush; + // cout << 400 << endl << flush; leni = result.size(); - cout << 500 << " result size " << result.size() << endl << flush; + // cout << 500 << " result size " << result.size() << endl << flush; for(int i=0; i Date: Sat, 21 Sep 2019 19:36:30 -0400 Subject: [PATCH 30/31] Added pystare.expand_intervals. --- contrib/python/test.py | 188 +++++++++++++++++++++++++++++++++-------- include/PySTARE.h | 4 +- src/PySTARE.cpp | 14 +++ src/PySTARE.i | 53 ++++++++++++ 4 files changed, 223 insertions(+), 36 deletions(-) diff --git a/contrib/python/test.py b/contrib/python/test.py index a95ea63..10b884a 100644 --- a/contrib/python/test.py +++ b/contrib/python/test.py @@ -48,30 +48,6 @@ print('5 intersected: ',[hex(i) for i in intersected]) print('') -vertices0 = numpy.zeros([3], dtype=numpy.int64) -vertices1 = numpy.zeros([3], dtype=numpy.int64) -vertices2 = numpy.zeros([3], dtype=numpy.int64) -centroid = numpy.zeros([3], dtype=numpy.int64) - -vertices0,vertices1,vertices2,centroid = pystare.to_vertices(indices) -print('vertices0: ',[hex(i) for i in vertices0]) -print('vertices1: ',[hex(i) for i in vertices1]) -print('vertices2: ',[hex(i) for i in vertices2]) -print('centroid: ',[hex(i) for i in centroid]) -print('') - -vertices0_lats,vertices0_lons = pystare.to_latlon(vertices0) -vertices1_lats,vertices1_lons = pystare.to_latlon(vertices1) -vertices2_lats,vertices2_lons = pystare.to_latlon(vertices2) -centroid_lats,centroid_lons = pystare.to_latlon(centroid) - -for i in range(len(vertices0_lats)): - print(i," vert0 lat,lon: ",vertices0_lats[i],vertices0_lons[i]) - print(i," vert1 lat,lon: ",vertices1_lats[i],vertices1_lons[i]) - print(i," vert2 lat,lon: ",vertices2_lats[i],vertices2_lons[i]) - print(i," centr lat,lon: ",centroid_lats[i],centroid_lons[i]) - print("") - indices1 = numpy.array([0,0,0], dtype=numpy.int64) pystare._to_compressed_range(indices,indices1) print('_indices1: ',[hex(i) for i in indices1]) @@ -129,14 +105,156 @@ print('cmp1: ',cmp1) print('') -tid = numpy.array([0x4c0000000000003],dtype=numpy.int64) -t0,t1,t2,tc = pystare.to_vertices(tid) -print('t0: ',hex(t0[0])) -print('t1: ',hex(t1[0])) -print('t2: ',hex(t2[0])) -print('tc: ',hex(tc[0])) -print('t0 ll : ',pystare.to_latlon(t0)) -print('t1 ll : ',pystare.to_latlon(t1)) -print('t2 ll : ',pystare.to_latlon(t2)) -print('tc ll : ',pystare.to_latlon(tc)) -print('') + +print('Expand intervals') +intervals_src = [\ + 0x2320000000000005,\ + 0x2324000000000005,\ + 0x2327ffffffffffff,\ + 0x3aa0000000000005,\ + 0x3aa7ffffffffffff,\ + 0x3aa8000000000004,\ + 0x3ab2000000000005,\ + 0x3ac7ffffffffffff,\ + 0x3ad0000000000005,\ + 0x3ae7ffffffffffff,\ + 0x3af2000000000005,\ + 0x3afa000000000005,\ + 0x3b40000000000004,\ + 0x3b4c000000000005,\ + 0x3b52000000000005,\ + 0x3b5a000000000005,\ + 0x3b5fffffffffffff,\ + 0x3b80000000000004,\ + 0x3b8a000000000005,\ + 0x3b8fffffffffffff,\ + 0x3b90000000000004,\ + 0x3b9fffffffffffff,\ + 0x3bc0000000000005,\ + 0x3bc7ffffffffffff,\ + 0x3bc8000000000004,\ + 0x3bd0000000000005,\ + 0x3bdfffffffffffff,\ + 0x3be2000000000005,\ + 0x3be8000000000004,\ + 0x3bf4000000000005,\ + 0x3bf8000000000005,\ + 0x3bfc000000000005,\ + 0x3bffffffffffffff,\ + 0x3e40000000000005,\ + 0x3e43ffffffffffff,\ + 0x3e46000000000005,\ + 0x3f20000000000005,\ + 0x3f24000000000005,\ + 0x3f27ffffffffffff,\ + 0x3f32000000000005,\ + 0x3f36000000000005,\ + 0x3f3a000000000005,\ + 0x3fa0000000000005\ +] + +intervals_expanded_src = [\ + 0x2320000000000005 ,\ + 0x2324000000000005 ,\ + 0x2326000000000005 ,\ + 0x3aa0000000000005 ,\ + 0x3aa2000000000005 ,\ + 0x3aa4000000000005 ,\ + 0x3aa6000000000005 ,\ + 0x3aa8000000000004 ,\ + 0x3ab2000000000005 ,\ + 0x3ab4000000000005 ,\ + 0x3ab6000000000005 ,\ + 0x3ab8000000000005 ,\ + 0x3aba000000000005 ,\ + 0x3abc000000000005 ,\ + 0x3abe000000000005 ,\ + 0x3ac0000000000005 ,\ + 0x3ac2000000000005 ,\ + 0x3ac4000000000005 ,\ + 0x3ac6000000000005 ,\ + 0x3ad0000000000005 ,\ + 0x3ad2000000000005 ,\ + 0x3ad4000000000005 ,\ + 0x3ad6000000000005 ,\ + 0x3ad8000000000005 ,\ + 0x3ada000000000005 ,\ + 0x3adc000000000005 ,\ + 0x3ade000000000005 ,\ + 0x3ae0000000000005 ,\ + 0x3ae2000000000005 ,\ + 0x3ae4000000000005 ,\ + 0x3ae6000000000005 ,\ + 0x3af2000000000005 ,\ + 0x3afa000000000005 ,\ + 0x3b40000000000004 ,\ + 0x3b4c000000000005 ,\ + 0x3b52000000000005 ,\ + 0x3b5a000000000005 ,\ + 0x3b5c000000000005 ,\ + 0x3b5e000000000005 ,\ + 0x3b80000000000004 ,\ + 0x3b8a000000000005 ,\ + 0x3b8c000000000005 ,\ + 0x3b8e000000000005 ,\ + 0x3b90000000000004 ,\ + 0x3b98000000000004 ,\ + 0x3bc0000000000005 ,\ + 0x3bc2000000000005 ,\ + 0x3bc4000000000005 ,\ + 0x3bc6000000000005 ,\ + 0x3bc8000000000004 ,\ + 0x3bd0000000000005 ,\ + 0x3bd2000000000005 ,\ + 0x3bd4000000000005 ,\ + 0x3bd6000000000005 ,\ + 0x3bd8000000000005 ,\ + 0x3bda000000000005 ,\ + 0x3bdc000000000005 ,\ + 0x3bde000000000005 ,\ + 0x3be2000000000005 ,\ + 0x3be8000000000004 ,\ + 0x3bf4000000000005 ,\ + 0x3bf8000000000005 ,\ + 0x3bfc000000000005 ,\ + 0x3bfe000000000005 ,\ + 0x3e40000000000005 ,\ + 0x3e42000000000005 ,\ + 0x3e46000000000005 ,\ + 0x3f20000000000005 ,\ + 0x3f24000000000005 ,\ + 0x3f26000000000005 ,\ + 0x3f32000000000005 ,\ + 0x3f36000000000005 ,\ + 0x3f3a000000000005 ,\ + 0x3fa0000000000005 \ +] +intervals = numpy.array(intervals_src,dtype=numpy.int64) +expected_expanded = numpy.array(intervals_expanded_src,dtype=numpy.int64) +intervals_expanded = numpy.zeros([len(expected_expanded)],dtype=numpy.int64) +expected_len = numpy.zeros([1],dtype=numpy.int64) + +print('len(intervals) ',len(intervals)) +print('len(expected_expanded) ',len(expected_expanded)) +print('len(intervals_expanded) ',len(intervals_expanded)) +print('len(expected_len) ',len(expected_len)) +intervals_len = len(intervals) +resolution = -1 +pystare._expand_intervals(intervals,resolution,intervals_expanded,expected_len) +# print('intervals_expanded: ',[hex(i) for i in intervals_expanded]) +print('expected_len: ',expected_len) +error_found = False +for i in range(len(intervals_expanded)): + if(intervals_expanded[i] != expected_expanded[i]): + error_found = True + print('_expanded_intervals error at i = ',i,' ',intervals_expanded[i],' != ',expected_expanded[i]) +if(not error_found): + print('_expanded_intervals ok') + +intervals_expanded_1 = pystare.expand_intervals(intervals,resolution) +for i in range(len(intervals_expanded_1)): + if(intervals_expanded_1[i] != expected_expanded[i]): + error_found = True + print('expanded_intervals error at i = ',i,' ',intervals_expanded_1[i],' != ',expected_expanded[i]) +if(not error_found): + print('expanded_intervals ok') diff --git a/include/PySTARE.h b/include/PySTARE.h index 96f51a8..a5876c4 100644 --- a/include/PySTARE.h +++ b/include/PySTARE.h @@ -40,7 +40,9 @@ void _to_vertices_latlon(int64_t* indices, int len, double* triangle_info_lats, void to_area(int64_t* indices, int len, double* areas); void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri); -void _to_hull_range(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); +void _to_hull_range (int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); +void _expand_intervals(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); + void _to_hull_range_from_latlon(double* lat, int len_lat, double* lon, int len_lon, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs); void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_t* indices_terminators ); diff --git a/src/PySTARE.cpp b/src/PySTARE.cpp index 3f63f07..bc17f7d 100644 --- a/src/PySTARE.cpp +++ b/src/PySTARE.cpp @@ -144,6 +144,20 @@ void from_intervals(int64_t* intervals, int len, int64_t* indices_starts, int64_ // } } +void _expand_intervals(int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs) { + STARE_SpatialIntervals si(indices, indices+len); + STARE_ArrayIndexSpatialValues result = expandIntervals(si,resolution); + if(len_ri < result.size()) { + cout << dec; + cout << "_expand_intervals-warning: range_indices.size = " << len_ri << " too small." << endl << flush; + cout << "_expand_intervals-warning: result size = " << result.size() << "." << endl << flush; + } + for(int i=0; i < (len_ri < result.size() ? len_ri : result.size()); ++i) { + range_indices[i] = result[i]; + } + result_size[0] = result.size(); +} + void _to_compressed_range(int64_t* indices, int len, int64_t* range_indices, int len_ri) { STARE_SpatialIntervals si(indices, indices+len); SpatialRange r(si); diff --git a/src/PySTARE.i b/src/PySTARE.i index a7f24ad..1ba6aae 100644 --- a/src/PySTARE.i +++ b/src/PySTARE.i @@ -243,6 +243,35 @@ $6 = (int64_t*) array_data((PyArrayObject*)out4); } +/* We use this to create a STARE array ... _to_expand_intervals +%typemap(in, numinputs=1) + (int64_t* in_array, int length, int resolution, int64_t* out_array1, int dmy1, int64_t* out_array2, int dmy2) + (PyObject* out1=NULL, PyObject* out2=NULL) +{ + int is_new_object=0; + npy_intp size[1] = { -1}; + PyArrayObject* array = obj_to_array_contiguous_allow_conversion($input, NPY_INT64, &is_new_object); + if (!array || !require_dimensions(array, 1)) SWIG_fail; + + size[0] = PyArray_DIM(array, 0); + + out1 = PyArray_SimpleNew(1, size, NPY_INT64); + if (!out1) SWIG_fail; + + size[0] = 1; + out2 = PyArray_SimpleNew(1, size, NPY_INT64); + if (!out2) SWIG_fail; + + $1 = (int64_t*) array_data(array); + $2 = (int) array_size(array,0); + $3 = (int) array_size(array,0); + $4 = (int64_t*) array_data((PyArrayObject*)out1); + $5 = (int) array_size(array,0); + $6 = (int64_t*) array_data((PyArrayObject*)out2); + $7 = (int) array_size(array,0); +} +*/ + /****************/ /* OUT typemaps */ /****************/ @@ -312,6 +341,16 @@ PyTuple_SetItem($result, 3, (PyObject*)out4$argnum); } +/* expand_intervals +%typemap(argout) + (int64_t* in_array, int length, int resolution, int64_t* out_array1, int dmy1, int64_t* out_array2, int dmy2) +{ + $result = PyTuple_New(2); + PyTuple_SetItem($result, 0, (PyObject*)out1$argnum); + PyTuple_SetItem($result, 1, (PyObject*)out2$argnum); +} + */ + %typemap(argout) (int64_t* in_array1, int length1, int64_t* in_array2, int length2, int64_t* out_array, int out_length) { @@ -326,6 +365,7 @@ /* Applying the typemaps */ + %apply (double * IN_ARRAY1, int DIM1) { (double* lat, int len_lat), (double* lon, int len_lon) @@ -392,6 +432,12 @@ (int64_t* indices, int len, double* triangle_info_lats, int dmy1, double* triangle_info_lons, int dmy2) } +/* _expand_intervals +%apply (int64_t* in_array, int length, int resolution, int64_t* out_array1, int dmy1, int64_t* out_array2, int dmy2) { + (int64_t* indices, int len, int resolution, int64_t* range_indices, int len_ri, int64_t* result_size, int len_rs) +} + */ + %pythonprepend from_utc(int64_t*, int, int64_t*, int) %{ import numpy datetime = datetime.astype(numpy.int64) @@ -412,6 +458,13 @@ def to_compressed_range(indices): range_indices = range_indices[:endarg] return range_indices +def expand_intervals(intervals,resolution,result_size_limit=1000): + result = numpy.full([result_size_limit],-1,dtype=numpy.int64) + result_size = numpy.full([1],-1,dtype=numpy.int64) + _expand_intervals(intervals,resolution,result,result_size) + result = result[:result_size[0]] + return result + def to_hull_range(indices,resolution,range_size_limit=1000): out_length = range_size_limit range_indices = numpy.full([out_length],-1,dtype=numpy.int64) From 7ca063ccf3773f98eb0db8ea958822490af42006 Mon Sep 17 00:00:00 2001 From: Michael Rilee Date: Sat, 21 Sep 2019 20:49:43 -0400 Subject: [PATCH 31/31] Notes & version update. --- NOTES | 16 ++++++++++++---- VERSION | 2 +- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/NOTES b/NOTES index 5ad3f8c..4502eed 100644 --- a/NOTES +++ b/NOTES @@ -1,12 +1,20 @@ -* 2019-09-X Version 0.8.0 +* 2019-09-21 Version 0.8.0 Many improvements to pystare. - Vertices & Centroids - ConvexHull -- SpatialRange lists -Some use is made of the STARE spatial index as a way to convey -location information. +- SpatialRange interval lists +- Expand intervals into spatial id values +- Intersection +- Visualization examples + +Some use is made of the STARE spatial index as a way to convey location information, +however there is some difficulty going back and forth to SpatialVectors. Note +that vertices and edges are not the easiest points to determine. + +There was a bug in expandIntervals due to the confusing nature of +EmbeddedLevelNameEncoding with its native and SciDB formats. Added a number of tests using Constraint. diff --git a/VERSION b/VERSION index 5db90e1..8adc70f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.8.0-candidate \ No newline at end of file +0.8.0 \ No newline at end of file