Tuning and hybrid parallelization of a genetic-based multi-point statistics simulation code

Tuning and Hybrid Parallelization of a Genetic-based Multi-Point Statistics Simulation Code Oscar Peredoa,b, Juli´ an M. Ortizb,c, Jos´ e R. Herrerod, Crist´ obal Samaniegoa aBarcelona Supercomputing Center (BSC-CNS), Department of Computer Applications in Science and Engineering, Edificio NEXUS I, Campus Nord UPC, Gran Capit´an 2-4, 08034 Barcelona, Catalunya, Spain. bAdvanced Laboratory for Geostatistical Supercomputing, Advanced Mining Technology Center, University of Chile. cDepartment of Mining Engineering, University of Chile. Av. Tupper 2069, Santiago, 837-0451, Chile. dComputer Architecture Department, Universitat Polit`ecnica de Catalunya (UPC-BarcelonaTech), Campus Nord, M`odul C6, Desp. 206, C/Jordi Girona 1-3, 08034 Barcelona, Catalunya, Spain. Abstract One of the main difficulties using multi-point statistical (MPS) simulation based on annealing techniques or genetic algorithms concerns the excessive amount of time and memory that must be spent in order to achieve convergence. In this work we propose code optimizations and parallelization schemes over a genetic-based MPS code with the aim of speeding up the execution time. The code optimizations involve the reduction of cache misses in the array accesses, avoid branching instructions and increase the locality of the accessed data. The hybrid parallelization scheme involves a fine-grain parallelization of loops using a shared-memory programming model (OpenMP) and a coarsegrain distribution of load among several computational nodes using a distributed-memory programming model (MPI). Convergence, execution time and speed-up results are presented using 2D training images of sizes 100×100×1 and 1000×1000×1 on a distributed-shared memory supercomputing facility. Keywords: Geostatistics, Stochastic simulation, Multi-point statistics, Code optimization, Parallel computing, Genetic algorithms 1. Multi-point statistics simulation1 Numerical modeling with geostatistical techniques aims at characterizing natural phenomena by summarizing and2 using the spatial correlation of collected data in order to measure the uncertainty at unsampled locations in space. As3 explained by Deutsch (2002), in simulation techniques, this spatial correlation is imposed into a model commonly con-4 structed on a regular lattice. The models must reproduce the statistical (histogram) and spatial distribution (variogram5 or other spatial statistics) and their quality is often judged in terms of the reproduction of geological features.6 Conventional techniques in geostatistics address the modeling using statistical measures of spatial correlation that7 quantify the expected dissimilarity (transition to a different category) between locations separated by a given vector8 distance, in reference to a given attribute, such as the facies, rock type, porosity, grade of an element of interest,9 among others. This is done using the variogram. Limitations of these techniques have been pointed out in that they10 only account for two locations at a time when defining the spatial structure (Krishnan and Journel (2003)). Much11 richer features can be captured by using multi-point statistics (MPS) that consider the simultaneous arrangement of12 the attribute of interest at several locations, providing the possibility to account for complex features, such as hierarchy13 between facies, delay effects, superposition or curvilinearity.14 MPS simulation aims at generating realizations that reproduce pattern statistics inferred from some training source,15 usually a training image. For example, in figure 1, left, we can see a training image based on sinuous channels with a16 simulated realization. These training images are used as a pattern database to generate simulations of the underlying17 image, as shown in figure 1, right. The simulations use those patterns with the aim that the training and simulated18 images share the same pattern histogram.19 There are several approaches to simulate accounting for MPS. Modifications of conventional methods to impose20 local directions of continuity using the variogram is a simple approach to impose some of the complex geological21 Preprint submitted to Parallel Computing December 12, 2013 Figure 1: Training image (left) and simulated realization (right) features (Xu (1996); Zanon (2004)). Object based methods and methods inspired in the genetic rules and physics22 of the deposition of sediments in different environments also seek to overcome the limitations of conventional cat-23 egorical simulation techniques, with significant progress (Deutsch and Wang (1996); Tjelmeland (1996); Pyrcz and24 Strebelle (2008)). Presently, the most popular method is a sequential approach based on Bayes’ postulate to infer25 the conditional distribution from the frequencies of multi-point arrangements obtained from a training image. This26 method, originally proposed by Guardiano and Srivastava (1993), and later efficiently implemented by Strebelle and27 Journel (2000), is called single normal equation simulation (snesim) (see also Strebelle (2002)). This method has28 been the foundation for many variants such as simulating directly full patterns (Arpat and Caers (2007); Eskandari29 and Srinivasan (2007)) and using filters to approximate the patterns (Zhang et al. (2006)). The use of a Gibbs Sampling30 algorithm to account directly for patterns has also been proposed (Boisvert et al. (2007); Lyster and Deutsch (2008)).31 A sequential method using a fixed search pattern and a ‘unilateral path’ also provides good results (Daly (2005); Daly32 and Knudby (2007); Parra and Ortiz (2009)). Other approaches available consider the use of neural networks (Caers33 and Journel (1998); Caers and Ma (2002)), updating conditional distributions with multi-point statistics as auxiliary34 information (Ortiz (2003); Ortiz and Deutsch (2004); Ortiz and Emery (2005)) or secondary variable (Hong et al.35 (2008)). Recently, a couple of new approaches focused on patching patterns directly to reduce computing time and36 impose larger scale structures, have been presented (Rezaee et al. (2013); Faucher et al. (2013)). These methods have37 a significant potential for practical applications. Alternatively, the problem can be addressed as an optimization one,38 using simulated annealing (Deutsch (1992)) or genetic algorithms (Peredo and Ortiz (2012)). The genetic approach is39 still under development, but essentially follows the same stochastic strategy as the annealing scheme. This work fo-40 cuses on code optimizations and parallelization of a genetic-based sequential code that simulates categorical variables41 to reproduce multi-point statistics. However, many of the techniques and ideas proposed here can be applied to other42 codes implementing similar simulation algorithms.43 In section 2 we explain the basic ideas about genetic algorithms, parallel architectures and programming models.44 After that, the main bottlenecks of the genetic-based simulation are detailed in section 3. A brief explanation of the45 actual implementation is presented in section 4, together with the proposed code optimizations and parallelization46 schemes, in sections 5 and 6 respectively. Finally, in the last sections we include the results obtained and final47 conclusions.48 2. Genetic algorithms and Parallel computing49 Genetic algorithms (GA) were developed in the 1970s with the work of Holland (1975) and in subsequent decades50 with De Jong (1980) and Goldberg (1989). Initially used to find good feasible solutions for combinatorial optimization51 problems, today they are used in various industrial applications, and recent advances in parallel computing have52 allowed their development and continuing expansion.53 In the canonical approach of GA, typically there is an initial population of individuals, where each individual is54 2 represented by a string of bits, as indivk=000110101, and a fitness function f itness(indivk) which represents the55 performance of each individual. The fitness function, or objective function, is the objective that must be minimized56 through the generations over all the individuals. A termination criteria must be defined in order to achieve the desired57 level of decrement in the fitness function. The main steps and operations performed in a canonical GA can be viewed58 in algorithm 1. Selection and restart are operations performed over the entire population. Selection extracts the best59 individuals and restart modifies part or the entire population in order to jump from local optimal values. Crossover60 and mutation are operations performed over particular individuals of the population. Crossover mixes the bits of two61 individuals according to a predefined set of cut points and mutation modifies specific random bits from one individual.62 Other operators can be found in the mentioned literature.63 Algorithm 1 Canonical genetic algorithm 1: INPUT:Nindividuals (population) 2: Evaluate a fitness function f itness in each individual 3: while termination criteria is not achieved do 4: {Breeding a new generation} 5: Sort the individuals by their fitness function value 6: if no improvement is measured in the population then 7: Restart: select some individuals and restart their bits 8: Sort the individuals by their fitness function value 9: end if 10: Selection: select the best individuals based on their fitness function 11: Crossover: breed new individuals crossing bits of individuals from the selection 12: Mutation: breed new individuals mutating some bits of individuals from the selection 13: Replace old individuals with new ones 14: Evaluate a fitness function f itness in each individual 15: end while 16: OUTPUT: best individual in the population In parallel computing architectures, as described by Culler et al. (1998), the two main models are distributed-64 memory and shared-memory, with their respective best known programming models MPI (Snir et al., 1998) and65 OpenMP (Chandra et al., 2001). In the first model, each processor has its own private memory and the data inter-66 changed between processors travels through a network in chunks of messages. The speed of this communication67 depends on the speed of the interconnection network. In the second model, each processor has access to a common68 memory through data coherence and data consistency methods.69 In order to use efficiently all the resources of parallel architectures, we need to explore algorithms that can exploit70 the parallelism and be able to adapt to future trends.71 Genetic algorithms receive the classification of embarrassingly parallel technique to solve problems. This clas-72 sification comes from the fact that separating the workload of the problem into several parallel tasks is trivial. This73 property motivates its investigation and application in the field of geostatistics, and particularly in MPS simulation.74 3. Bottlenecks of genetic-based MPS simulation75 We can see in algorithm 1 that the evaluation of the function f itness is performed #generations ×#individuals,76 and strong evidence indicates that this function is the most time consuming routine (a profiling using the gprof tool,77 Graham et al. (2004), tells us that for sufficiently large training images, more than 96% of the execution time is spent78 in this routine). Its calculation is based on an object called template. A template T consists in a set of coordinates that79 defines cell positions. Another interpretation is that a template is basically a pattern of memory accesses. Patterns(T)80 represents all the possible patterns that can be generated from a template Tgiven kpossible categories. If a template81 is defined as T={(1,1), (2,1), (1,2), (2,2)}(4 nodes) and the number of categories is k=2, the pattern database82 Patterns(T) will have 24elements. Complex geometries can be used to define the template and its corresponding83 pattern database, for example, in figure 2 we have a template defined as T={(1,1), (5,1), (9,1), (4,4), (5,4), (6,4),84 3 (1,5), (4,5), (5,5), (6,5), (9,5), (4,6), (5,6), (6,6), (1,9), (5,9), (9,9) }. This template is disconnected and its memory-85 access pattern is very irregular. This irregularity induces a slowdown in the overall performance when we need to86 traverse all its nodes. In contrast, the template of figure 3 has a regular memory-access pattern, if we access to the87 first element T(1,1), additionally the CPU will bring to the cache memory contiguous elements for free (each CPU88 has its own cache line size, for example 64 bytes, which means that each cache line can store 16 consecutive integers89 of 4 bytes). In Fortran, the CPU will bring a column line, in C/C++, a row line (see Hennessy and Patterson (1990)90 for more details about the CPU memory hierarchy). Figure 2: Template of 17 nodes with a complex geometry (irregular memory accesses) Figure 3: Template of 18 nodes with a simple geometry (regular memory accesses) 91 Given the pattern database Pattern(T), two main tasks must be performed:92 •First we have to count the frequency of appearances of each pattern in the training image and store them in an93 appropriate structure.94 •After that, for each individual in the population, and in each generation, we have to count the frequency of95 appearances of those patterns and calculate the following equation (other possible equations can be viewed in96 Peredo and Ortiz (2011)):97 f itness(indivk)=X p∈Patterns(T) Opf reqT I(p)−f reqindivk(p)2(1) 4 with Opa weight factor for each pattern, f reqT I(p) and f reqindivk(p) the number of appearances of pattern pin98 a training image and individual indivkrespectively.99 In both tasks, we have to handle with patterns located at boundary nodes. A buffer zone of halo nodes, with size100 equal to h=max{width(T),height(T)}, is added at the boundaries of the training image and realizations. For the101 training image, an extension of the original image is added keeping the geological continuity and statistical properties102 (alternative, for sufficiently large images, a reduction of size hin each side of the images can be applied, keeping the103 removed space as buffer zone). For each individual in the population, the buffer zone is filled with random categorical104 values between 0 and k−1.105 For the first task, following ideas from Straubhaar et al. (2011), we store the frequencies of the patterns that appear106 in the training image in a list L. In this list each element is a pair (d,f) where d=(s1,...,s|T|) is the pattern (data107 event), stored as an array of integers of length |T|with si∈ {0,...,k−1}(kcategories), and fis the frequency of108 appearance, stored as an integer. With this structures defined, the algorithm used in this task corresponds to algorithm109 2 using an empty frequency list L.110 For the second task, the accounting process is performed using algorithm 2. In this algorithm, we compare the Algorithm 2 Fitness function calculation: f itness(indivk) 1: INPUT: individual indivk, training image frequency list L, template T 2: sum ←0 3: Laux ← L 4: for each node (i,j) from indivkdo 5: Extract pattern located in node (i,j) using template Tand store it in an array localPattern of length |T| 6: Search localPattern in the list Laux 7: if localPattern exists in Laux then 8: Laux[localPattern]← Laux[localPattern]−1 9: end if 10: end for 11: for each pattern pin Laux do 12: sum ←sum +Laux[p]∗ Laux[p] 13: end for 14: OUTPUT:sum 111 histograms of the individual and the training image. A considerable bottle-neck for this calculation is the access to112 the list Lwhich stores the training image histogram, because for each extracted pattern from an individual, a search113 must be performed over it in order to see if this pattern exists in the histogram of the training image or not. A pro-114 posed solution to this problem is to store the elements (d,f) of the list using a lexicographical order in the patterns d.115 This order allows to search the existence of a pattern in the list using a binary search with an average and worst case116 performance of order O(log2n) comparisons, with nthe length of L.117 118 In the next section we will explain the implementation issues for this two tasks using code examples in Fortran 90119 as programming language, in order to see the data structures that are used and the proposed optimizations.120 4. Implementation121 4.1. Storage of pattern frequencies from training image122 The routines involved in the storage and management of the list Lare encapsulated in a Fortran module called123 patternOperations. This module consists of a set of global variables and routines performed over arrays of inte-124 gers. Part of its code structure is depicted in code 1.125 126 module patternOperations127 implicit none128 integer (4) :: npatterns129 5 type patternType130 integer(4), pointer :: pattern(:)131 integer (4) :: frequency132 end type patternType133 type( patternType ), pointer :: patternList (:)134 135 contains136 subroutine patternInsertion (...)137 subroutine patternSearch (...)138 subroutine patternComparison (...)139 subroutine printPatternList (...)140 end module patternOperations141 142143 Code 1: Module patternOperations.f90 Initially, when we store the pattern histogram of the training image, we use the global array patternList to144 keep track of the different patterns patternList(i)%pattern(1:tem nodes) with their respective frequencies145 patternList(i)%frequency. A first scan to the training image must be performed in order to fill patternList.146 Using the routine patternInsertion, each time a new pattern is found, the memory space used by patternList is147 re-allocated (adding one new element with frequency equal to 1) and all the elements previously inserted together with148 the new pattern are re-ordered according to a lexicographical order. If an existing pattern is found, the frequencies are149 updated. The lexicographical order is as follows: given two arrays of integers of the same length A=(a1,...,an) and150 B=(b1,...,bn) we will say that Ais greater than Bif and only if ∃k∈ {1,...,n}such that ak>bkand ∀isatisfying151 i<k,ai=biholds. For example, given this arrays A=(0,0,0,1,0,0) and B=(0,0,0,0,0,0), this order will indicate152 that Ais greater than B.153 After the filling step, patternList will be scanned each time we call the routine patternSearch. This routine154 gives the index position in the global array patternList where the searched pattern is located, if there is a match. If155 there is no match, it returns -1. This routine is implemented through a simple iterative binary search, which uses the156 routine patternComparison to compare patterns.157 The routine patternComparison basically traverses each pair of pattern’s nodes until they have different values158 and keeps track of which pattern has the greater one, returning 0 if they are equal, 1 if the first array is larger than the159 second one, and -1 otherwise. The parameters of this routine are two arrays with the pattern values and the length of160 those arrays (must be equal on both).161 As we explained before patternSearch and patternComparison are intensively used by the routine that cal-162 culates the fitness function. In the worst case, if the training image pattern list Lhas size nand the template has t163 nodes, a search in the list will perform t×O(log2n) comparisons. For very large templates, this pattern search is very164 time consuming and different data structures must be used in order to get a reasonable execution time in the search165 task.166 4.2. Calculate f itness(indivk)167 Following algorithm 2, which explains all the steps involved in the calculation of the fitness function, in this168 subsection we explain its associated routine presented in code 2. In this routine, the input indivk is a 2D integer169 array with an image loaded (an individual in the genetic algorithm terminology). The inputs tem nodes,tem rows,170 tem cols,tem coord rows and tem coord cols correspond to the specific template that we are using. For exam-171 ple, in figure 2, tem nodes =17, tem rows =tem cols =9 and tem coord rows and tem coord cols are arrays172 that store respectively the row and column coordinates of the nodes in template T.173 174 subroutine fitnessFunction (indivk , rows , cols ,175 tem_rows , tem_cols , tem_nodes ,176 tem_coord_rows ,tem_coord_cols ,value)177 use patternOperations !! contains global variables npatterns and patternList178 179 integer (4) , intent (in ) :: indivk (rows , cols )180 integer (4) , intent (in ) :: rows , cols181 integer (4) , intent (in ) :: tem_nodes182 integer (4) , intent (in ) :: tem_rows , tem_cols183 integer (4) , intent (in ) :: tem_coord_rows ( tem_nodes )184 6 integer (4) , intent (in ) :: tem_coord_cols ( tem_nodes )185 integer (4) , intent ( out ):: value186 integer (4) :: ii , irow ,icol ,inode ,pos187 integer (4) :: freq_aux ( npatterns )188 integer (4) :: localPattern ( tem_nodes )189 190 do ii =1, npatterns191 freq_aux ( ii )= patternList (ii)% frequency192 end do193 value =0194 do icol = 0,cols - tem_cols195 do irow = 0, rows - tem_rows196 do inode = 1, tem_nodes197 localPattern ( inode ) = indivk (198 irow + tem_coord_rows ( inode ) ,199 icol + tem_coord_cols ( inode )200 )201 end do202 call patter nSearch ( tem_nodes , localPattern , pos )203 if (pos /= -1) then204 freq_aux ( pos )= freq_aux ( pos ) -1205 end if206 end do207 end do208 do ii =1, npatterns209 value = value + freq_aux ( ii )* freq_aux ( ii )210 end do211 end subroutine fitnessFunction212 213214 Code 2: Subroutine fitnessFunction We can see that in this implementation, the values of Op, the weight factors described in equation (1), are equal to 0215 if f reqT I(p)=0 and equal to 1 otherwise (we only take into account patterns that are present in the training image).216 5. Code optimization217 We have applied several code optimization techniques to the previous described routines, in order to better exploit218 the CPU resources of the sequential execution. These optimizations can be grouped as: increase data locality, improve219 stack memory usage, code specialization of fitness routine, branch and load reductions.220 5.1. Increasing data locality of the main data structures221 The routines patternSearch and patternComparison are based on the global data structure patternList.222 The first modification consists in adapting this structure to the column-major order of the Fortran language in order to223 exploit the data and temporal locality.224 The column-major order in Fortran is related to the way in which the CPU accesses the data stored in memory.225 In this order, the matrices are accessed using the address row +(col −1) ∗numrows with numrows fixed. The cache226 lines that are moved from the main memory to the cache memory consist in contiguous memory addresses of fixed227 size. Leaving col fixed and traversing first all values of row, we can minimize the accesses to non-contiguous memory228 addresses (Hennessy and Patterson (1990)). Modifying the data structures in order to increase this kind of accesses229 reduces the cache data misses, reducing the overall execution time.230 The new structure is simply a 2D array in which the row size is the number of nodes in the template, tem nodes,231 and the column size is the number of patterns found in the training image.232 233 module patternOperations234 implicit none235 integer (4) :: npatterns236 integer (4) , pointer :: patternList (: ,:) !! ( tem_nodes ) X ( patterns in training image )237 integer (4) , pointer :: frequency (:) !! ( patterns in training image )238 239 7 contains240 subroutine patternInsertion (...)241 subroutine patternComparison (...)242 subroutine patternSearch (...)243 subroutine printPatternList (...)244 end module patternOperations245 246247 Code 3: Module patternOperations.f90 with re-designed data structures 5.2. Using the stack memory to store local arrays248 In Fortran, using the Intel’s compiler ifort, we can use the option -auto, which causes all local, non-SAVED249 variables to be allocated on the run-time stack, including fixed-length arrays. The default is -auto-scalar, saving250 all the scalar variables in the run-time stack. The main advantage is that the access time to the run-time stack is251 faster than the time access to the run-time heap (which stores the dynamically allocated memory), which decreases252 the execution time (Intel Corporation (2006)). A negative side of this option is related with the stack size. The stack253 has a maximum size fixed before the execution and if that size is exceeded, a stack overflow error can be obtained.254 Also, if we use several threads with OpenMP, each thread has its own stack memory space, so the amount of global255 stack space is increased proportionally to the number of threads. In order to avoid a stack overflow, we can calculate256 exactly how much data we need to allocate before compile time and see if it can fit in the stack. If it is too big, we257 can add more space to the stack, for example using the ulimit command in Linux operating systems or setting the258 environment variable OMP STACK SIZE with an appropriate value.259 5.3. Specialization of fitness function to an input template260 Given an input template, we can specialize our fitness function routine in order to exploit the data accesses pro-261 vided by the template. For example, in the fitness function of code 2 using the template of figure 2, the variables262 icol+1,...,icol+9 are calculated several times, but in reality we only need to calculate them one time per row-263 iteration and keep their values stored in an auxiliary variable. This allows for avoiding the accesses to the coordinate264 arrays tem coord rows and tem coord cols (see code 4) and keeping the cache memory clean for other data.265 If we denote by t=tem nodes,n=cols −tem cols and m=rows −tem rows, the total number of memory266 accesses performed by the fitness routine (only taking into account the arrays tem coord rows,tem coord cols267 and indivk) is 3 ×t×n×m. Using this optimization the total number of memory accesses for the same arrays268 is t×n×mwith a reduction of 3x less memory accesses than the original scenario. In the modified code depicted269 in 4, using the template described in figure 2, loop unrolling and common subexpression eliminations are included270 in order to eliminate the accesses of the arrays tem coord rows and tem coord cols and re-utilize the values of271 icol+1,...,icol+9.272 273 subroutine fitnessFunction (...)274 ...275 integer (4) :: rowplus1 , rowplus4 , rowplus5 , rowplus6 , rowplus9276 integer (4) :: colplus1 , colplus4 , colplus5 , colplus6 , colplus9277 ...278 do icol = 0, cols - tem_cols279 colplus1 = icol +1280 colplus4 = icol +4281 colplus5 = icol +5282 colplus6 = icol +6283 colplus9 = icol +9284 do irow = 0, rows - tem_rows285 rowplus1 = irow +1286 rowplus4 = irow +4287 rowplus5 = irow +5288 rowplus6 = irow +6289 rowplus9 = irow +9290 localPattern (1) = indivk ( rowplus1 , colplus1 )291 localPattern (2) = indivk ( rowplus5 , colplus1 )292 localPattern (3) = indivk ( rowplus9 , colplus1 )293 8 ...294 localPattern (15) = indivk ( rowplus1 , colplus9 )295 localPattern (16) = indivk ( rowplus5 , colplus9 )296 localPattern (17) = indivk ( rowplus9 , colplus9 )297 call patter nSearch ( tem_nodes , localPattern , pos )298 if (pos /= -1) then299 freq_aux ( pos )= freq_aux ( pos ) -1300 end if301 end do302 end do303 ...304 end subroutine fitnessFunction305 306307 Code 4: Subroutine fitnessFunction specialized to the input template of figure 2 5.4. Branch reduction308 The routine patternComparison, described in code 5, which is the most intensive in terms of execution time,309 can be re-designed in order to reduce the number of branch instructions (control flow, for example if-then-else or310 while loops) executed inside a loop.311 312 subroutine patternComparison ( length , onePattern1 , onePattern2 , value )313 !314 ! Compare two patterns .315 ! If value = 0, both patterns are equals316 ! If value = 1, the first pattern is bigger317 ! If value = -1, the second pattern is bigger318 !319 integer (4) , intent ( in ) :: length320 integer (4) , intent ( in ) :: onePattern1 ( length )321 integer (4) , intent ( in ) :: onePattern2 ( length )322 integer (4) , intent ( out ) :: value323 integer (4) :: ii324 325 value = 0326 ii = 0327 do while ( value == 0 . and . ii < length )328 ii = ii + 1329 if ( onePattern1 (ii) > onePattern2 (ii) ) then330 value = 1331 elseif ( onePattern1 ( ii) < onePattern2 ( ii ) ) then332 value = -1333 end if334 end do335 end subroutine patternComparison336 337338 Code 5: Subroutine patternComparison When a branch instruction is processed by the CPU, some cycles may be lost due to an incorrect branch prediction339 or an expensive condition evaluation (Hennessy and Patterson (1990)). For that reason, we know that reducing the340 number of branches and relaxing their boolean conditions are good practices in order to reduce the overall execution341 time. A first version of the modified routine can be viewed in code 6.342 343 subroutine patternComparison ( length , onePattern1 , onePattern2 , value )344 ! Compare two patterns .345 ! If value = 0, both patterns are equals346 ! If value = 1, the first pattern is bigger347 ! If value = 2, the second pattern is bigger348 integer (4) , intent ( in ) :: length349 integer (4) , intent ( in ) :: onePattern1 ( length )350 integer (4) , intent ( in ) :: onePattern2 ( length )351 integer (4) , intent ( out ) :: value352 9 100×100 1000×1000 Processes×Threads Time (seconds) Speed up Processes×Threads Time (seconds) Speed up 1×1 3.466 1x 1×1 312.133 1x 1×2 1.852 1.87x 1×12 55.62 5.61x 2×2 0.934 3.71x 2×12 27.53 11.33x 4×2 0.556 6.23x 4×12 14.03 22.24x 8×2 0.295 11.74x 8×12 7.26 42.99x 16×2 0.213 16.27x 16×12 4.02 77.64x Table 4: Coarse-grained parallelization without code optimizations: average generation step (30 generations) of the genetic algorithm using a population of 1000 individuals, each one a 100×100 and 1000×1000 images respectively, using the template in figure 2 100×100 1000×1000 Processes×Threads Time (seconds) Speed up Processes×Threads Time (seconds) Speed up 1×1 3.466 1x 1×1 312.133 1x 1×1+code-opt 2.032 1.70x 1×12+code-opt 42.911 7.10x 2×1+code-opt 1.027 3.37x 2×12+code-opt 22.102 14.12x 4×1+code-opt 0.561 6.17x 4×12+code-opt 11.354 27.49x 8×1+code-opt 0.310 11.18x 8×12+code-opt 5.916 52.76x 16×1+code-opt 0.221 15.68x 16×12+code-opt 3.085 101.17x Table 5: Coarse-grained parallelization with code optimizations: average generation step (30 generations) of the genetic algorithm using a population of 1000 individuals, each one a 100×100 and 1000×1000 images respectively, using the template in figure 2 According to the results from tables 2 and 3, the fine-grain parallelization allows to accelerate the execution590 considerably only when the size of the training image and realizations are large, in our case a size of 1000×1000591 reaching a speedup of 17.09x with 12 threads in the best case. With a size of 100×100 no profit from the parallelization592 was observed (except in the non-optimized code with 2 threads), due to the small amount of work that each thread593 performs compared with the overhead introduced by thread management. The best speedup result, 2.03x, was obtained594 using the optimized code with only one thread.595 The coarse-grain parallelization accelerates the execution proportionally to the number of processes involved in596 it. This feature is obtained after distributing the workload evenly among the processes. In our case, the workload is597 represented by the population, which is divided among the processes. According to tables 4 and 5, the best results598 achieved in terms of speedup were obtained using the hybrid parallelization scheme with optimized code, MPI and599 OpenMP on large training and realization images of size 1000×1000. The thread selection policy used to defined how600 many threads will run in each test was based in the best results obtained in the corresponding fine-grain tests from601 tables 2 and 3, namely, 2 and 12 threads for small and large images respectively, using non-optimized code, and 1 and602 12 threads for small and large images respectively, using optimized-code.603 7.2. Convergence604 In order to test the convergence of the method using the parallel optimized code, we choose a fixed set of param-605 eters for the genetic algorithm and with those parameters we observe the behaviour of the fitness function and the606 realization obtained, starting from a randomly generated population of individuals. Those parameters are:607 •population size: 640 individuals608 •mutation rate: described in the next paragraph609 •crossover percentage: 50%, corresponding to the number of individuals selected to perform a crossover with610 other individual611 •restart percentage: 10%, corresponding to the percentage of individuals that will be restarted, in each restart612 step613 •number of cut-points: 10%, corresponding to the number of cut-points in the crossover operation614 •number of mutation nodes: described in the next paragraph615 16 In order to accelerate the convergence, a multi-point mutation strategy was implemented. In this strategy, a616 random node and a random category were selected (in our tests, we have only 2 categories). Using an influence radius617 rcalculated a priori (in our case r=0.02 ×min(#rows,#columns)) all nodes that fall within the circle centered in618 the random node with radius rare mutated (changed) into a new selected category (leaving the conditionant nodes619 without modification). Additionally, a cyclic-cooling scheme was implemented in order to control the variability of620 the population, by means of reducing the mutation ratio by a factor λ∈(0,1) each 1000 generations. Starting from a621 mutation rate m=1, if λmgoes below a threshold, the ratio is restarted from m=1 and a new mutation cycle begins.622 Samples of realization images and convergence plots of function 1 are included in figures 6-7. The convergence623 plots show the relative decrease in percentage of the fitness function with respect to the initial value (the best fitness624 value obtained in the first population of individuals). As we can see in these figures, the simulated images are not equal625 to the training image, which is not bad, because one of the objectives of the simulation process using MPS methods626 is to obtain simulated images that fits the underlying statistics of the training image but not being necessarily equal.627 If conditional data is used (nodes with information from the training image which are not modified in the simulation628 process) the resulting simulated image will be more similar to the training image with the corresponding match of629 the underlying statistics. The tests that we are considering do not use any conditional data. Therefore they must be630 considered as a worst case scenario in terms of convergence.631 Figure 6: Simulated realization of size 100×100 (left) and convergence plot showing 17000 generations Figure 7: Simulated realization of size 1000×1000 (left) and convergence plot showing 17000 generations 8. Conclusions632 The proposed hybrid parallelization using an optimized code has shown reasonable speedup results, according to633 the time measurements reported in section 7. The differences between non-optimized and optimized code execution634 are considerable and justify our research in this recent topic. A negative aspect has to do with the loss of generality635 of some of the proposed optimizations. Applying routine specialization or branch/load reduction introduces several636 17 modifications in the corresponding code, were specific values of the size of the template are used. If the size or637 geometry of the template changes, we need to modify those routines in order to adapt them to the new template’s638 values. If the complete code has to be used in an agile framework allowing the utilization of user-defined templates,639 this behaviour can generate inconveniences because each new template needs its own optimized routines, with their640 corresponding creation, compilation and inclusion in the main application. A possible solution to this problem can be641 the adoption of auto-tuning techniques to automatically deploy new kernels according to the geometry of the template642 in use. Examples of auto-tuning techniques applied to stencil optimization for Finite Differences PDE solvers can be643 reviewed in Datta et al. (2008, 2009).644 The strategies proposed to accelerate the convergence allow us to get fast realizations, using reasonable small645 populations (40 individuals per process) and a small restart percentage (10% of the population is restarted). However,646 a further research on the acceleration of convergence by tuning the genetic parameters (population size, mutation647 and crossover rates, percentage of selection, percentage of restarted population, number of cross points and mutation648 points, among others) is left open for future research. Also, several other research topics can be explored, among649 them we can mention: dynamic mutation rate (using the full annealing scheme), non-linear crossovers (using external650 information to mix two realizations with some physical or geological interpretation), or selective restart (allowing to651 reset only individuals that meet certain properties).652 This work focuses on 2D training images and realizations. However 3D models are used in real geostatistical653 scenarios. In order to adapt our code to the 3D scenario, several further optimizations and modifications must be654 done, but these are left as future work. Among the most relevant ones we can mention the specialization of the655 fitness evaluation to 3D templates, use of efficient data structures to manage large 3D images and 1D individuals,656 new methods to perform crossover and mutations in these individuals, and modifications of the fine and coarse grain657 strategies to the new 3D scenario. Another related future work could be the application of the fitness calculation658 optimized routines into the simulated annealing scheme of simulation as described in Deutsch (1992), using the MPI659 implementation described in Peredo and Ortiz (2011) or the standard routines implemented in the GSLIB library from660 Deutsch and Journel (1992).661 Finally, another possible research area is related to explore new computer architectures with this algorithm. Among662 the possible alternatives are NVIDIA’s GPUs (using two programming models, CUDA and OpenACC), Intel’s MICs663 (several cores in one chip and accelerators working toghether in Intel architectures) or energy-efficient new supercom-664 puters that will be available in the next years.665 References666 Arpat, B., Caers, J., 2007. Stochastic simulation with patterns. Math. Geology 39, 177–203.667 Boisvert, J.B., Lyster, S., Deutsch, C.V., 2007. Constructing training images for veins and using them in multiple-point geostatistical simulation.668 33rd International Symposium on Application of Computers and Operations Research in the Mineral Industry, APCOM 2007, E. J. Magri (ed.)669 , 113–120.670 Caers, J., Journel, A.G., 1998. Stochastic reservoir simulation using neural networks trained on outcrop data. SPE Annual Technical Conference671 and Exhibition, New Orleans, LA, September 1998. Society of Petroleum Engineers. SPE paper # 49026 , 321–336.672 Caers, J., Ma, X., 2002. Modeling conditional distributions of facies from seismic using neural nets. Mathematical Geology 34, 143–167.673 Cant´ u-Paz, E., 1998. A survey of parallel genetic algorithms. Calculateurs paralleles, Reseaux et Systems Repartis 10.674 Chandra, R., Dagum, L., Kohr, D., Maydan, D., McDonald, J., Menon, R., 2001. Parallel programming in OpenMP. Morgan Kaufmann Pub. Inc.,675 San Francisco, CA, USA.676 Culler, D., Singh, J., Gupta, A., 1998. Parallel Computer Architecture: A Hardware/Software Approach. Morgan Kaufmann. 1st edition. The677 Morgan Kaufmann Series in Computer Architecture and Design.678 Daly, C., 2005. Higher order models using entropy, markov random fields and sequential simulation. Geostatistics Banff2004, Leuangthong, O.679 and Deutsch, C.V., eds., , Springer , 215–224.680 Daly, C., Knudby, C., 2007. Multipoint statistics in reservoir modelling and in computer vision. Petroleum Geostatistics A.681 Datta, K., Murphy, M., Volkov, V., Williams, S., Carter, J., Oliker, L., Patterson, D., Shalf, J., Yelick, K., 2008. Stencil computation optimization682 and auto-tuning on state-of-the-art multicore architectures, in: Proceedings of the 2008 ACM/IEEE conference on Supercomputing, IEEE Press,683 Piscataway, NJ, USA.684 Datta, K., Williams, S., Volkov, V., Carter, J., Oliker, L., Shalf, J., Yelick, K., 2009. Auto-tuning the 27-point stencil for multicore, in: In Proc.685 iWAPT2009: The Fourth International Workshop on Automatic Performance Tuning.686 De Jong, K.A., 1980. Adaptive system design: a genetic approach. IEEE Transactions on Systems, Man, and Cybernetics SMC-10, 566–574.687 Deutsch, C., Wang, L., 1996. Hierarchical object-based stochastic modeling of fluvial reservoirs. Mathematical Geology 28, 857–880.688 Deutsch, C.V., 1992. Annealing techniques applied to reservoir modeling and the integration of geological and engineering (well test) data. Doctoral689 dissertation. Ph.D. thesis. Stanford University.690 Deutsch, C.V., 2002. Geostatistical Reservoir Modeling. Oxford University Press, New York.691 18 Deutsch, C.V., Journel, A.G., 1992. GSLIB : geostatistical software library and user’s guide. Oxford University Press, New York, Oxford.692 Eskandari, K., Srinivasan, S., 2007. Growthsim - a multiple point framework for pattern simulation. Petroleum Geostatistics A.693 Faucher, C., Saucier, A., Marcotte, D., 2013. A new patchwork simulation method with control of the local-mean histogram. Stochastic Environ-694 mental Research and Risk Assessment 27, 1–21.695 Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley.696 Graham, S.L., Kessler, P.B., McKusick, M.K., 2004. gprof: a call graph execution profiler. SIGPLAN Not. 39, 49–57.697 Guardiano, F., Srivastava, M., 1993. Multivariate geostatistics: beyond bivariate moments. Geostatistics Troia 1, 133–144.698 Hennessy, J.L., Patterson, D.A., 1990. Computer Architecture: A Quantitative Approach. Morgan Kaufmann, San Mateo, CA.699 Holland, J., 1975. Adaptation in natural and artificial systems. U. of Michigan Press.700 Hong, S., Ortiz, J.M., Deutsch, C.V., 2008. Multivariate density estimation as an alternative to probabilistic combination schemes for data integra-701 tion. Geostats 2008 - Proceedings of the Eighth International Geostatistics Congress, J.M. Ortiz and X. Emery (eds.), Gecamin Ltda., Santiago,702 Chile 1, 197–206.703 Intel Corporation, 2006. Fortran Compiler Use of Temporaries. http://software.intel.com/file/20415 (Aug. 2012).704 Intel Corporation, 2011. Intel c Fortran Composer XE 2011. http://software.intel.com/en-us/articles/intel-composer-xe/ (Sep.705 2011).706 Krishnan, S., Journel, A.G., 2003. Spatial connectivity: From variograms to multiple-point measures. Mathematical Geology 35, 915–925.707 Labarta, J., Girona, S., Pillet, V., Cortes, T., Cela, J.M., 1995. A parallel program development environment, in: In proceedings of the 2th. Int.708 Euro-Par Conference, Springer. pp. 665–674.709 Lyster, S., Deutsch, C.V., 2008. Mps simulation in a gibbs sampler algorithm. Geostats 2008 - Proceedings of the Eighth International Geostatistics710 Congress, J.M. Ortiz and X. Emery (eds.), Gecamin Ltda., Santiago, Chile 1, 79–88.711 Ortiz, J.M., 2003. Characterization of high order correlation for enhanced indicator simulation, Unpublished doctoral dissertation. Ph.D. thesis.712 University of Alberta.713 Ortiz, J.M., Deutsch, C.V., 2004. Indicator simulation accounting for multiple-point statistics. Mathematical Geology 36, 545–565.714 Ortiz, J.M., Emery, X., 2005. Integrating multiple point statistics into sequential simulation algorithms. Geostatistics Banff2004, Leuangthong,715 O., and Deutsch, C.V., eds., Springer , 969–978.716 Parra, A., Ortiz, J.M., 2009. Conditional multiple-point simulation with a texture synthesis algorithm. IAMG 09 Conference, Stanford University .717 Peredo, O., Ortiz, J.M., 2011. Parallel implementation of simulated annealing to reproduce multiple-point statistics. Computers & Geosciences 37,718 1110 – 1121.719 Peredo, O., Ortiz, J.M., 2012. Multiple-point geostatistical simulation based on genetic algorithms implemented in a shared-memory supercom-720 puter, in: Abrahamsen, P., Hauge, R., Kolbjrnsen, O. (Eds.), Geostatistics Oslo 2012. Springer Netherlands. volume 17 of Quantitative Geology721 and Geostatistics, pp. 103–114.722 Pyrcz, M.J., Strebelle, S., 2008. A showcase of event-based geostatistical models, in: Ortiz, J.M., Emery, X. (Eds.), Geostas Oslo 2012. Gecamin723 Ltda., Santiago, Chile. volume 2, pp. 1143–1148.724 Rezaee, H., Mariethoz, G., Koneshloo, M., Omid Asghari, O., 2013. Multiple-point geostatistical simulation using the bunch-pasting direct725 sampling method. Computers & Geosciences 54, 293 – 308.726 Snir, M., Otto, S., Huss-Lederman, S., Walker, D., Dongarra, J., 1998. MPI-The Complete Reference, Volume 1: The MPI Core. MIT Press,727 Cambridge, MA, USA. 2nd. (revised) edition.728 Straubhaar, J., Renard, P., Mariethoz, G., Froidevaux, R., Besson, O., 2011. An improved parallel multiple-point algorithm using a list approach.729 Mathematical Geosciences , 1–24.730 Strebelle, S., 2002. Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology 34, 1–21.731 Strebelle, S., Journel, A.G., 2000. Sequential simulation drawing structures from training images. 6th International Geostatistics Congress, Cape732 Town, South Africa. Geostatistical Association of Southern Africa. .733 Tjelmeland, H., 1996. Stochastic models in reservoir characterization and Markov random fields for compact objects, Unpublished doctoral734 dissertation. Ph.D. thesis. Norwegian University of Science and Technology.735 Xu, W., 1996. Conditional curvilinear stochastic simulation using pixel-based algorithms. Mathematical Geology 28, 937–949.736 Zanon, S., 2004. Advanced aspects of sequential Gaussian simulation. Master’s thesis. University of Alberta.737 Zhang, T., Switzer, P., Journel, A., 2006. Filter-based classification of training image patterns for spatial simulation. Mathematical Geology 38,738 63–80.739 19