openModeller
Version 1.4.0
|
00001 00028 #include <openmodeller/CallbackWrapper.hh> 00029 #include <openmodeller/Projector.hh> 00030 #include <openmodeller/Log.hh> 00031 #include <openmodeller/Environment.hh> 00032 #include <openmodeller/Model.hh> 00033 00034 #include <openmodeller/env_io/Header.hh> 00035 #include <openmodeller/env_io/Map.hh> 00036 #include <openmodeller/env_io/Raster.hh> 00037 #include <openmodeller/env_io/GeoTransform.hh> 00038 00039 #include <openmodeller/AreaStats.hh> 00040 00041 #include <openmodeller/Exceptions.hh> 00042 00043 #ifdef MPI_FOUND 00044 #include "mpi.h" 00045 #endif 00046 00047 #include <utility> 00048 using std::pair; 00049 00050 #ifdef MPI_FOUND 00051 #define N_X 10 00052 #define size_block 30000 00053 #endif 00054 00055 #ifndef MPI_FOUND 00056 /******************/ 00057 /*** create Map ***/ 00058 bool 00059 Projector::createMap( const Model& model, 00060 const EnvironmentPtr& env, 00061 Map *map, 00062 AreaStats *areaStats, 00063 CallbackWrapper *callbackWrapper ) 00064 { 00065 // Retrieve possible adjustments and/or additions made 00066 // on the effective header. 00067 Header hdr = map->getHeader(); 00068 00069 #ifndef GEO_TRANSFORMATIONS_OFF 00070 if ( ! hdr.hasProj() ) { 00071 00072 throw; 00073 } 00074 #endif 00075 00076 // Normalize the environment. 00077 model->setNormalization( env ); 00078 00079 if ( areaStats ) { 00080 00081 areaStats->reset( areaStats->getPredictionThreshold() ); 00082 } 00083 00084 MapIterator fin; 00085 MapIterator it = map->begin(); 00086 00087 int pixels = 0; 00088 int pixelcount = hdr.ydim * hdr.xdim; 00089 int pixelstep = pixelcount/20; 00090 bool abort = false; 00091 00092 Coord lg; 00093 Coord lt; 00094 Scalar val; 00095 00096 while ( it != fin ) { 00097 00098 // Call the abort callback function if it is set. 00099 if ( callbackWrapper && pixels%pixelstep == 0 ) { 00100 00101 try { 00102 00103 abort = callbackWrapper->abortionRequested(); 00104 00105 if ( abort ) { 00106 00107 Log::instance()->info( "Projection aborted." ); 00108 00109 if ( ! map->deleteRaster() ) { 00110 00111 Log::instance()->warn( "Could not delete map file." ); 00112 } 00113 00114 return false; 00115 } 00116 } 00117 catch( ... ) {} 00118 } 00119 00120 pair<Coord,Coord> lonlat = *it; 00121 00122 pixels++; 00123 ++it; 00124 00125 lg = lonlat.first; 00126 lt = lonlat.second; 00127 00128 Sample const &amb = env->get( lg, lt ); 00129 00130 // Read environmental values and find the output value. 00131 if ( amb.size() == 0 ) { 00132 00133 // Write noval on the map. 00134 map->put( lg, lt ); 00135 00136 val = -1; // could be used in a log 00137 } 00138 else { 00139 00140 val = model->getValue( amb ); 00141 00142 if ( val < 0.0 || val > 1.0 ) { 00143 00144 std::string msg = Log::format( "Suitability for point (%f, %f) is outside the range: %f", lg, lt, val ); 00145 throw AlgorithmException( msg.c_str() ); 00146 } 00147 00148 if ( areaStats ) { 00149 00150 areaStats->addPrediction( val ); 00151 } 00152 00153 // Write value on map. 00154 map->put( lg, lt, val ); 00155 } 00156 00157 // Call the callback function if it is set. 00158 if ( callbackWrapper && pixels%pixelstep == 0 ) { 00159 00160 float progress = pixels/(float)pixelcount; 00161 00162 if ( progress > 1.0 ) { 00163 00164 progress = 1.0; 00165 } 00166 00167 try { 00168 00169 callbackWrapper->notifyModelProjectionProgress( progress ); 00170 } 00171 catch( ... ) {} 00172 } 00173 } 00174 00175 // Call the callback function if it is set. 00176 if ( callbackWrapper ) { 00177 00178 try { 00179 00180 callbackWrapper->notifyModelProjectionProgress( 1.0 ); 00181 } 00182 catch ( ... ) {} 00183 } 00184 00185 map->finish(); 00186 00187 return true; 00188 } 00189 00190 #else 00191 /***************************/ 00192 /*** create Map Parallel ***/ 00193 bool 00194 Projector::createMap( const Model& model, 00195 const EnvironmentPtr& env, 00196 Map *map, 00197 AreaStats *areaStats, 00198 CallbackWrapper *callbackWrapper ) 00199 { 00200 00201 /*********************struct buff *************************************/ 00202 typedef struct 00203 { 00204 double lt; 00205 double lg; 00206 Scalar val; 00207 } buff; 00208 /**********************************************************************/ 00209 00210 buff X[N_X][size_block]; 00211 MPI_Status status; 00212 MPI_Request req_dado[N_X]; 00213 MPI_Request inf_N[N_X]; 00214 int ix; 00215 int n_elem[N_X]; 00216 int env_dado[N_X],total_pixels; 00217 int pixel_next,interval_pixels[3]; 00218 int count_tag,max_count_tag; 00219 int myrank,mysize; 00220 int j, N; 00221 00222 // Retrieve possible adjustments and/or additions made 00223 // on the effective header. 00224 MPI_Comm_rank(MPI_COMM_WORLD,&myrank); 00225 00226 Header hdr = map->getHeader(); 00227 00228 #ifndef GEO_TRANSFORMATIONS_OFF 00229 if ( ! hdr.hasProj() ) 00230 { 00231 throw; 00232 } 00233 #endif 00234 Log::instance()->debug( "Parallel version of createMap\n" ); 00235 // Normalize the environment 00236 model->setNormalization( env ); 00237 00238 if ( areaStats ) { 00239 areaStats->reset(); 00240 } 00241 00242 MapIterator fin; 00243 MapIterator my_it = map->begin(); 00244 00245 j = 0; //count the lines on struct array 00246 00247 MPI_Comm_rank(MPI_COMM_WORLD,&myrank); 00248 MPI_Comm_size(MPI_COMM_WORLD,&mysize); 00249 00250 int pixels, init_pixel, lim_pixel; 00251 count_tag=10; 00252 total_pixels=hdr.ydim*hdr.xdim; 00253 00254 MPI_Barrier(MPI_COMM_WORLD); 00255 00256 if (myrank==0) { //sif rank 0 00257 00258 max_count_tag=10+total_pixels/size_block; 00259 if ((total_pixels%size_block)!=0) 00260 max_count_tag++; 00261 00262 while (count_tag<max_count_tag) { 00263 MPI_Recv(&N,1,MPI_INT,MPI_ANY_SOURCE,count_tag,MPI_COMM_WORLD,&status); 00264 MPI_Recv(X[0],N*sizeof(buff),MPI_UNSIGNED_CHAR,status.MPI_SOURCE,count_tag,MPI_COMM_WORLD,&status); 00265 00266 for (j=0; j< N; j++) { 00267 if (X[0][j].val==100000) 00268 map->put( X[0][j].lg, X[0][j].lt); 00269 else { 00270 map->put(X[0][j].lg,X[0][j].lt,X[0][j].val ); 00271 } 00272 } 00273 count_tag++; 00274 } 00275 } //end if rank 0 00276 else 00277 if (myrank==mysize-1) { 00278 int k,req_pixels; 00279 pixel_next=0; 00280 interval_pixels[2]=10; //count_tag 00281 00282 while (pixel_next<(total_pixels)) { 00283 lim_pixel=pixel_next+size_block; 00284 00285 if (lim_pixel>total_pixels) 00286 lim_pixel=total_pixels; 00287 00288 interval_pixels[0]=pixel_next; 00289 interval_pixels[1]=lim_pixel; 00290 pixel_next=lim_pixel; 00291 MPI_Recv(&req_pixels,1,MPI_INT,MPI_ANY_SOURCE,1,MPI_COMM_WORLD,&status); 00292 MPI_Send(interval_pixels,3,MPI_INT,status.MPI_SOURCE,2,MPI_COMM_WORLD); 00293 interval_pixels[2]++; 00294 } 00295 00296 for (k=1;k<(mysize-1);k++) { 00297 interval_pixels[0]=total_pixels; 00298 MPI_Recv(&req_pixels,1,MPI_INT,MPI_ANY_SOURCE,1,MPI_COMM_WORLD,&status); 00299 MPI_Send(interval_pixels,2,MPI_INT,status.MPI_SOURCE,2,MPI_COMM_WORLD); 00300 } 00301 } 00302 else { 00303 ix=0; 00304 for (int k=0;k<N_X;k++) 00305 env_dado[k]=0; 00306 MPI_Send(&ix,1,MPI_INT,(mysize-1),1,MPI_COMM_WORLD); 00307 MPI_Recv(interval_pixels,3,MPI_INT,(mysize-1),2,MPI_COMM_WORLD,&status);; 00308 init_pixel=interval_pixels[0]; 00309 lim_pixel=interval_pixels[1]; 00310 00311 while (init_pixel<(total_pixels)) { 00312 my_it.nextblock(init_pixel); 00313 00314 j=0; 00315 00316 for (pixels=init_pixel; pixels<lim_pixel; pixels++) { 00317 00318 pair<Coord,Coord> lonlat = *my_it; 00319 ++my_it; 00320 Coord lg = lonlat.first; 00321 Coord lt = lonlat.second; 00322 00323 { //begin block 00324 Sample const &amb = env->get( lg, lt ); 00325 // TODO: use mask to check if pixel should contain prediction 00326 // Read environmental values and find the output value. 00327 if ( amb.size() == 0 ) { 00328 // Write noval on the map. 00329 00330 /* map->put( lg, lt );********************************************************/ 00331 /*substituir pela função, coloco lat, long, 255*******************************/ 00332 /*****escrevo nada na struct, 255*********************************************/ 00333 00334 X[ix][j].lt=lt; 00335 X[ix][j].lg=lg; 00336 X[ix][j].val=100000; 00337 00338 } 00339 else { 00340 Scalar val = model->getValue( amb ); 00341 00342 if ( val < 0.0 ) 00343 val = 0.0; 00344 else 00345 if ( val > 1.0 ) 00346 val = 1.0; 00347 00348 if ( areaStats ) { 00349 areaStats->addPrediction( val ); 00350 } 00351 // Write value on map. 00352 /********** map->put( lg, lt, val );substituir pela função***************/ 00353 X[ix][j].lt=lt; 00354 X[ix][j].lg=lg; 00355 X[ix][j].val=val; 00356 } 00357 } //end block 00358 j++; 00359 } //end do for 00360 00361 n_elem[ix]=j; 00362 00363 MPI_Isend(&n_elem[ix],1,MPI_INT,0,interval_pixels[2],MPI_COMM_WORLD,&inf_N[ix]); 00364 //MPI_Isend(&j,1,MPI_INT,0,interval_pixels[2],MPI_COMM_WORLD,&inf_N[ix]); 00365 MPI_Isend(X[ix],j*sizeof(buff),MPI_UNSIGNED_CHAR,0,interval_pixels[2],MPI_COMM_WORLD,&req_dado[ix]); 00366 env_dado[ix]=1; 00367 ix=(ix+1)%N_X; 00368 00369 if (env_dado[ix]==1) { 00370 MPI_Wait(&inf_N[ix],&status); 00371 MPI_Wait(&req_dado[ix],&status); 00372 env_dado[ix]=0; 00373 } 00374 00375 MPI_Send(&ix,1,MPI_INT,(mysize-1),1,MPI_COMM_WORLD); 00376 MPI_Recv(interval_pixels,3,MPI_INT,(mysize-1),2,MPI_COMM_WORLD,&status);; 00377 init_pixel=interval_pixels[0]; 00378 lim_pixel=interval_pixels[1]; 00379 } // end while 00380 } // end else myrank == 0 00381 00382 MPI_Barrier(MPI_COMM_WORLD); 00383 00384 // Call the callback function if it is set. 00385 00386 if ( callbackWrapper ) { 00387 00388 try { 00389 00390 callbackWrapper->notifyModelProjectionProgress( 1.0 ); 00391 } 00392 catch ( ... ) {} 00393 } 00394 00395 00396 return true; 00397 } 00398 #endif 00399