An Implementation of Markov Random Fields - Index of
Contents
1. lt steps gt 10 lt steps gt real number in 0 1 for annealing lt annealing gt schedule lt steps gt positive integer grid updates per temperature update 12 ME An Implementation of Markov Random Fields Annealing is useful mainly for MCMC simulations If the simulation method is set to sequen tial and annealing is nevertheless used then this simulation will be performed at the temperature Tow i 4 4 Grid size The grid is assumed to be cell centered Cartesian The parameters to be specified describe grid size in x y and z direction Parameters in model file Description Details Valid values lt grid gt Grid size the number of cells in xz lt nx gt 1 2 3 lt nx gt 100 lt nx gt y and z direction lt ny gt 100 lt ny gt lt nz gt 1 lt nz gt lt grid gt 4 5 Neighbourhood In a Markov random field simulation setting it may be convenient to think about the doublet neighbourhood interaction potential in two slightly different ways In the first of these the interaction potential and neighbourhood are defined independently Think of this as a situation where there is a physical system with interactions strong or weak at arbitrary inter particle lags Superpositioned on this we impose a neighbourhood structure The resulting MRF model consists of all interactions for which the interacting particles are within one anothers neighbourhood All other interactions are set to zero2. B Parsing Xml files with NRIib An Implementation of Markov Random Fields J o oNN O O O O 10 10 12 12 13 13 14 16 16 17 17 18 18 18 19 19 21 22 23 23 23 24 26 1 Introduction The project of which this report is a part aims at comparing and further develop different meth ods for multipoint statistics Multipoint methods have been introduced as a tool for modelling geological structures and they seem to provide better representations than classical variogram methods Markov random fields constitute one type of multipoint methods To a large extent a comparison between Markov random fields and other multipoint methods has to be based on simulations In particular this is true when comparing how well the differ ent methods can reproduce patterns and essential characteristics of a given training image It is therefore of great importance to establish a suitable software implementation of MRFs We have implemented a program that uses traditional Markov chain Monte Carlo MCMC methods for simulating Markov random fields Sequential simulation has been implemented as an optional choice although the implementation for this has not yet been developed as far as the MCMC implementation The implementation is done in C 2 MRF the implemented models This section provides the definition of a Markov random field MRF as well as a basic expression for its joint probability distribution The section a
3. lhe 13 Ayi E Ly Ly 1 1 0 1 Ly 1 Ly 14 Azij Lz L 1 1 0 1 Lz 1 Lz 15 The positions are listed according to the rule looping over x is done before looping over y is done before looping over z The first distance for which a potential value is listed is Ax Ayij Azi lig Ly Lz For each distance there is a list of K elements where K is the number of facies Each of these elements corresponds to a two particle facies configuration z zj ki kj The order of these configurations is assumed to correspond to the nested loops for kl 0 K 1 for k2 0 K 1 list potential value G end end Combining the distance and facies conventions the file should list the potential according to the order of the following nested loops for z Lz Lz for y Ly Ly for x Lx Lx for ki 0 K 1 for k2 0 K 1 list potential value G end end end end end An external field can also be included in the input file It should be listed as the potential value for z y z O and kz k 1 in the loop above This format of the input file is useful if the interaction potential is a known mathematical function In this case the file is easily generated by using the nested for loops described above If the potential is not a known mathematical function it may be easier or more convenient to generate the file on the format used by the
4. lt type gt potts lt type gt lt beta gt 0 5 lt beta gt lt input file gt filename dat lt input file gt lt interaction gt With lt type gt general lt type gt an input file must be provided lt beta gt is ignored and the neighbourhood is implicitly defined through an input file With lt type gt two particle lt type gt an input file must be provided and lt beta gt is ignored With lt type gt generalized potts lt potts gt an input file for spatial part must be provided lt beta gt determines the strength of the facies dependency of the interaction With lt type gt potts lt potts gt any input file for spatial part will be ignored lt beta gt determines the strength of the facies dependency of the interaction File formats are described in Section 5 lt type gt general two particle generalized potts or potts lt beta gt R lt input file gt Any existing text file specific format 2 The relationship between the sequential probability distribution and p z z_ is discussed in detail in the report 2 An Implementation of Markov Random Fields ME 15 The path of the input potential file should be absolute or relative to the folder lt topdir gt de scribed in Section 4 9 4 7 Lithology The number of facies to be used in the simulation is communicated to the program by listing arbitrary facies names The interpretation of these names in the program is tha
5. for convergence checking Default name measurementlog txt lt result gt File where the final realization of the grid is written when simulation is over Default name result dat lt temp gt Files where the temporary grid configuration is periodically written Each writ ing is done to a unique file that is automatically labelled with a number by the program The number is prepended to the given file name Default name tmpresult dat Log frequencies lt writefreq gt The number of iterations between each measurement and writing of mean facies lt gridfreq gt The number of iterations between each writing of instantaneous grid config uration Parameters in model file Valid values lt output gt lt measurements gt measurementlog txt Any file names may be used but lt measurements gt keep the file extensions as they are lt result format sgems gt result dat lt result gt here txt dat dat Formats may lt temp format sgems gt tmpresult dat lt temp gt be sgems or storm lt gridfreq gt 500000 lt gridfreq gt lt gridfreq gt lt writefreq gt 1 2 lt writefreq gt 10000 lt writefreq gt lt output gt 4 9 File directories Two names specify to which folder the simulation results are stored The parameters holding these two names are called lt topdir gt and lt prefix gt The folder with the name given by lt prefix gt is a subfolder under the folder
6. in the MRF model This way of thinking about the doublet is convenient for exploring how a change in the neighbourhood structure affects the results of the MRF simulation provided the underlying interactions are fixed One example is to study the effect of neglecting interactions at long lags in a two particle interaction model The other way to think about the doublet is to let the neighbourhood be implicitely defined by the interaction potential That is if two particles interact via a non zero potential then they are per definition in one anothers neighbourhood These two ways of thinking about the doublet are of course completely interchangable in the sense that any model that is set up using one point of view can easily be set up using also the other point of view Which viewpoint that is seen as the most natural or intuitive is merely a matter of taste and has no effect on the simulated results Our implementation of Markov random fields supports explicitely both these points of view provided that the interactions are limited to two particle interactions For general interaction models only the second point of view is supported The different interaction models are described in Section 4 6 This section describes the parameters that are to be specified if the neighbourhood is defined independently of the interaction potential i e adopting the first viewpoint described above Two neighbourhood types are implemented specified by the paramet
7. interaction mode general i e the format described in Section 5 1 If so the simulation itself must of course be run using the general mode 5 3 Generalized Pott s interactions If the interaction type is generalized potts the expression for the interaction is according to Eq 8 written Fi 24 ziyi b zizi fi 22 ME An Implementation of Markov Random Fields In this case the input file is only assumed to list the spatial dependency f Apart from this restric tion the convention for the order of the elements is analogous to what was said for the general two particle interaction in the previous section Hence the order in which the file should list the spatial potential is for z Lz Lz for y Ly Ly for x Lx Lx list potential value f end end end In addition at the very beginning of the file the three integers Lz Ly L must be given These are used by the program to identify the index of the correct file element for a given lag l 5 4 Closing remarks on interaction models For a new user of the program it is probably most convenient to stick to the interaction model general and the corresponding format of input potential file This interaction model has the sim plest file format and as it can handle all cases including the two particle specific cases described earlier it may be a useful attitude to concentrate on this kind of interaction model We also expect that further extensions of
8. is 1 Parameters in model file Description Details Valid values lt lithology gt Names and target global volume lt facies gt fractions lt name gt Any text lt name gt sand lt name gt lt name gt and lt fraction gt must be Target values lt fraction gt 0 7 lt fraction gt specified for each facies lt fraction gt should lt facies gt be R and must lt facies gt sum to 1 lt name gt shale lt name gt lt fraction gt 0 3 lt fraction gt lt facies gt lt lithology gt lt fraction control gt Fraction control The parameter S R 1 that determines the change in Recommended lt fraction control gt external field given that target value is 1 Set to 0 value and measured volume if fraction control fraction are unequal is not to be used 4 8 Output Several output files are automatically generated by the simulation program These have default names that may be changed if desired In addition there are two parameters describing the fre quency of logging 3 For correct interpretation of simulation results the input potentials of Section 4 6 should reflect this ordering in the sense that z 0 means the first facies listed z 1 the next and so on 16 ME An Implementation of Markov Random Fields Output files There are three kinds of log files lt measurements gt File where the instantaneous mean facies in the grid is periodically appended throughout the simulation May be used
9. pointer to a particular element in the xml file Below we list some illustrative parts of the Parse function in the FaciesSimulation class The idea is to parse the whole file and then report errors or warnings right before we return To this end we use a string variable to collect error and warning messages and a status variable to keep track of the keyword status for each keyword as well as the overall parsing status For direct decendants of the action element access of individual input data is done by use of Xm1Parser GetChildArgument or XmlParser RequireChildArgument depending on whether the element is optional or necessary Optional elements have a default value In this example we assume that default values and list of allowed values are defined in a Definitions class Let us first try to read in the annealing keywords lt annealing gt lt max temperature gt lt max temperature gt lt temperature reduction factor gt 0 9 lt temperature reduction factor gt lt steps gt 10 lt steps gt lt annealing gt 26 INRES An Implementation of Markov Random Fields 11 13 15 17 19 21 23 25 27 29 31 33 35 39 const Xml Element annealing element element FirstChildElement annealing if annealing element kw status XmlParser GetChildArgument annealing element max temperature annealing temperature errors Default GetAnnealingTemperature status std max stat
10. the program will focus on using this model and then the interaction models two particle generalized potts and potts will phase out 6 Conclusion Markov random fields have been implemented in a C framework called MRF User input is taken from an Xml model file In this report we have documented the overall structure of the program the syntax and parsing of the model file with the help of NRlib Hence the document is both a user guide and a programmers guide The MRF software is a research tool under development and features and specifications may change without notice References 1 Harald H Soleng NRlib A C toolbox Note SAND 10 05 Norwegian Computing Center Oslo Norway November 2005 2 Heidi Kjonsberg and Ingeborg Ligaarden Second order markov mesh models described as markov random fields Note SAND 07 06 Norwegian Computing Center Oslo Norway 2006 An Implementation of Markov Random Fields NE 23 lt A Example of model file lt xml version 1 0 encoding IS0 8859 1 7 gt lt xml stylesheet type text css href http intern nr no sand css cohiba css gt lt mrf version 0 1 gt lt title gt Example of MRF model file lt title gt lt Changed by Harald H Soleng 25 Sep 2006 gt lt name gt Test model lt name gt lt messages gt lt logfile gt lt level gt 5 lt level gt lt name gt logs mrf log lt name gt lt logfile gt lt screen gt lt level gt 5 l
11. 41 P zizizi fi 8 Here 0 non spatial factor of the interaction and f is a function that determines spatial dependence 2 4 IS the Kronecker delta function 8 is a parameter that gives the strength of the That is if the two cells have the same facies their interaction energy is p times a spatial function that depends on their relative position If the two cells have different facies their interaction energy is zero We refere to the special case of f 1 as the simple Pott s model and f 1 as the generalized Pott s model The most general model above number 1 in the list does of course include the special cases listed as second and third The reason why we have implemented the special cases of general two particle interactions and Pott s model simple and generalized separately is partly due to the historic development of the program partly to speed up the simulations whenever possible Input parameters to the simulation program determines which model to consider This is further explained in Section 3 3 Program structure This section describes the implementation of the Markov random field The implementation is done in C with extensive use of the in house proprietary library NRlib 1 3 1 Environment A model file provides input parameters to the simulation with the model file being in xml format The model file specifies all necessary parameters except the exact form of the Markov random field potential which in mo
12. An Implementation of Markov Random Fields Overview and user s manual for version 0 1 www photos com 2006 SAND 06 06 Heidi Kjonsberg Harald H Soleng 17th October 2006 Norsk Regnesentral NORWEGIAN COMPUTING CENTER Norwegian Computing Center Norsk Regnesentral Norwegian Computing Center NR is a private independent non profit foundation established in 1952 NR carries out contract research and development projects in the areas of information and communication technology and applied statistical modeling The clients are a broad range of industrial commercial and public service organizations in the national as well as the international market Our scientific and technical capabilities are further developed in co operation with The Research Council of Norway and key customers The results of our projects may take the form of reports software prototypes and short courses A proof of the confidence and appreciation our clients have for us is given by the fact that most of our new contracts are signed with previous customers Norsk Regnesentral Bes ksadresse Telefon telephone Internett internet Norwegian Computing Center Office address 47 22 85 25 00 Www nr no Postboks 114 Blindern Gaustadall en 23 Telefaks telefax E post e mail NO 0314 Oslo Norway NO 0373 Oslo Norway 47 22 69 76 60 nr nr no Title An Implementation of Markov Random Fields Authors Heidi Kjonsberg Harald H Soleng Date 17th Octob
13. configuration for MCMC simula constant lt initialize gt tions lt realizations gt 1 lt realizations gt The desired number of indepen 1 2 dent realizations 4 2 Path The path specifies in which order the grid cells are visited It can be sequential cell number 1 2 3 or random If a random path is chosen the MCMC simulation will use a reshuffled path for each grid update For sequential simulations each cell is visited only once using the specified path to run through the grid Parameters in model file Description Details Valid values lt path gt sequential lt path gt Specifies the path to be used sequential random 4 3 Annealing Annealing can be used in the MCMC simulations It slowly brings the system down from a start ing temperature Tmaz to the default temperature T 1 The temperature develops in steps according to the formula T 41 Tr 1 1 Three parameters must be specified Tmax Starting temperature nr the number of grid updates per temperature determines how much the temperature is decreased in each step see the formula above Parameters in model file Valid values lt annealing gt lt max temperature gt non negative lt max temperature gt 10 lt max temperature gt real number Zero means no lt temperature reduction factor gt 0 9 annealing lt temperature reduction factor gt lt temperature reduction factor gt
14. d by the imple mented program For this kind of simulation the initial grid is empty and the cells are se quentially assigned a facies throughout the simulation Each cell in the grid is visited only once and the assignment of facies is based on the conditional probability p z z where ni is the so called sequential neighbourhood of cell i Our implementation allows the user to specify 7 see Section 4 5 for further details The probability p z z is assumed to have a form similar to the two particle version of Section 2 3 This is further explained in Section 4 6 lt iterations gt The number of iterations per cell Will be reset to 1 if sequential simulation is chosen lt initialize gt lt type gt The parameter lt type gt determines the start configuration of the grid in the case of MCMC simulations If specified to random a random initial configuration is 10 ME An Implementation of Markov Random Fields Timer amp timer Singleton lt Timer gt Instance Program amp program Singleton lt Pro let dram Instance program Initialize MRF 1 Markov Random Fields Runs a MRF simulation _ 2005 _ Norwegian Computing Center Options options Singleton lt Options gt Instance options Initialize options Read argc argv const std list lt std string gt args options GetArguments std string model_file args size gt 0 model file args fro
15. e 0 for v Definitions GetSimulationModes begin 97 v Definitions GetSimulationModes end v mode VE ev exe 99 break 101 cast to enum 103 method static_cast lt Definitions SimulationMode gt mode Consider the result keyword below Assume that it is required by the program lt output gt lt measurements gt measurementlog txt lt measurements gt lt result format sgems gt result dat lt result gt lt temp format sgems gt tmpresult dat lt temp gt lt gridfreq gt 500000 lt gridfreq gt lt writefreq gt 10000 lt writefreq gt lt output gt In order to read it we first have to locate the output element const Xml Element output element wo element FirstChildElement output if output_element w errors std string ERROR output not found in action 109 GetName n status Xml ModelFailure m elisiomm Output element is present us kw status XmlParser RequireChildArgument output element 15 result facies result file 17 errors status std max kw status status 119 std string text 121 kw_status XmlParser An Implementation of Markov Random Fields ME RequireChildAttribute output element 123 result format 125 text errors 127 true lowercase amp Definitions GetFileFormatNames 129 std vector lt std string gt const_iterator v un
16. e interaction must be given the format is described in Section 5 3 The neighbour hood is specified by the parameters in Section 4 5 Simple Pott s model potts The interaction is given by Eq 8 but with f 1 Any input file for interaction potential will be ignored The interactions are analogous to those of the Ising model but an Ising model not necessarily only for nearest neighbours The neighbourhood is specified by the parameters in Section 4 5 It is possible to use the interaction type general also to simulate the Pott s models as well as more general two particle interaction models just specify the input file in accordance with the desired model Sequential simulation can be carried out for two particle interaction models i e for the param eter lt type gt set to two particle generalized potts or potts The sequential conditional probabil ity p 2 mi is then assumed to have the form shown on the right hand side of Eq 7 alternatively as specified by Eq 8 Be aware that this does not imply that the probability distribution p z z_ is on the same form For sequential simulations the interaction type general has not yet been implemented If this is nevertheless attempted an error message will be written to screen and to the simulation log and the program will abort Parameters in model file Description Details Valid values lt interaction gt
17. e same thing 0 s are otherwise arbitrary in the sense that the corresponding p defined through Eq 5 will be a Markov random field for any choice of F functions provided they fullfill the clique condition In particular a function F may be zero even when the corresponding cells provide a clique 2 3 Implemented MRF models Our implementation is done under the assumption of stationarity That is we assume the condi tional probability distribution p z z 7 0 to be independent of the actual location of cell i only depending on the facies configuration of cell i s neighbouring cells Under this assumption we have implemented 1 General MRF The form of the conditional probability function is given by plalz j x wn Fa DP Paa 6 i I k Notice that the external field F is independent of the position i and only the relative po sitions of the cells are used to label the interaction parameters This is a consequence of the stationarity assumption 8 ME An Implementation of Markov Random Fields 2 General two particle interactions In this case all interactions of order higher than 2 are as sumed to be zero The conditional probability distribution thus has the general form p zi z j Oj x exp Fas Fila 20 7 l 3 Pott s model simple and generalized This is a special form of the two particle interaction model where the interactions of Eq 7 are assumed to be of the form Filzi Zi
18. en box This means the sequential neighbourhood depends on the path Normally the pa rameter lt max gt should be set to a number larger than the number of cells in the grid in which case all visited cells that are within the box are included in the sequential neigh bourhood Another extreme is to use a box larger than the grid and define the sequential neighbourhood to consist of all the last lt max gt visited cells no matter where in the grid they are If neighbourhood lt type gt is sphere Sperical neighbourhoods are supported for the MCMC simulation but not for sequential simulation an error message and program abortion results if this is attempted Radius of sphere lt radius gt The parameters of this section are neglected by the program if the interaction method is set to general In that case the input potential itself implicitly defines the neighbourhood Parameters in model file Description Details Valid values lt neighbourhood gt lt type gt box lt type gt lt type gt neighbourhood type lt type gt box sphere lt nx gt 10 lt nx gt lt nx gt lt shiftx gt box size and lt nx gt lt shiftx gt lt ny gt 10 lt ny gt shift of center x direction 0 1 2 lt nz gt 0 lt nz gt lt shiftx gt 0 lt shiftx gt lt shifty gt 0 lt shifty gt lt shiftz gt 0 lt shiftz gt lt radius gt 1 lt radius gt lt radius gt radius of spherical lt radius gt 0 1 2 nei
19. er lt type gt The param eter can take two values box or sphere For each of these there are additional parameters that need to be specified They are described in the following If neighbourhood lt type gt is box Box neighbourhoods are supported both for the MCMC simulation and for sequential simulation Periodic boundary conditioned are used if the neighbourhood type is box 1 The second method requires a new input potential file for each change of neighbourhood Hence if the main purpose of the simulations is to study the effect of changing the neighbourhood for an otherwise fixed potential the first method may imply less work on producing input files An Implementation of Markov Random Fields NE 13 Box size in x y and z direction lt nx gt lt ny gt lt nz gt Shift of center in x y and z direction lt shiftx gt lt shifty gt lt shiftz gt If shifts are set to 0 then the box is centered at the cell in question with lt nx gt cells in both positive and negative x direction and similarly for lt ny gt and lt nz gt For sequential simulations there is also a parameter lt max gt that allows the user to limit the number of last visited cells to be checked for whether they are part of the given box or not That is of all previously visited cells each of the last lt max gt cells are con sidered being part of the sequential neighbourhood provided they are inside the giv
20. er 2006 Publication number SAND 06 06 Abstract The main purpose of these notes is to provide documentation of version 0 1 of the Markov random field implementation done for the Multipoint project The report is mainly meant for internal use in the project serving as a user manual for running the implemented program and as a guide for programmers wanting to modify it Keywords Markov random field implementation Target group Multipoint project participants Availability Internal Project Multipoint Project number 808002 Research field Reservoir characterization Number of pages 30 Copyright Norwegian Computing Center Contents 1 Introduction 2 MRF the implemented models 2 1 Markov random field definition eG te oe 2 2 2 General expression based on the Hammersley Cilifford 2 3 Implemented MRF models 3 Program structure 3 1 Environment 3 2 Overall structure 3 3 Actions 4 input parameters LL 2 4 1 Overall simulation specification 42 Path 4 3 Annealing 44 Grid size 4 5 Neighbourhood 4 6 Interaction type 4 7 Lithology 4 8 Output 4 9 File directories 4 10 Miscellaneous 5 Format of input files 5 1 General interactions 5 1 1 File format en 5 1 2 Internal labeling of cells in maximal clique 5 1 3 Examples of potential file 5 2 Two particle interactions 5 3 Generalized Pott s interactions 5 4 Closing remarks on interaction models 6 Conclusion References A Example of model file
21. ghbourhood lt max gt 3001 lt max gt lt max gt the number of cells to be lt max gt 0 1 2 checked for being inside box if lt neighbourhood gt sequential simulation 4 6 Interaction type The interactions can for MCMC simulations be of four types 14 General interactions general In this case there is no restriction on the joint probability distribution other than that given by Eq 6 An input file that specifies the potential must be given The file specifies all interactions and implicitely also the neighbourhood The file format is described in Section 5 1 For this interaction type the neighbourhood parameters of Section 4 5 are ignored even when if they are explicitely listed in the model file ME An Implementation of Markov Random Fields General two particle interactions two particle In this case the joint probability distribution is restricted to have the form given by Eq 7 An input file that specifies the potential must be given The file format is described in Section 5 2 All interactions listed in the input file will be used provided that the interacting particles are within each other neighbourhood with the neighbourhood being defined by the parameters of Section 4 5 Generalized Pott s model generalized potts This is a special form of the two particle in teraction model and the interaction is given by Eq 8 An input file that specifies the spatial part f of th
22. ined in globalcommands h and globalcommands cpp The first template argument given to the NRlib UI TaskManager initializer function defines the global commands Subsequent ac tions can be declared by adding new template arguments in this function An action class e g FaciesSimulation should inherit from NRlib UI Action and should contain a copy constructor a Parse and an Execute function 4 Input parameters The input parameters for simulations are described in the following grouping parameters that naturally belong together 4 1 Overall simulation specification This section describes the parameters that determine which simulation method to use the num ber of iterations the initial grid configuration and the number of independent realizations to be made lt method gt Two types of simulation methods have been implemented Markov Chain Monte Carlo simulations MCMC and Sequential Simulations Markov Chain Monte Carlo simulations The main implementation uses standard MCMC simulations using a single site Gibb s sampler The simulation starts with some arbitrary grid configuration then updates each cell a number of times Each update is done using the conditional probability p z za Theoretically as the number of updates increases the probability distribution for the grid configurations converges towards the joint probability p 21 22 heey ZN Sequential simulations Sequential simulations are to some extent supporte
23. lso presents the probability distributions that are supported by our implementation 2 1 Markov random field definition Consider a finite regular grid in arbitrary many spatial dimensions and let the one dimensional index i label the cells of the grid The set of all cells is S 1 2 N In order to define a Markov random field it is necessary to introduce the concept of a neighbourhood system Let 0 be a collection of cells where the label i points to one specific cell By definition the set 0 0 i S is a neighbourhood system if i i not neighbour to oneself 1 and j i amp i O symmetry 2 The set of cells 0 is then the neighbourhood of cell i A Markov random field on the grid is specified with respect to the neighbourhood system By definition a probability distribution p 21 z2 zy that satisfies p zi z p z zj j Oi 3 is a Markov random field provided that it also fulfills the positivity assumption Vi p 2i gt 0 gt p z1 22 2n gt 0 4 The positivity assumption says that all combinations of single cell configurations are valid Equa tion 3 expresses that the conditional distribution p z z_ depends only on the cells in the neigh bourhood of cell i In this sense the interaction between the cells is limited by the neighbourhood system There exists a wide range of different neighbourhood systems Some common examples im plicitely assuming the metric
24. lt topdir gt Both folders are created if they do not already exist Parameters in model file Valid values lt topdir gt topdir lt topdir gt Any folder name lt prefix gt sim name lt prefix gt Any folder name 4 10 Miscellaneous Seed for the random number generator may be provided In addition there are two parameters that specify the level of detail with which simulation log file information and information to screen is given Log file information is written to the file mrf log lt seed gt The seed for the programs random number generator If not specified seed is com puted based on the computer s clock lt logfile gt lt level gt Information level to be used for writing to file 4 Note that the Mersenne twister implementation gives different random number sequences on 32 bit and 64 bit plat forms An Implementation of Markov Random Fields ME 17 lt screen gt lt level gt Information level to be used for writing to screen Parameters in model file Valid values lt seed gt 345 lt seed gt Any positive integer or noth ing lt messages gt Valid information levels lt logfile gt 0 no messages given lt level gt 5 lt level gt 1 severe error messages lt name gt logs mrf log lt name gt 2 warnings lt logfile gt 3 information messages lt screen gt 4 detailed information lt level gt 5 lt level gt 5 debugging messages lt screen gt 6 requires DEBUG compile lt me
25. ml ModelFailure else Get all children const Xml Elements facies element lithology element gt FirstChildElement facies for facies_element facies_element facies_element gt NextSiblingElement facies std string name double fraction 0 kw status XmlParser GetChildArgument facies element name name errors facies push_back name status std max status kw_status kw_status XmilParser GetChildArgument facies element Htracthionwe fraction errors target facies volume fraction push back fraction status std max status kw status Sometimes a program parameter has a predefined list of allowed values In the program we may want to represent these with enums In the input file for the sake of user friendlieness we must use strings Hence we have a list of allowed input string values 28 p Read method stdi string text kw_status XmlParser GetChildArgument element method text read_string errors Default GetMethodName default true Lowercase flag An Implementation of Markov Random Fields 5 amp Definitions GetSimulationModes 87 status std max status kw status 89 Now we have an allowed string Loop over the allowed values 91 and find the corresponding enum value represented by an unsigned integer 93 std vector lt std string gt const_iterator v 95 unsigned int mod
26. nt siden es model a model file name TaskManager model model Initialize lt MRF GlobalCommands MRF FaciesSimulation gt 0 model Execute std cout lt lt timer Report lt lt std endl std cerr lt lt Screen Boldface Exception caught lt lt exception what lt lt std endl return EXIT_SUCCESS return EXIT_FAILURE Figure 1 Flowchart of main in MRF cpp An Implementation of Markov Random Fields NE S n generated by the progam if set to constant all cells have z 0 in the initial configuration For sequential simulations the simulation starts with an empty grid and the parameter has no effect lt realizations gt The desired number of independent realizations should be specified Each independent realization is achieved by starting from an independent initial configuration and then run the specified number of iterations for each cell on the grid For sequential simulation the number of iterations is 1 Parameters in model file Description Details Valid values lt method gt mcmc lt method gt Specifies whether MCMC or se memc or quential simulation is to be used sequential lt iterations gt The number of iterations per cell 1 2 3 1000 Will be reset to 1 if sequential sim lt iterations gt ulation is chosen lt initialize gt lt type gt Specifies the initial grid random lt type gt random lt type gt
27. of this file Be aware that this example only illustrates the file format and is not guaranteed to represent a stationary MRF model 6 02 9 0 5 6 1 5 0 1 1 0 6 1 5 3 0 6 1 2 6 0 1 8 0 4 9 2 0 1 Figure 4 Clique configurations corresponding to the more general file example 5 2 Two particle interactions This format of the input file is to be used whenever the parameter specifying the interaction type is set to two particle see Section 4 6 At the very beginning of the file the integers L Ly Lz K must be given These are used by the program to identify the index of the correct file element for a given relative position Ar j of the two interacting particles The interpretation of these parameters are the same as for interaction type general see Section 5 1 but the actual neighbourhood used to specify the MRF model may be different see Section 4 5 An Implementation of Markov Random Fields ME 21 For general two particle interactions both spatial dependency and dependency on facies must be listed Recall from Eq 7 that the general form of the probability distribution is assumed to be pleier j i x exp Fe z 5 Filzi 241 12 l The two particle interactions F j i zj are assumed to be listed for a number of relative positions Arij Az j Ayij Azi that constitute a box That is we assume DE EEE Lr 1 1 0 1 Lo
28. r Ly Ly Lz K 10 The interpretation of these integers is that the size of the neighbourhood is 2L 1 2Ly 1 2L 1 with a maximum clique size Lz 1 Ly 1 L 1 The number of facies is K The four integers are used by the program to interprete the remaining lines of the input file in a consistent way All remaining lines of the input file should have the following format 18 ME An Implementation of Markov Random Fields Inew n N fai D 0 where N is the number of cells in the considered clique fa gives 1D indices that determine the locations of the clique cells relative to an internal labelling convention of the maximum clique fi is the facies indicators of the cells with each f 0 1 K 1 0 is the parameter that gives this clique configuration s contribution to the potential see Eq 9 The two integers Inew and 7 provide extra information that allows the program to be more effi cient Their interpretation is 1 if this row is the first row in a sequence of one or more rows Inew 0 otherwise n gives the number of subsequent rows with the same set ai N 1 There is no harm done by setting Lew 1 and 1 in all lines but the simulation will be somewhat slower 5 1 2 Internal labeling of cells in maximal clique The maximal clique has dimensions Lz 1 Ly 1 1 Figure 2 illustrates the labelling of cells i e the possible values in the set a _ for a 2D case The e
29. signed int mode 0 s for v Definitions GetFileFormatNames begin v Definitions GetFileFormatNames end v mode 133 ie oy SS Gera break 135 137 facies_result_format_ static_cast lt Definitions FaciesFormat gt mode status std max kw_status status 139 At the end of the parsing function we check the status and print out error or warning messages as required Depending on the status we return true for success and false for failure i Check final parsing status and write out any errors if status gt Xml ModelWillDo w log Message Logger SEVERE lt lt errors return false failure u E else if status gt Xml PerfectModel 15 log Message Logger WARNNG lt lt errors 147 return true success 30 ME An Implementation of Markov Random Fields
30. ssages gt flag 5 Format of input files This section describes the format in which the potential functions should be listed in order for the program to interprete them correctly The potential file lists the potential values in a specific order and it can for instance be made by using MatLab Whether the listed values result from evaluating some known mathematical function or from estimating parameters from a training image is up to the user 5 1 General interactions This format of the input file is to be used whenever the parameter specifying the interaction type is set to general see Section 4 6 In this case the interactions are only restricted by the Hammersley Clifford theorem A con venient expression for the conditional probability is then p 2ilze x expt Y ac 9 C YEC zc The outer sum is over all cliques such that cell i is in the clique the inner sum is over all config urations zc of the clique C The parameter 0 gives the strength with which this configuration contributes to the overall potential If 0 0 the configuration has no effect on the potential Eq 9 gives the key to interpreting the potential file that is described in the following Only pa rameters 0z 0 will be listed in the file for the obvious reason of saving memory 5 1 1 File format The potential file should start with 4 integers that give information about the neighbourhood and number of facies The integers are in orde
31. st cases is specified through its own input file see Sections 4 6 and 5 for more details on this It is not necessary to recompile the program in order to run a different interaction model or use different parameters Simply provide a new model file and or a new file for the input potential Appendix A gives an example of a model file and appendix B explaines how the program parses the file 3 2 Overall structure The main function in MRF is found in the MRF cpp Classes from NRlib are used to run stop watch for timing purposes hold program data such as version number copyright holder and so on read and keep command line options initialize the task manager with tailor made actions 1 2 3 4 read model file name from command line 5 6 parse the model file 7 execute the model AD An Implementation of Markov Random Fields NR 8 catch exceptions and report error messages 9 report time usage and finish A flowchart for this function is shown in Fig 1 3 3 Actions The implementation is designed so that it it possible to define different actions to be taken by the program So far the only action that is fully implemented is simulation One more action allowing for the calculation of likelihood and different probability functions for a given input image is presently being developed but will not be described here Global parameters valid for all actions are read and kept by GlobalCommands declared and def
32. t level gt lt screen gt lt messages gt lt topdir gt lt topdir gt lt prefix gt real lt prefix gt lt seed gt 345 lt seed gt lt action type simulate facies gt lt method gt mcmc lt method gt lt annealing gt lt max temperature gt 10 lt max temperature gt lt temperature reduction factor gt 0 9 lt temperature reduction factor gt lt steps gt 10 lt steps gt lt annealing gt lt fraction control gt 1 lt fraction control gt lt iterations gt 1000 lt iterations gt lt realizations gt 1 lt realizations gt lt grid gt lt nx gt 100 lt nx gt lt ny gt 100 lt ny gt lt nz gt 1 lt nz gt lt grid gt pp 24 INRES An Implementation of Markov Random Fields lt lithology gt lt facies gt lt name gt sand lt name gt lt fraction gt 0 55 lt fraction gt lt facies gt lt facies gt lt fraction gt 0 45 lt fraction gt lt name gt shale lt name gt lt facies gt lt lithology gt lt initialize gt lt type gt random lt type gt lt initialize gt lt path gt random lt path gt lt interaction gt lt type gt general lt type gt lt beta gt 0 5 lt beta gt lt input file gt dSeq_10_10 dat lt input file gt lt interaction gt lt neighbourhood gt lt active only if interaction type is not general gt lt type gt box lt type gt lt nx gt 10 lt nx gt lt ny gt 10 lt ny gt lt nz gt 0 lt nz gt lt shiftx gt 0 lt shiftx gt lt shift
33. t they correspond to the cell s facies indicator as z 0 1 2 respective to the order of the facies listing The program is so far used and tested only for binary facies values i e two listed names but is made so as to accept up to 256 of facies types For each facies the user of the program can specify a target value for the volume fraction of this facies If using the fraction control option explained shortly the program will then automatically adjust its parameters so that the measured volume fractions converge towards the target values Thus for each facies the following parameters should be specified lt name gt Facies name lt fraction gt Target value of volume fraction to be used in fraction control see next paragraph The volume fraction of each facies can be controlled throughout the simulation The algorithm used is equivalent to adjusting the external field of the model the adjustment being based on the target value per facies in combination with measurements done throughout the simulations The algorithm is adopted from another project at NR Sand where its utility has been tested It is described in the following lt fraction control gt the change parameter S If we define a as a S 1 f fo where f is measured volume fraction and fo is target value the conditional probability for this facies is changed with the factor e42 The value S 0 means fraction control is not performed Recommended value
34. to be Euclidean are An Implementation of Markov Random Fields NE 7 Ist order neighbourhood only grid cells that are nearest neighbours to cell i belong to 0 2nd order neighbourhood only grid cells that are nearest or next nearest neighbours to cell i belong to Oi the whole grid 0 S i The last example shows that any probability distribution that fullfills the positivity assumption is a Markov random field The typical application of Markov random field theory assumes that each neighbourhood is relatively small compared to the grid size 2 2 General expression based on the Hammersley Clifford There exists a profound theorem that gives a general expression for the form of a Markov random field the Hammersley Clifford theorem The theorem states that a Markov random field can be expressed as rosa 0 om 2 F z 5 i i Pag isz 1 lt i lt N 1 lt i lt j lt N 5 TE 5 5 gery zero Any set of cells where all set members are neighbours is called a clique Hence the Hammersley Clifford theorem expresses the Markov random field as a series expansion over 1 cliques 2 cliques up to and including the largest clique provided by the given neighbourhood system This shows that any Markov random field is a Gibbs field as the Gibbs field is defined in terms of a clique expansion of this kind The contrary is also true any Gibbs field is a Markov random field That is Markov random fields and Gibbs fields are th
35. us kw status kw_status XmlParser GetChildArgument annealing element temperature reduction factor annealing_epsilon_ errors Default GetAnnealingEpsilon status std max status kw status kw_status XmlParser GetChildArgument annealing element steps annealing_steps_ errors Default GetAnnealingSteps status std max status kw_status else using default values from constructor For all of the annealing parameters there are default values defined in default h These values are used in the constructor and hence if no lt annealing gt element is found default values are used Likewise if for example lt steps gt is missing its default is used Sometimes we need to read a group of data with variable number of input elements lt lithology gt lt facies gt lt name gt sand lt name gt lt fraction gt 0 55 lt fraction gt lt facies gt lt facies gt lt fraction gt 0 45 lt fraction gt lt name gt shale lt name gt lt facies gt lt lithology gt Note that the order of subelements for name and fraction is arbitrary For nested child elements An Implementation of Markov Random Fields 27 const Xml Element lithology_element element FirstChildElement lithology lithology_element errors std string ERROR lithology not found in action GetName n status X
36. xtention to 3D is similar The 1D index a is related to the 3D lags lz ly l by a ly ly Le 1 Le 1 Ly 1 11 where l 0 1 Lx and similarly for l and 1 5 1 3 Examples of potential file In the following we give an example of the potential file on the format corresponding to general interactions The example represents the 2D nearest neighbour Ising model with an interaction strength 0 0 5 The file should then look like this 1102 12201000 5 02201110 5 12202000 5 02202110 5 Figure 3 illustrates the clique configurations corresponding to this file The figure shows from left to right lines 2 5 of the file The numbers indicate the facies of the cells and the cells with no number do not belong to the clique considered in each case The interaction parameter 0 is in each case indicated above the figure A slightly more complex example of potential file is given by An Implementation of Markov Random Fields NE 19 lt lag I Figure 2 Convention used for indexing cells in max clique 2D illustration 8 0 5 8 0 5 8 0 5 8 0 5 Oo 0 1 Figure 3 Clique configurations corresponding to the file example for the Ising model 20 ME An Implementation of Markov Random Fields 2102 14201000 2 04201141 0 5 04201011 5 04201101 5 124013411 1 0 3 024013400 11 2 1360123451010100 1 036012345 101010 4 03601234510110 1 2 0 Figure 4 illustrates the clique configurations
37. y gt 0 lt shifty gt lt shiftz gt 0 lt shiftz gt lt radius gt 1 lt radius gt lt max gt 3001 lt max gt lt neighbourhood gt lt output gt lt measurements gt measurementlog txt lt measurements gt lt result format sgems gt result dat lt result gt lt temp format sgems gt tmpresult dat lt temp gt lt gridfreq gt 50 lt gridfreq gt lt writefreq gt 10 lt writefreq gt lt output gt lt action gt lt mrf gt An Implementation of Markov Random Fields 25 B Parsing Xml files with NRlib The MRF package takes user input from an ascii file using the Extensible Markup Language XML An xm file has a strict tree structure consisting of nested elements some of which have at tributes and content An element typically consists of two tags a start tag and an end tag pos sibly surrounding text and other elements The start tag consists of a name surrounded by angle brackets like lt step gt the end tag consists of the same name surrounded by angle brackets but with a forward slash preceding the name like lt step gt For output to screen and log file the MRF model file syntax is as follows lt messages gt lt logfile gt lt level gt 5 lt level gt lt name gt logs mrf log lt name gt lt logfile gt lt screen gt lt level gt 5 lt level gt lt screen gt lt messages gt The parsing function takes a reference to its appropriate action element as input It could be regarded as a
Download Pdf Manuals
Related Search