modules/saupePredictor.cpp

   1 #include "saupePredictor.h"
   2 using namespace std;
   3
   4 IStdEncPredictor::IOneRangePredictor* MSaupePredictor
   5 ::newPredictor(const NewPredictorData &data) {
   6 //      ensure the levelTrees vector is long enough
   7         int level= data.rangeBlock->level;
   8         if ( level >= (int)levelTrees.size() )
   9                 levelTrees.resize( level+1, (Tree*)0 );
  10 //      ensure the tree is built for the level
  11         Tree *tree= levelTrees[level];
  12         if (!tree)
  13                 tree= levelTrees[level]= createTree(data);
  14         ASSERT(tree);
  15
  16         int maxPredicts= (int)ceil(maxPredCoeff()*tree->count);
  17         if (maxPredicts<=0)
  18                 maxPredicts= 1;
  19
  20         OneRangePredictor *result=
  21                 new OneRangePredictor(data,settingsInt(ChunkSize),*tree,maxPredicts);
  22
  23
  24         #ifndef NDEBUG
  25         maxpred+= tree->count*(data.allowRotations?8:1)*(data.allowInversion?2:1);
  26         result->predicted= &predicted;
  27         #endif
  28
  29         return result;
  30 }
  31
  32 MSaupePredictor::Tree* MSaupePredictor::createTree(const NewPredictorData &data) {
  33 //      compute some accelerators
  34         const ISquareEncoder::LevelPoolInfos::value_type &poolInfos= *data.poolInfos;
  35         int domainCount= poolInfos.back().indexBegin;
  36         int level= data.rangeBlock->level;
  37         int sideLength= powers[level];
  38         int pixelCount= powers[2*level];
  39 //      create space for temporary domain pixels
  40         KDReal *domPix= new KDReal[ domainCount * pixelCount ];
  41 //      init domain-blocks' from every pool
  42         KDReal *domPixNow= domPix;
  43         int poolCount= data.pools->size();
  44         for (int poolID=0; poolID<poolCount; ++poolID) {
  45         //      check we are on the place we want to be; get the current pool, density, etc.
  46                 ASSERT( domPixNow-domPix == poolInfos[poolID].indexBegin*pixelCount );
  47                 const ISquareDomains::Pool &pool= (*data.pools)[poolID];
  48                 int density= poolInfos[poolID].density;
  49                 if (!density) // no domains in this pool for this level
  50                         continue;
  51                 int poolXend= density*getCountForDensity( pool.width, density, sideLength );
  52                 int poolYend= density*getCountForDensity( pool.height, density, sideLength );
  53         //      handle the domain block on [x0,y0] (for each in the pool)
  54                 for (int x0=0; x0<poolXend; x0+=density)
  55                         for (int y0=0; y0<poolYend; y0+=density) {
  56                         //      compute the average and standard deviation (data needed for normalization)
  57                                 Real sum, sum2;
  58                                 pool.getSums(x0,y0,x0+sideLength,y0+sideLength).unpack(sum,sum2);
  59                         //      the same as  = 1 / sqrt( sum2 - sqr(sum)/pixelCount  ) );
  60                                 Real multiply= 1 / sqrt( sum2 - sqr(ldexp(sum,-level)) );
  61                                 if ( !finite(multiply) )
  62                                         multiply= 1; // it would be better not to add the domains into the tree
  63                         //      if inversion is allowed and the first pixel is below zero then invert the block
  64                                 //if ( data.allowInversion && *domPixNow < 0 )
  65                                 //      multiply= -multiply;
  66                                 Real avg= ldexp(sum,-2*level);
  67                                 MatrixWalkers::AddMulCopy<Real> oper( -avg, multiply  );
  68                         //      copy every column of the domain's pixels (and normalize them on the way)
  69                                 using namespace FieldMath;
  70                                 KDReal *linCol= domPixNow;
  71                                 KDReal *linColEnd= domPixNow+pixelCount;
  72                                 for (int x=x0; linCol!=linColEnd; ++x,linCol+=sideLength)
  73                                         transform2( linCol, linCol+sideLength, pool.pixels[x]+y0, oper );
  74                         //      move to the data of the next domain
  75                                 domPixNow= linColEnd;
  76                         }
  77         }
  78         ASSERT( domPixNow-domPix == domainCount*pixelCount ); // check we are just at the end
  79 //      create the tree from obtained data
  80         Tree *result= Tree::Builder
  81                 ::makeTree( domPix, pixelCount, domainCount, &Tree::Builder::chooseApprox );
  82 //      clean up temporaries, return the tree
  83         delete[] domPix;
  84         return result;
  85 }
  86
  87 namespace NOSPACE {
  88         template<class T> struct AddMulCopyTo2nd {
  89                 T toAdd, toMul;
  90
  91                 AddMulCopyTo2nd(T add,T mul)
  92                 : toAdd(add), toMul(mul) {}
  93
  94                 template<class R1,class R2> void operator()(R1 f,R2 &res) const
  95                         { res= (f+toAdd)*toMul; }
  96                 void innerEnd() const {}
  97         };
  98
  99         struct SignChanger {
 100                 template<class R1,class R2> void operator()(R1 in,R2 &out) const
 101                         { out= -in; }
 102         };
 103 }
 104 MSaupePredictor::OneRangePredictor::OneRangePredictor
 105 ( const NewPredictorData &data, int chunkSize_, const Tree &tree, int maxPredicts )
 106 : chunkSize(chunkSize_), predsRemain(maxPredicts)
 107 , firstChunk(true), allowRotations(data.allowRotations), isRegular(data.isRegular)
 108 {
 109 //      compute some accelerators, allocate space for normalized range (+rotations,inversion)
 110         int rotationCount= allowRotations ? 8 : 1;
 111         heapCount= rotationCount * (data.allowInversion ? 2 : 1);
 112         points= new KDReal[tree.length*heapCount];
 113 //      if the block isn't regular, fill the space with NaNs (to be left on unused places)
 114         if (!isRegular)
 115                 fill( points, points+tree.length*heapCount, numeric_limits<KDReal>::quiet_NaN() );
 116 //      compute normalizing transformation and initialize a normalizator object
 117         Real multiply= sqrt(data.pixCount) / data.rnDev;
 118         AddMulCopyTo2nd<Real> colorNorm( -data.rSum/data.pixCount, multiply );
 119 //      compute SE-normalizing accelerator
 120         errorNorm.initialize(data);
 121         {
 122                 int sideLength= powers[data.rangeBlock->level];
 123                 Block localBlock(0,0,sideLength,sideLength);
 124         //      create normalized rotations
 125                 for (int rot=0; rot<rotationCount; ++rot) {
 126                 //      fake the matrix for this rotation
 127                         MatrixSlice<KDReal> currRangeMatrix;
 128                         currRangeMatrix.allocate( sideLength, sideLength, points+rot*tree.length );
 129                 //      fill it with normalized data
 130                         using namespace MatrixWalkers;
 131                         walkOperateCheckRotate( Checked<const SReal>(data.rangePixels,*data.rangeBlock)
 132                         , colorNorm, currRangeMatrix, localBlock, rot );
 133                 }
 134         }
 135 //      create inverse of the rotations if needed
 136         if (data.allowInversion) {
 137                 KDReal *pointsMiddle= points+tree.length*rotationCount;
 138                 FieldMath::transform2( points, pointsMiddle, pointsMiddle, SignChanger() );
 139         }
 140 //      create all the heaps and initialize their infos (and make a heap of the infos)
 141         heaps.reserve(heapCount);
 142         infoHeap.reserve(heapCount);
 143         for (int i=0; i<heapCount; ++i) {
 144                 PointHeap *heap= new PointHeap( tree, points+i*tree.length, !data.isRegular );
 145                 heaps.push_back(heap);
 146                 infoHeap.push_back(HeapInfo( i, heap->getTopSE() ));
 147         }
 148 //      build the heap from heap-informations
 149         make_heap( infoHeap.begin(), infoHeap.end() );
 150 }
 151
 152 MSaupePredictor::Predictions& MSaupePredictor::OneRangePredictor
 153 ::getChunk(float maxPredictedSE,Predictions &store) {
 154         ASSERT( PtrInt(heaps.size())==heapCount && PtrInt(infoHeap.size())<=heapCount );
 155         if ( infoHeap.empty() || predsRemain<=0 ) {
 156                 store.clear();
 157                 return store;
 158         }
 159 //      get the number of predictions to make (may be larger for the first chunk)
 160         int predCount= chunkSize;
 161         if (firstChunk) {
 162                 firstChunk= false;
 163                 if (heapCount>predCount)
 164                         predCount= heapCount;
 165         }
 166 //      check the limit for prediction count
 167         if (predCount>predsRemain)
 168                 predCount= predsRemain;
 169         predsRemain-= predCount;
 170 //      compute the max. normalized SE to predict
 171         float maxNormalizedSE= errorNorm.normSE(maxPredictedSE);
 172 //      make a local working copy for the result (the prediction), adjust its size
 173         Predictions result;
 174         swap(result,store); // swapping is the quickest way
 175         result.resize(predCount);
 176 //      generate the predictions
 177         for (Predictions::iterator it=result.begin(); it!=result.end(); ++it) {
 178                 pop_heap( infoHeap.begin(), infoHeap.end() );
 179                 HeapInfo &bestInfo= infoHeap.back();
 180         //      if the error is too high, cut the vector and exit the cycle
 181                 if ( bestInfo.bestError > maxNormalizedSE ) {
 182                         result.erase( it, result.end() );
 183                         infoHeap.clear(); // to be able to exit more quickly in the next call
 184                         break;
 185                 }
 186         //      fill the prediction and pop the heap
 187                 ASSERT( 0<=bestInfo.index && bestInfo.index<heapCount );
 188                 PointHeap &bestHeap= *heaps[bestInfo.index];
 189                 it->domainID= isRegular
 190                         ? bestHeap.popLeaf<false>()
 191                         : bestHeap.popLeaf<true>();
 192                 it->rotation= allowRotations ? bestInfo.index%8 : 0; // modulo - for the case of inversion
 193         //      check for emptying the heap
 194                 if ( !bestHeap.isEmpty() ) {
 195                 //      rebuild the infoHeap heap
 196                         bestInfo.bestError= bestHeap.getTopSE();
 197                         push_heap( infoHeap.begin(), infoHeap.end() );
 198                 } else { // just emptied a heap
 199                         infoHeap.pop_back();
 200                 //      check for emptying the last heap
 201                         if ( infoHeap.empty() )
 202                                 break;
 203                 }
 204         }
 205 //      return the result
 206         swap(result,store);
 207
 208         #ifndef NDEBUG
 209         *predicted+= store.size();
 210         #endif
 211
 212         return store;
 213 }