From 8f27e65bddd7d4b8515ce620fb485fdd78fcdf89 Mon Sep 17 00:00:00 2001 From: Constantin Jucovschi Date: Fri, 24 Apr 2009 07:20:22 -0400 Subject: Initial commit --- rasodmg/interesttiling.cc | 809 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 809 insertions(+) create mode 100644 rasodmg/interesttiling.cc (limited to 'rasodmg/interesttiling.cc') diff --git a/rasodmg/interesttiling.cc b/rasodmg/interesttiling.cc new file mode 100644 index 0000000..aad8841 --- /dev/null +++ b/rasodmg/interesttiling.cc @@ -0,0 +1,809 @@ +/* +* This file is part of rasdaman community. +* +* Rasdaman community is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Rasdaman community is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with rasdaman community. If not, see . +* +* Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Peter Baumann / +rasdaman GmbH. +* +* For more information please see +* or contact Peter Baumann via . +/ +/** + * SOURCE: interesttiling.cc + * + * MODULE: rasodmg + * CLASS: r_Interest_Tiling + * + * COMMENTS: + * None +*/ + +#include "rasodmg/interesttiling.hh" +#include "rasodmg/alignedtiling.hh" +#include "rasodmg/dirdecompose.hh" +#include "rasodmg/dirtiling.hh" +#include "raslib/rminit.hh" + +#include +#include +#include +#include + +#ifdef AIX +#include +#endif + +// Uncoment the _VISUALIZE_2D_DECOMP_ line to generate ppm files the +// visualization of the domain decomposition done by the algoritm +// #define _VISUALIZE_2D_DECOMP_ + +// Uncoment the folling line for some debug printfs on this class +// #define _DEBUG_INTERESTTILING_ + +#ifdef _VISUALIZE_2D_DECOMP_ +#include "tools/visualtiling2d.hh" +#endif + +// This is a data structure for internal use within interesttiling.cc. +// It is defined here because of problems with ptrepository. +// ------------- +// This structure consists of a r_Minterval and a counter. Also, it has +// some auxiliary constructors. + +class Classified_Block + { + public: + int intersection_count; // Intersection count + r_Minterval block; // The actual block + + // Default constructor + Classified_Block(int count = 0); + + // Constructor with the actual block and the counter + Classified_Block(const r_Minterval& b, int count = 0); + + // Same data structure (this operator is needed because of dlist) + bool operator==(const Classified_Block& other) const; + + // Different data structure (this operator is needed because of dlist) + bool operator!=(const Classified_Block& other) const; + + // Friend of std::ostream (this operator is needed because of dlist) + friend std::ostream& operator<<(std::ostream& os, const Classified_Block& block); + }; + +const char* +r_Interest_Tiling::description = "dimensions, areas of interest, tile size (in bytes) and tile size limit [NOLIMIT|REGROUP|SUBTILING|REGROUPSUBTILING] (ex: \"2;[0:9,0:9];[100:109,0:9];100;REGROUPSUBTILING\")"; + +const char* r_Interest_Tiling::tilesizelimit_name_nolimit ="NOLIMIT"; +const char* r_Interest_Tiling::tilesizelimit_name_regroup ="REGROUP"; +const char* r_Interest_Tiling::tilesizelimit_name_subtiling ="SUBTILING"; +const char* r_Interest_Tiling::tilesizelimit_name_regroupandsubtiling ="REGROUPSUBTILING"; +const char* r_Interest_Tiling::all_tilesizelimit_names[r_Interest_Tiling::NUMBER]={ + tilesizelimit_name_nolimit, + tilesizelimit_name_regroup, + tilesizelimit_name_subtiling, + tilesizelimit_name_regroupandsubtiling + }; + +r_Interest_Tiling::Tilesize_Limit +r_Interest_Tiling::get_tilesize_limit_from_name(const char* name) +{ + + if(!name) { + RMInit::logOut << "r_Interest_Tiling::get_tilesize_limit_from_name(" << (name?name: "NULL") << ")" << endl; + return r_Interest_Tiling::NUMBER; + } + + unsigned int i=r_Interest_Tiling::NUMBER; + + for (i=0; i<(unsigned int)r_Interest_Tiling::NUMBER; i++) + { + if (strcasecmp(name, all_tilesizelimit_names[i]) == 0) + break; + } + return (r_Interest_Tiling::Tilesize_Limit)i; +} + +const char* +r_Interest_Tiling::get_name_from_tilesize_limit(Tilesize_Limit tsl) +{ + static const char* unknown="UNKNOWN"; + unsigned int idx = (unsigned int)tsl; + + if (idx >= (unsigned int)r_Interest_Tiling::NUMBER) + return unknown; + + return all_tilesizelimit_names[idx]; +} + +r_Interest_Tiling::r_Interest_Tiling(const char* encoded) throw (r_Error) + : r_Dimension_Tiling(0,0) +{ + if(!encoded) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << (encoded?encoded: "NULL") << ")" << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + std::vector vectInterestAreas; + r_Interest_Tiling::Tilesize_Limit tileSizeLimit; + r_Dimension tileD=0; + r_Bytes tileS=0, lenToConvert=0; + const char *pStart=NULL, *pRes=NULL, *pTemp=NULL, *pEnd=NULL; + char *pToConvert=NULL; + + + pStart=encoded; + pEnd=pStart+strlen(pStart); + pTemp=pStart; + pRes=strstr(pTemp,COLON); + + if(!pRes) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pTemp << "\"." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //deal with dimension + lenToConvert=pRes-pTemp; + pToConvert=new char[lenToConvert+1]; + memcpy(pToConvert,pTemp, lenToConvert); + pToConvert[lenToConvert]='\0'; + + tileD=strtol(pToConvert, (char**)NULL, DefaultBase); + if (!tileD) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pToConvert << "\"." << endl; + delete[] pToConvert; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + if (tileD < 0) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pToConvert << "\", is negative." << endl; + delete[] pToConvert; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //skip COLON && free buffer + delete[] pToConvert; + if(pRes != (pEnd-1)) + pRes++; + else + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, end of stream." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //parse interest areas + pTemp=pRes; + pRes=strstr(pTemp,COLON); + + if(!pRes) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + while(pRes) + { + //is interest areas? + if(*pTemp != *LSQRBRA) + break; + + //copy parsed interest area + lenToConvert=pRes-pTemp; + pToConvert=new char[lenToConvert+1]; + memcpy(pToConvert, pTemp, lenToConvert); + pToConvert[lenToConvert]='\0'; + + //add to vector + try + { + r_Minterval a(pToConvert); + vectInterestAreas.push_back(a); + } + catch(r_Error& err) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest area from \"" << pToConvert << "\"." << endl; + RMInit::logOut << "Error " << err.get_errorno() << " : " << err.what() << endl; + delete [] pToConvert; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //skip COLON + delete[] pToConvert; + if(pRes != (pEnd-1)) + pRes++; + else + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, end of stream." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //try next item + pTemp=pRes; + pRes=strstr(pTemp, COLON); + } + + if(vectInterestAreas.empty()) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, no interest areas specified." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //deal with tile size + lenToConvert=pRes-pTemp; + pToConvert=new char[lenToConvert+1]; + memcpy(pToConvert, pTemp, lenToConvert); + pToConvert[lenToConvert]='\0'; + + tileS=strtol(pToConvert, (char**)NULL, DefaultBase); + if (!tileS) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size from \"" << pToConvert << "\"." << endl; + delete[] pToConvert; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + if (tileS < 0) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size from \"" << pToConvert << "\", is negative." << endl; + delete[] pToConvert; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + //skip COLON && free buffer + delete[] pToConvert; + if(pRes != (pEnd-1)) + pRes++; + else + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size limit, end of stream." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + pTemp=pRes; + + tileSizeLimit=r_Interest_Tiling::get_tilesize_limit_from_name(pTemp); + if(tileSizeLimit==r_Interest_Tiling::NUMBER) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size limit from \"" << pTemp << "\"." << endl; + throw r_Error(TILINGPARAMETERNOTCORRECT); + } + + iareas = vectInterestAreas; + ts_strat = tileSizeLimit; + dimension = tileD; + tile_size = tileS; +} + +r_Interest_Tiling::r_Interest_Tiling(r_Dimension dim, const std::vector& interest_areas, r_Bytes ts, Tilesize_Limit strat) throw (r_Error) + : r_Dimension_Tiling(dim, ts), + iareas(interest_areas), + ts_strat(strat) +{ + for (std::vector::iterator it = iareas.begin(); it != iareas.end(); it++) + if (it->dimension() != dimension) + { + RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << dim << ", " << interest_areas << ", " << ts << ", " << strat << ") the interest area domain " << *it << " does not match the dimension of this tiling scheme (" << dimension << ")" << endl; + throw r_Edim_mismatch(dimension, it->dimension()); + } +} + +r_Interest_Tiling::~r_Interest_Tiling() +{ +} + +r_Tiling* r_Interest_Tiling::clone() const +{ + r_Tiling* copy = new r_Interest_Tiling(dimension, iareas, tile_size, ts_strat); + return copy; +} + +void r_Interest_Tiling::print_status(std::ostream& os) const +{ + os << "r_Interest_Tiling[ "; + r_Dimension_Tiling::print_status(os); + os << " interest areas = " << iareas << ", tiling strategy = " << ts_strat << " ]"; +} + +r_Tiling_Scheme +r_Interest_Tiling::get_tiling_scheme() const + { + return r_InterestTiling; + } + +static int r_Range_comp(const void *elem1, const void *elem2) +{ + r_Range e1 = *((r_Range*) elem1); + r_Range e2 = *((r_Range*) elem2); + + if (e1 == e2) + return 0; + + if (e1 < e2) + return -1; + else + return +1; +} + +std::vector* +r_Interest_Tiling::make_partition(const r_Minterval& domain) const +{ + r_Dimension dim = domain.dimension(); + int total = 2 * iareas.size(); + + // We need one decomp from each dimension + std::vector* part = new std::vector(dim); + + // We have at most (number of interest areas + 2) intervals + r_Range* intervals = new r_Range[total + 2]; + + // Create iterator for interest areas + std::vector::const_iterator it = iareas.begin(); + int total_iareas = iareas.size(); + + + // For all dimensions + for (r_Dimension i = 0; i < dim; i++) + { + it = iareas.begin(); // Reset iterator + intervals[0] = domain[i].low(); // Input lower domain limit + intervals[total+1] = domain[i].high(); // Input higher domain limit + + + for (int j = 1; j < total + 1; j += 2, ++it) // For all possible intervals + { + if ((*it)[i].low()-1 <= domain[i].low()) // Input low iarea limit + intervals[j] = domain[i].low(); + else + intervals[j] = (*it)[i].low()-1; + + intervals[j+1] = (*it)[i].high(); // Input higher iarea limit + } + + // Sort the table + qsort((void*) intervals, total+2, sizeof(r_Range), r_Range_comp); + + // Create partition using the limits table + for (int k=0; k* +r_Interest_Tiling::group(std::vector& blocks, r_Bytes typelen, Blocks_Type btype) const +{ + r_Bytes tilesize = get_tile_size(); + int joins = 0; + bool group_blocks = true; + + // The list of threated blocks + std::vector* treated = new std::vector; + + // An iterator for the blocks list + std::vector::iterator blocks_it = blocks.begin(); + + // The current block + r_Minterval current_block; + + // For all the blocks in list + while (!blocks.empty()) + { + // Get first block from list + current_block = blocks.back(); + blocks.pop_back(); + + //this is neccessary when the compiler optimizes out the .end() check + unsigned int numberOfLevels = blocks.size(); + for (blocks_it = blocks.begin(); blocks_it != blocks.end(); blocks_it++) + { + if (numberOfLevels == 0) + { + RMInit::logOut << "r_Interest_Tiling::group() the for loop was incorrectly optimized. breaking the loop." << endl; + break; + } + r_Minterval aux = *blocks_it; + + // In principle two blocks can't be merged + group_blocks = false; + + // If they can be merged + if (current_block.is_mergeable(aux)) + { + std::vector::iterator ia_it = blocks.begin(); + + switch (btype) + { + case BLOCKS_A: + + group_blocks = true; + + // Check if the two blocks belong exaclty to the same iareas + for (; ia_it != blocks.end(); ia_it++) + { + if (aux.intersects_with(*ia_it) != current_block.intersects_with(*ia_it)) + { + group_blocks = false; + break; + } + } + + break; + + case BLOCKS_B: + + for (; ia_it != blocks.end(); ia_it++) // For all iareas + { + // Find the one this block intersects + if (current_block.intersects_with(*ia_it)) + { + // Check if the other area intersects it + if (!aux.intersects_with(*ia_it)) + group_blocks = false; + else + group_blocks = true; + + break; + } + } + + if (!group_blocks) + break; + + group_blocks = false; + + case BLOCKS_C: // Falls in (this is, also applies to B); + + // Only on this two strategies, tilesize should be looked at + if ((ts_strat == REGROUP) || (ts_strat == REGROUP_AND_SUBTILING)) + { + // If the resulting size isn't larger than tilesize + if ((current_block.cell_count()+aux.cell_count()) * typelen + < get_tile_size()) + group_blocks = true; + } + else + group_blocks = true; + + break; + } + } + + // take care of the iterator advance, if is possible + if (group_blocks) + { + // Merge them + blocks.erase(blocks_it); // Automatically advances to next block + // take care of the size of the blocks + numberOfLevels--; + current_block.closure_with(aux); + ++joins; + } + else + numberOfLevels--; + } + + // Update the treated list with the current block + treated->push_back(current_block); + } + + std::vector* result = 0; + + // If there was joins, the algoritm must be repeted + if (joins != 0) + { + result = group(*treated, typelen, btype); + delete treated; + } + else + { + result = treated; + } + + // Return the result + return result; +} + +std::vector* +r_Interest_Tiling::compute_tiles(const r_Minterval& domain, r_Bytes typelen) const throw (r_Error) +{ + r_Dimension num_dims = domain.dimension(); // Dimensionality of dom +if (domain.dimension() != dimension) + { + RMInit::logOut << "r_Interest_Tiling::compute_tiles(" << domain << ", " << typelen << ") dimension (" << dimension << ") does not match dimension of object to tile (" << num_dims << ")" << endl; + throw r_Edim_mismatch(dimension, num_dims); + } +if (typelen > tile_size) + { + RMInit::logOut << "r_Interest_Tiling::compute_tiles(" << domain << ", " << typelen << ") tile size (" << tile_size << ") is smaller than type length (" << typelen << ")" << endl; + throw r_Error(TILESIZETOOSMALL); + } + +#ifdef _VISUALIZE_2D_DECOMP_ // User wants a visual + static int count; // Number of decomps + ++count; // Update num decomps + // of the 2D decomp. + Visual_Tiling_2D* vis; + if (domain.dimension() == 2) + { + // Create an object for visualization + char fname[80]; + sprintf(fname, "2D_decomp_int_%d.ppm", count); + vis = new Visual_Tiling_2D(domain, fname); + } + +#endif + + + // *** Main algoritm *** + + // The result + std::vector* result = new std::vector; + + // Create a partition for dir tiling + std::vector* part = make_partition(domain); + + // Perform dirtiling + r_Dir_Tiling *dir_tiling= NULL; + std::vector* dir_domain=NULL; + + try + { + dir_tiling=new r_Dir_Tiling(num_dims, *part, tile_size, r_Dir_Tiling::WITHOUT_SUBTILING); + + // Compute all tiles (directional tiles) + dir_domain = dir_tiling->compute_tiles(domain, typelen); + + delete dir_tiling; + } + catch(r_Error& err) + { + delete result; + delete part; + if(dir_tiling) + delete dir_tiling; + throw; + } + + // Get an interator for interest areas + //std::vector::iterator interest_area; + + // Create a list for holding the classifed blocks + std::vector part_domain; + + // Finds how many intersections exist between a block an the interest areas + for (std::vector::iterator dir_block = dir_domain->begin(); dir_block != dir_domain->end(); dir_block++) + { + Classified_Block b(*dir_block, 0); + + for (std::vector::const_iterator interest_area = iareas.begin(); interest_area != iareas.end(); interest_area++) + { + if (b.block.intersects_with(*interest_area)) + ++b.intersection_count; + } + + part_domain.push_back(b); + } + + // Lists used for grouping blocks + std::vector Out; + std::vector In_Unique; + std::vector In_Common; + + // Divide blocks into lists according to their number of intersections + for (std::vector::iterator class_block = part_domain.begin(); class_block != part_domain.end(); class_block++) + { + switch ((*class_block).intersection_count) + { + case 0: + Out.push_back((*class_block).block); + break; + case 1: + In_Unique.push_back((*class_block).block); + break; + default: + In_Common.push_back((*class_block).block); + break; + } + } + + + // Group blocks + std::vector* Blocks_A = group(In_Common, typelen, BLOCKS_A); + std::vector* Blocks_B = group(In_Unique, typelen, BLOCKS_B); + std::vector* Blocks_C = group(Out, typelen, BLOCKS_C); + + + std::vector::iterator it_A = Blocks_A->begin(); + std::vector::iterator it_B = Blocks_B->begin(); + std::vector::iterator it_C = Blocks_C->begin(); + + // If may be necessary to perform sub-tiling + if ((ts_strat == SUB_TILING) || (ts_strat == REGROUP_AND_SUBTILING)) + { + // Variable to hold result of sub-tiling + std::vector* subtiles = 0; + + // We need an array to hold the 3 iterators of the resulting blocks + std::vector::iterator* blocks_it[3]; + std::vector* blocks_vec[3]; + + // Put iterators on the list + blocks_it[0] = &it_A; + blocks_it[1] = &it_B; + blocks_it[2] = &it_C; + + blocks_vec[0] = Blocks_A; + blocks_vec[1] = Blocks_B; + blocks_vec[2] = Blocks_C; + + // For all the lists (Blocs_A, Blocks_B and Blocks_C) + for (int j=0; j<3; j++) + { + // Reset iterator to the begin + (*blocks_it[j]) = blocks_vec[j]->begin(); + + // Tile each block if necessary + for (; *blocks_it[j] != blocks_vec[j]->end(); (*blocks_it[j])++) + { + if ((**blocks_it[j]).cell_count()*typelen > get_tile_size()) + { + // Create a specification of a regular n-dim cube grid + r_Minterval specs(num_dims); + for (r_Dimension i=0; i::iterator it = subtiles->begin(); it != subtiles->end(); it++) + { + result->push_back(*it); + +#ifdef _VISUALIZE_2D_DECOMP_ + if (domain.dimension() == 2) + (*vis) << (*it); +#endif + } + + delete subtiles; + } + else // No subtiling needed + { + // Insert block as it is + result->push_back(**blocks_it[j]); + +#ifdef _VISUALIZE_2D_DECOMP_ + if (domain.dimension() == 2) + (*vis) << (**blocks_it[j]); +#endif + } + } + } + +#ifdef _VISUALIZE_2D_DECOMP_ + if (domain.dimension() == 2) + { + vis->set_pen(255, 255, 0); + + interest_area.seek_begin(); + while (interest_area.not_done()) + { + (*vis) << (*interest_area); + ++interest_area; + } + } +#endif + + } + else + { + // The result is just the sum of all blocks + + result->insert(result->end(), Blocks_A->begin(), Blocks_A->end()); + result->insert(result->end(), Blocks_B->begin(), Blocks_B->end()); + result->insert(result->end(), Blocks_C->begin(), Blocks_C->end()); + + // Visualization + +#ifdef _VISUALIZE_2D_DECOMP_ + if (domain.dimension() == 2) + { + vis->set_pen(255, 255, 0); + for (it_A.seek_begin(); it_A.not_done(); it_A++) + (*vis) << *it_A; + + vis->set_pen(0, 255, 0); + for (it_B.seek_begin(); it_B.not_done(); it_B++) + (*vis) << *it_B; + + vis->set_pen(0, 0, 255); + for (it_C.seek_begin(); it_C.not_done(); it_C++) + (*vis) << *it_C; + } +#endif + } + + // *** End of the main algorithm *** + + // *** Clean up *** + + delete part; + delete dir_domain; + delete Blocks_A; + delete Blocks_B; + delete Blocks_C; + +#ifdef _VISUALIZE_2D_DECOMP_ + if (domain.dimension() == 2) + delete vis; +#endif + + // Return result + + return result; +} + + + +std::ostream& operator<<(std::ostream& os, const Classified_Block block); + + +Classified_Block::Classified_Block(int count) + : intersection_count(count) + { + } + +Classified_Block::Classified_Block(const r_Minterval& b, int count) + : intersection_count(count), + block(b) + { + } + + +bool Classified_Block::operator==(const Classified_Block& other) const +{ + return (this->block == other.block) + && (this->intersection_count == other.intersection_count); +} + +bool Classified_Block::operator!=(const Classified_Block& other) const +{ + return !(*this == other); +} + + +std::ostream& operator<<(std::ostream& os, const Classified_Block block) +{ + os << "CBlock(" << block.intersection_count << "x: " << block.block << ")"; + + return os; +} -- cgit