- Email: [email protected]
Procedures for creating a benchmarking data set
Computers & Geosciences Vol. 14. No. 2. pp. 139-150, 1988 Printed in Great Britain. All rights reserved
0098-3004/88 $3.00 + 0.00 Copyright ~ 1988 Pergamon Prt~ pie
PROCEDURES FOR CREATING A BENCHMARKING DATA SET MARK A. HERKOMM'ER PETROSPEC Computer Geology and Geophysics. 4005 Burning Tree Lane. Garland. TX 75042. U.S.A. (Received 6 May 1987) Abstract--The testing of mapping (and other spatially oriented) software usually requires a test-data set. When benchmarking a program, it is worthwhile to use a data set that simulates the various types of data distributions observed in the real world: random (random wells), cluster (oil fields), grid (mineral surveys). or linear (seismic surveys). A FORTRAN V program, BENCHMARK, is presented that samples a 2-D speatial function with these data distributions. Key Words: Data distribution. Software testing, Spatial analysis.
INTRODUCTION When tcsting ncw software, it may be useful to have a benchmarking data set with known characteristics. The program presented in this paper, BENCHMARK, is designcd to produce data sets of the form Y, Y, Z. Data sets of this type are used for gridding and contouring, 2-D fihcring, and othcr spatially oricntcd applications. BENCHMARK was crcated originally to assist in testing of the accuracy of computer-contouring programs. Computer-mapping programs employ many different algorithms for handling the spatial input data. Many of these are designed for special data distributions. Broadly speaking, there are three basic data distributions: random, clustered, and regular (Davis, 1973). When benchmarking software, it is worthwhile to use data sets that reflect these types of data distribution, Each of the aforementioned data distributions has a real world analog in mapping applications. The random distribution is similar to random well locations, as seen in infantile or mature petroleum basins. The clustered distribution in the oil-well field groups are observed in juvenile to middle-age petroleum provinces. The regular grid sampling method is used in mining operations. A special (fourth) situation of the clustered distribution is data clustered along lines. Seismic survey data are a well-known example of this linear type of data distribution.
THE BENCHMARK PROGRAM B E N C H M A R K is a F O R T R A N V program (Appendix 1) designed to simulate the four types of data distributions sampling from a known 2-D surface. The input to the program is an output filename to receive the data, the range to sample the X and Y data from type of sampling mode (random, line, clustered, or regular), and the number of sample points. The Z values are determined by evaluating a single. continuous 2-D function at each of the selected X 139
and Y locations in the "map" area. For the function defined with B E N C H M A R K , the map area is 2 units by 2 units square. The minimum X a n d Yvalues used are 0.0, making the maximum X and Y values 2.0. For the random well-!ocation simulation, X and Y values were located in the range of 0.0-1.0 using RAN, a pseudorandom number generator (Press and others, 1986). They are scaled to the preselected size of the map area. The X and Y values then are put into the function Z F U N C and a Z value calculated. The seismic survey shotpoint simulation data set is better described as regularly spaced control points gathered along random lines. The initial and terminal points (X, Y pairs) of the line are random numbers. The length of the line is calculated and if it is determined that there would be < 5 or > 500 shot points on the line, it is rejected and new initial and terminal points are selected. After a line is accepted, the initial point of line is written to the data set as shot point No. 1. Subsequent points are offset from the previous point by a user specified distance along the line. Points are added along the line until it is populated fully, at which time a new line is selected. The oil-field grouping simulation data set have the characteristic of being clustered. The clusters are synthesized by first selecting a center for the cluster at random. The center of the cluster is cluster point No. I. Then the balance of the control points are determined at random having a deviation of no more than the search distance in X and Y cluster's center. The general shape of a cluster is square. The mineral-survey grid simulation is a rectangular grid of X and Y locations. The grid is created from dividing the Xand Y ranges into uniform sample intervals as determined by the number of samples along the X and Y axes. Samples are taken from the surface based on the X and Y locations. Output from B E N C H M A R K is a file containing the sampled surface. The record format of each distribution is similar so different distributions can be mixed. The general output record format is:
140
M.A.
HERKOMMER
Z.O,
x 1.9 1.6 x
2.6 X 1.9 1.6X X
1.8
1.0 -0.6 0.0 x x x-0.8
1.6
~.,~ 1.4
1.2
-2-2x
-o.fV
-1.5x
-xl.9
x 0.0 x
-1.4 x-1.3
-2.3 -2.5 x
O.9 X 0.4
~.9
0.Tx
x
4.6 X
x 0.5 2.2 x
1.7 x x
5.9 x 5.9 x
~7 4.0 x
-O.8 x
2.0 x
4.2 x
-3.3 x
0.8
0.0
-2.1x -4.4 Ix 0,2
-2.2 1 0.4
I~ 0,6
2.0
1.3x
x2.0 0.8 x
-2.0 x
x
2.2 x x 1.3
-3.5
x-1,3
3.5 x 4.8 x
x-5.2
x-3-2 x 0.2
xi.4 -O.2 x
-3.1 x
-2.9 x
6.3 x
2.4 x 2.3 x xl.8 2.2 x 2.0 1.9 x # 1.~ x
1.8 x
0.9 x
4.t
6.8 X
-0.1 xX 0.3
2,2 x
0.8
0.4
0.2x
x 1.1
1.8
4.6
1.4 1.Z x
1.4
2.0
0.6
2.1 X
x1.1
0.1 x 1.6~
x
2.2 2.1 X X 1.8 X
0.8 x
1.7 x 1.0
2.0 X
1.6 x
x -5.7
-6.9 x
-6.3 x
4,9 x 3.~x
-Z6 x
-5.8 x
x -7-3 I
1
l
1
0.8
1.0
1.2
x -3"91
1.4
1.6
I 1.8
I 2.0
Figure I. Random distribution--simulated well locations. pointtlsymblltyllclass]lindiv[l--x value--II--)' valuc--lJ--z value value--l, ] where point is the sample point number, symb is the symbol code, ty is the type, class is the group number of the particular type, indiv is the individual number withing the particular group, x value is the X location, ). value is the Y location, and z value is the sampled surface value at the determined location. Appendix Z shows an output file from B E N C H M A R K that contains all four distribution types. The random distribution record format is (2i6, 3x, 'P', i7.7, 7x, 3glS.7e2). The data values it contains are: point number, symbol code 'P' type, point number (zero padded), X, Y,.and Z. Figure ! shows a map with 100 random locations posted. The linear distribution has the format (2i6, 3x, 'L', i7.7, i7, 3g15.7e2). These data values are: point number, symbol code, 'L' type, line number (zero padded), "shot" point number on the line, X, Y, and Z. Figure 2 shows a map with 100 simulated seismic line data with a sample rate of 0.10 units along each line. The cluster distribution record format is (2i6, 3x, 'C', i7.7. i7, 3g15.7e2). The data values on its record are: point number, symbol code, 'C' type, cluster number (zero padded), cluster point number in cluster, X, Y, and Z. Figure 3 shows a map with 100 clustered data points with each cluster having 20 points and a search distance of 0.15 units. The regular grid distribution has the format: (2i6, 3x, 'G', 2i7, 3g15.7e2). These data values are: point number, symbol code. "G' type. X grid node index, Y
grid node index, X, Y, and Z. Figure 4 shows a map with a grid of 10 x 10 grid nodes. T I l E A R B I T R A R Y S U R F A C E U S E D BY BENCIlMARK
The arbitrary surface sampled by BENCHMARK is a single, continuous 2-D function. It is the sum of two discrete functions; an 8th-order polynomial and interference wave forms intersecting at an angle of 90.0- (Appendix i). Users of BENCHMARK can substitute their own function if they wish. Figure 5 shows the surface contoured on a grid of 101 x 101 grid nodes. The 8th-order polynomial part of the function was constructed in several stages. First, a set of 45 random X, Y, and Z coordinates was generated. The Xand Y values of the random numbers were in the range of - 1.0-3.0, and the range of the Z values were in the range of 0.0-2.0. The random numbers then were fitted with an 8th-order polynomial function by the method of least squares, the coefficients of which are shown in the program listing (Appendix I, subroutine ZFUNC). Because the test map area considers X and Y values sampled only in the range of 0.0-2.0, there are only minor "'edge effects" in the surface. A simple interference wave form function was added to the 8th-order polynomial surface to simulate noise and lower the degree of autocorrelation within
Creating a benchmarking
d a t a set
141
2.0
1.8
1.6
1.4
1.2
1.0
1.0 Xl. 0 x -2.1
2.2 x 1.0 x
2.4 x
1.9 x
x2.0
1.6 x 1.3 x -2.3 x 1.2 x 1.6 -2.1 0.9 x x x..1. 5 0.6x 1.9 x x x 2.2 _ 0 . -9-0 . ; ¢ 0.5 x 2.4 x 0.2 x 0.1 ~0.7 x2.5 x 1.1 xO.9 x 1.7 1.5 1.8 2.4x ~1.4 x x 2.3 .62 xl.O x 6~ ~(8 1.8x 2.1 X3. x
4.0 ~
0.8
1.6 ~
4.9 x 5.3 x
0.6
-0.2 X -1.3 x ~Om4 -2.3 x x 1.3 -3.1 x x - 3x. 6
0.2 -0.3 I~ 0.2
3.1 x -3.7 1.5 x x_3.9 ~5.7. -1.8 x 0.0 x xO.O 0.8 x -1.4 x lzr I l 0.6 0.8 0.4
1.4 x 0.5 1.6 x x 1.9 0.9 x x 1.4 2.0 x
2.41"~ _,11"" 2., 2.~2~.
1.1 x 0.4 x
5.1 x 4.3 x
0.4
0.0
2.4 x
Z.2 x
I 1.0
x 2.02.0 1.6 X x xl.5 1.8 0.8 x X 0.6 x 1.6 x - 0 . 6 x x- 0 . 2 -1.9 -2-0 x x -3.5 -2.3 x x -5.0 -3.2 x x -6.4 x -3.8 x ! I [ 1.2 1.4 1.6
F i g u r e 2. L i n e a r d i s t r i b u t i o n - - s i m u l a t e d
I 1.0
I 2.0
seismic survey.
2.0
1.0
1.6
0.4 x -0.8 -1.8d~1.7 X x - 0 . 9
~, x-2"~-2.,
-oi2
-2.0 1.4
1.3 x-O.7 x x-0.2 -0.4 x 0.5
-1.5 x -I.1
1.2
1 . ~ "1 2.3
1.0
1.2 x 2.3 x x 1k4 xl.6 2.1 x 5.4 :3.7 1.3 1.0 x 2.0 2.0 5.5#x5.6 x x x x x x 1.2 Xl.9 x x3.7 - 0 . 4 - 0 . 9 0.3 12~k 1.6 x5-6 3.9 -0.5~F x 0.SXoX91.3 x 6x8 5x'O x -I 2 x -2.3 -27 1.4 2.3 1.1 _.. -2.9 x ;(" x x 1.9 x x , . 5 " ~ X x ..x32 1 4 2"x32.6 x2.9 ; ~ 0.6 -3.4 ' " g. ~¢ X3.4 • x-O.9 X , = -0.5 16~2.4"29 ~z'0"3-"'~ - ~ . 4 "2.3 xx x 3.7 -1.4 -4.3X~4"? x x x 3.4 ~35Xx -4.7 -1.5 1.1 XxX4.2 • x - 4 . 6 -7.4 4~0
0.8
0.6
0.4
O.Z
-4,8
0.0
I 0.2
1 0.4
I 0.6
-4.7. [ 0.8
I 1.0
[ 1.2
F i g u r e 3. C l u s t e r d i s t r i b u t i o n - - s i m u l a t e d
I 1.4
oil fields.
[ 1.6
I 1.8
I 2.0
M. A. I-IERKOMM]LR
142
2.0x5.1
2.9
x
3.1
3,0
2.8
2.6
2.5
1.9
1.2
2.9
1.0 2.5
0.5
1.1 X
1.7 X
2.0 X
2.3 X
2.4 X
1,9 X
1.3 X
2.0
1.6 I-1.9 X
-2.3 X
-1.0 X
0.1 X
0.9 X
1.5 X
2.0 X
1.'/' X
1.1 X
1.4
1,4 ~ 3 . 8 X
-2.6 X
--0.9 X
0.1 X
0.8 X
1.5 X
2.2 X
2.1 X
1.3 X
0.2
1.2 ~ 1 . 8 x
0.1 x
1.4 X
1.5 X
1.4 X
1.7 X
2.4 X
2.5 X
1.7 X
-0.3
1"0!.2
4.4x
4.7 X
3.1 X
1.6 X
1.0 X
1.6 X
2.:) X
2.0 X
-0.4
6.4 X
3.2 X
0.0 X
-1,3 X
-0.6 X
0.7 X
1.9 X
2.1
5.0 x
1,4 x
-2.9 x
-4,8 x
-3.8 x
-1,2 x
1.8 x
4.9
2.9 X
-1.4 X
-6.1 X
-0.1 X
-6.4 X
-2.6 X
2.2 X
7.9
l°i ° o4
r2i ' 0.6
I-7i4 o.0
%31 ,6
~o ,o
98~ 20
1.8
2.0
x
x
v
K
X
X
X
X
X
t
0.8 t . 6 x
7.2 x
0.6 1 . 9
6.5
0,41
x
0.1
1.6
0.2i
0.0
x
x
x
x
x
X
x
x
x
I - 24.0 -5.6 x ix 0.2
j - 9 i, ~ - ~ 8 , ,o ,2 ,4
Figure 4. Grid distribution--simulated mineral survey.
2.0 1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.2
0.4
0.6
0.8
1.0
L2
1.4
1.6
Figure 5. Benchmark surface used to sample distributions.
Creating a benchmarking data set the surface (Harbaugh, Doveton, and Davis, 1977). Adding the wave function makes it more difficult to model the surface. Only those data sets using at least 10,000 control points would be able the "see" this overlain interference surface. To create the "noise", two wave patterns are intersected in phase at 90.0 ° with a wavelength of 0.04 units and an amplitude of 0.05. This is seen as the last line in the equation for the function in Appendix I.
143
F O R T R A N may have difficulty with the number of continuation lines in ZFUNC. This can be overcome by breaking the function into pieces and summing the parts. The program requires approximately 64 K bytes of memory to execute. Although run times may vary, on an IBM PC/AT, 10,000 points can be generated in about 4 min. REFERENCES
COMPUTER REQUIREMENTS FOR USING BENCHMARK B E N C H M A R K was written originally on a VAX 11/780 in ANSI standard F O R T R A N V with minor VAX extensions. The version presented here has been ported to an IBM PC/AT and compiled with R y a n McFarland FORTRAN. Some other versions of
Davis, J. C.. 1973. Statistics and data analysis in geology: John Wiley & Sons, New York. 550p. Harbaugh. J. W., Doveton, J. H., Davis, J. C.. 1977, Probability methods in oil exploration: John Wiley & Sons, New York. 267p. Press, W. H., Flannery, B. P.. Teukolsky. S. A., and Vetteding, W. T., 1986. Numerical recipes: Cambridge Univ. Press, Cambridge, 818p.
APPENDIX 1
Program BENCHMARK c cnnmmnlaiun~nnl|ummnilamammllimmlamnmnmmmaeallmul|mmN|iim|anaugmmmnulnn
c c c c c c c
program to select data points randomly from a given surface data selection is in 4 modes; random points and points along random lines, points in random clusters, and points on a grid. written
by Mark A. Herkommer
c
modified
for pc on 4-28-87
on 6-6-83
Cm.minmnnmmanmmnn.uumnennminmimmnmmminnmm.wmi.mmmmmmmm...mnnnm.iwenlnl~ c
character cyn*l,coutfil*40 integer*2
cbell
parameter
(cbell - 7)
c c common /range/ C
. . . . . . . . . .
"
. . . . . . .
xmin, xmax, xdif, ymln, ymax, ydlf
"
100
wrlte(*,1100) &'Welcome to BENCHMARK -- a Fortran V language program used to', &'generate data sets for benchmarklng purposes. Given a user', &'defined 3-d function, benchmark samples the surface at random', &'locatlons. The program uses 4 sampllng techniques, selectlng', &'random points (mode-0), points along random lines (mode-l),', &'points clustered around a random point (mode-2), or points', &'arranged in the pattern of a regular grid (mode-3). For', &'"polnts-along-random-lines" mode, the sampling distance along', &'llne selected must be specified. Lines will have more than', &'5 and less than 500 points. For "polnts-clustered-around', &'a-random-polnt", the cluster distance and number of points per', &'per cluster must be specified. For "points-on-a-regular-grldO',' &'an x and y increment must be selected. The range for the data', &'range to be selected is input as xmin, ymin, xmax, and ymax.' 1100 format(6x,a) wrlte(*,1100) ' ', &'The output file from benchmark is a sequential data file with', &'each record having the followlng format:'," ', &'point number (16), llne/cluster Id (aS), x (g15.7),'// &'y (g15.7), z (g15.7}'
c c ..... get parameters for data point selection c write(*,1200} 1200 format(//," enter the output filename read(*,'(a40)')coutfll open(unlt - 10,file - coutfil,status - 'unknown')
: ')
144
M . A . H£RKOM.'~R
xmin - aread('enter the x-mlnimum ymln - aread('enter the y-minimum xmax - aread('enter the x-maximum ymax - aread('enter the y-maximum xdif - xmax - xmin ydlf - ymax - ymin i£ ( x d l f . l e . O . O .or. y d l f . l e . O . O ) & stop 'minlmum/maximum r a n g e error"
=') =') =') =')
c c ..... sampling mode c 200 wrlte(*,1210) 1210 f o r m a t ( / , ' e n t e r t h e s a m p l i n g m o d e ; ' , t 3 0 , ' O - random points',/ & t30,'1 - points along random lines',/ & t30,'2 - points clustered around a random point',/ & t 3 0 , ' 3 - p o i n t s on a r e g u l a r g r i d ' , / ) c mode = aread('enter the sampling mode c if(mode.eq.O) then npnt a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eg.1) t h e n d i s t s a m - a r e a d ( ' e n t e r the s a m p l e d i s t a n c e a l o n g l i n e s npnt = a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eq.2) t h e n srad a r e a d ( ' e n t e r the s e a r c h d i s t a n c e for all c l u s t e r s n c l u p n t - a r e a d ( ' e n t e r the n u m b e r of p o i n t s p e r c l u s t e r npnt a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eq.3) t h e n nx a r e a d ( ' e n t e r n u m b e r of s a m p l e s a l o n g t h e x a x i s ny a r e a d ( ' e n t e r n u m b e r of s a m p l e s a l o n g t h e y a x i s If (nx.le.l .or. ny.le.l) & s t o p ' s e l e c t s a m p l e r a t e s g r e a t e r t h a n I' x i n c - x d i f / ( n x - I) y i n c - y d i f / ( n y - i) npnt - nx*ny c else write(*,'(/,a)') ' p l e a s e s e l e c t O, I, 2, or 3...' g o t o 200 endif isymb a r e a d ( ' e n t e r t h e s y m b o l c o d e for t h e p o i n t s c c ..... write out selection parameters c wrlte(*,1250)coutfll 1250 f o r m a t ( / / / / ' d a t a p o i n t s e l e c t i o n p a r a m e t e r s ' , / & , ............................... ,,/ &' o u t p u t f i l e n a m e : ',a40) call awritef('x-minlmum - ',xmin) call awritef('y-mlnimum - ',ymin) call awrltef('x-maximum = ',x max) call awritef('y-maximum = ',y max) if(mode.eq.O) then write(*,1260) 1260 format(/,' sampling mode - ', & 'random points') e l s e i f (mode.eq.l) t h e n write(*,1270) 1270 format(/, ' sampling mode - ', & ' p o i n t s a l o n g r a n d o m lines') call awriteE('sample distance along lines - ',distsam) ishtpnt - 0 illne - 0 c e l e e l f (mode.eq.2) t h e n write(*,1280) I ', format(/, ' sampling mode 1280 & 'points clustered around a random point') call awrltef('search d i s t a n c e for a l l c l u s t e r s = call awrltei('number of p o i n t s p e r c l u s t e r = iclupnt - 0 icluster - 0
elseif
(mode.eq.3)
then
',srad) ',nclupnt)
=')
=')
=') =')
=') =') =')
-') =')
=')
Creatingabenchmarkingdata
1290
145
~t
write(*,1290) format(/, ' s a m p l i n g m o d e " ', & . ' p o i n t s o n a r e g u l a r grid') call a w r i t e i ( ' n u m b e r of s a m p l e s a l o n g x a x i s call a w r l t e i ( ' n u m b e r of s a m p l e s a l o n g y a x i s ix - 0 iy - 0 endif call s w r l t e i ( ' n u m b e r of p o i n t s to be s a m p l e d call a w r l t e i ( ' s y m b o l c o d e for p o i n t s -
- ",nx) - ',ny)
',npnt) ',isymb)
c c . . . . . c h e c k to m a k e s u r e d a t a is ok c write(*,1300) 1300 f o r m a t ( / , ' are t h e s e v a l u e s ok ? ') read(*,'(al)')cyn If(cyn.eq.'N'.or.cyn.eq.'n') then close(10) g o t o 100 endif c c ..... sample the surface c do 300 i ~ 1 , n p n t c if (mode.eq.0) t h e n c c ..... random points c call s a m p n t ( x v a l , y v a l , z v a l ) wrlte(10,'(216,3x,''P'',i7.7,Tx,3gl5.7e2)') & l,lsymb, i , x v a l , y v a l , z v a l c e l s e i f (mode.eq.1) then c c . . . . . p o i n t s a l o n g r a n d o m lines c call s a m l l n ( x v a l , y v a l , z v a l , d i s t s a m , i l l n e , l s h t p n t ) wrlte(10,'(2i6,3x,''L'',iT.7,17,3915.7e2)') & l,lsymb,iline,ishtpnt,xval,yval,zval c e l s e i f (mode.eq.2) t h e n c c ..... points clustered around a random point c call samclu(xval,yval,zval,srad,nclupnt,lcluster,iclupnt) wrlte(10,'(216,3x,''C'',iT.7,17,3g15.7e2)')
&
i,lsymb, icluster,lclupnt,xval,yval,zval
c e l s e i f (mode.eq.3) t h e n c c . . . . . p o i n t s on a r e g u l a r g r i d c c a l l s a m g r d ( x v a l , y v a l , z v a l , x l n c , y l h c , l x , ly,nx,ny) write(10,'(216,3x,''G'',217,3glS.7e2)') & i,lsymb, lx, i y , x v a l , y v a l , z v a l endif 300 continue c wrlte(*,1400) cbell 1400 f o r m a t ( I x , a l , ' s a m p l l n g complete.',//, &' d o y o u w a n t s a m p l e a g a i n w i t h the s a m e p a r a m e t e r s ? ') read(*,'(al)') cyn If(cyn.eq.'Y'.or.cyn.eq.'y')goto 200 c close(10) wrlte(*,1500) 1500 f o r m a t ( 1 3 ( / ) , ' g o o d - b y e and t h a n k y o u for the u s a g e . ' , 1 0 ( / ) ) end c C l l l ] l l l l l l l l l l l l l l l l l l l l l l l l l a U l l l l l i l l l l l l i i l l l l l i i l l l l l l l l l l l l l l l l l l
C function aread(c) c h a r a c t e r c*(*) wrlte(*,1000)c 1000 f o r m a t ( I x , a , ' ') read(*,'(f20.0)')aread return
end
146
M . A . HERKOMMER c C
i
l
l
l
| a a n e m n
N i l
n m e | g a a u m
i i
e | a m m
l
nlnllemm
i
i
l
i
mm
l
i
l
l
i
mnlmmn
i i
| m
n m n n m n
| n | n m l |
c
subroutine awritef(c,val) character c.32
write(*,lOOO)c,val 1000
format(/,Ix,a32,f20.6) return end
c c subroutine a w r i t e i ( c , i v a l ) character c.32 write(*,lOOO)c,ival i000 f o r m a t ( / , i x , a 3 2 , i 2 0 ) return end c c subroutine sampnt(xval,yval,zval) c c ....................................................................... c subroutine samples random points on surface c (in an e f f o r t to s i m u l a t e r a n d o m w e l l l o c a t i o n data) C .......................................................................
C common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
c data
iseed
/1/
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c ..... generate the random point and evaluate c xval - ran(iseed)*xdif + xmin yval - ran(iseed)*ydif + ymin zval - z f u n c ( x v a l , y v a l ) c c . . . . . r e t u r n to m a i n p r o g r a m c return end c
-- . . . . . . . . . . . . . . . . . . . . . . . . . . .
z
C ~ = m m m m m m m m m m m m m m m m n m m m m m m m m m m m m m m m m m m l m m m m n m m m m m m m m m m n m m m m m m m m ~ m m ~ = = a
c subroutine samlin(xval,yval,zval,distsam, iline,ishtpnt) c c ....................................................................... c s u b r o u t i n e s a m p l e s p o i n t s at r e g u l a r i n t e r v a l s a l o n g r a n d o m lines c (in an e f f o r t to s i m u l a t e s e i s m i c l l n e data) C .......................................................................
C common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
c data
iseed
/1/
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c . . . . . t e s t if we n e e d to s t a r t a n e w s a m p l i n g l i n e c if ( i s h t p n t . l e . O .or. i s h t p n t . g e . n s h t p n t ) then c c ..... new line--get initial and terminal points c iline - iline + 1 I00 xval - ran(iseed)*xdif + xmin yval - ran(iseed)*ydif + ymin xend - ran(iseed)*xdif + xmin yend - ran(iseedl*ydif + ymin xlindif - xend - xval ylindiF - yend - yval Ishtpnt - 0 c c . . . . . f i n d o u t h o w m a n y p o i n t s w i l l be a l o n g t h i s l i n e c r e j e c t l i n e if it h a s l e s s t h a n 20 o r m o r e t h a n 500 s h o t p o i n t s c nshtpnt - int(sqrt(xlindif**2 + yllndif**2l/distsam) if ( n s h t p n t . l t . 5 .or. n s h t p n t . g t . 5 0 0 ) g o t o I00 c c . . . . d e t e r m i n e t h e i n c r e m e n t s in x a n d y d i r e c t i o n s
Creating a benchmarking data set
147
c xinc yinc
- xllndlf/nshtpnt - yllndlf/nshtpnt
c endif c c ..... determine c
the
z value
at the
selected
zval - zfunc(xval,yval) c c ..... set x and ¥ values for next c along this line c xval - Kval + xinc yval - yval + yinc ishtpnt - ishtpnt + 1 c c ..... return to main program c return end c
point
point
end
increment
number
of points
c subroutine
samclu(xval,yval,zval,srad,nclupnt,
icluster,lclupnt)
c c c
subroutine samples points (in a n e f f o r t t o s i m u l a t e
C - - - - - - "
. . . . .
" - - - - "
clustered ell field
. . . . . . . . . . . . . . . . . . .
"
. . . . .
"
around a random point distributions of well
. . . . . . . . . . . . . . . . . . . . . .
"
data) . . . . .
"------
c common
/range/
xmln,xmax,xdif,ymin,ymax,ydlf
c data C
iseed
. . . . . . . . . . . . . . .
/I/ "
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
c . . . . . t e s t if w e n e e d t o s t a r t a n e w s a m p l i n g c l u s t e r c if (iclupnt.le.O .or. i c l u p n t . g e . n c l u p n t ) then c c ..... new cluster--get the center of the cluster c xcenter - ran(Iseed)*(xdlf - 2.*srad) + (xmln + erad) ycenter - ran(Iseed)*(ydlf - 2.*srad) + (ymln + srad) xval - xcenter yval - ycenter Iclupnt ~ 0 Icluster - Icluster + 1 else c c ..... still the same old cluster c xval - ran(Iseed)*(2.*srad) + (xcenter - srad) yval - ran(iseed)*(2.*srad) + (ycenter - erad) endif c c ..... determine the z value at the selected point c zval - zfunc(xval,yval) c c ..... increment the number of points in this cluster c iclupnt - iclupnt + 1 c c ..... return to main program c return end c C
subroutine
eamgrd(xval,yval,zval,xinc,yinc, ix,iy,nx,ny)
c C
c c C
. . . . . . . . . . . . . . . . . . . . .
"
. . . . .
"
. . . . .
subroutine samples points (in a n e f f o r t t o s l m u l a t e . . . . . . . . . . . . . . . . . . . . .
"
. . . . .
"
. . . . . . .
"
. . . . . . . . . . . . . .
on a regular grid regular drill hole
" - - - - ' - - - - "
. . . . . . . .
-- . . . . .
------
c common
/range/
xmin, xmax, xdif, ymln, ymax, ydlf
--
. . . . . . . . . . . . . .
sampling . . . . . . . . . . .
methods) "
. . . . . . . .
148
M . A . H£RKOMM~R c data iseed /i/ c ............................................. c ..... test if we need to start a new grid c if ( i x . l e . O .or. ix.ge.nx) then
- ........................ row
c c ..... new grid row--get the center of the cluster c xval - xmin yval = ymin + iy*yinc ix = 0 iy = iy + I else c c ..... still the same old row c xval = xmin + ix*xinc endif c c ..... determine the z value at the selected point c zval = zfunc(xval,yval) c c ..... increment the grid in the x direction c ix = ix + 1 c c ..... return to main program c return end c C function
zfunc(x,y)
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
zfunct
returns
the
z-value
at
the
x,y
c c c c c c c c c
The surface was generated by random numbers o n X, ¥, a n d Z. The X and ¥ values of the random numbers were in the range of -i.0 to 3.0, and the range of the Z values were in the range of 0.0 to 2.0. Since the test map area considers X and ¥ values sampled only in the range of 0.0 to 2.0, there are no significant "edge effects" in the surface.
c c c c c
and
A noise factor is sumulated by an intersecting at 90.0 degrees with an
amplitude
Z where
pi
of
0.05.
0.025*sln(X
* 50
-
3.14...,
C . . . . . . . . . . . . . . . . . . . .
and
The
sin
is
interference a wavelength
surface
* pi)
location
+
* sin(¥
triggnometric
the
surface.
wave pattern of 0.04 units
is d e f i n e d
0.025
the
on
as: * 50 Sine
* pi) function.
-- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c real*8
pi
parameter
(pi
-
3.141592653589793238462643)
c • common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C zfunc 1 1 2 2 2 3 3 3 3 4 4 4 4 4 5 5 5
+ + + + + + + + +
-.2396609820441817d+02 0.ii07854299397272d+03 0.7971325018951808d+02 0.1073037791619365d+03 0.2724598506594927d+03 0.1754729851507735d+02 0.1331986343733640d+03 0.2608997594233118d+03 0.2012544672059821d+03 0.9876180255358249d+02 0.2296436705657264d+03 0.7478698842146132d+02 0.2424236688359301d+03 0.8488048414202117d+02 0.6543174966233367d+02 0.6674810997854953d+02 0.2848169230914072d+03 0.8317310483769853d+02
* * * * * * * * * * * * * * * * *
x y x**2 x x**3 x**2 x x**4 x**3 x**2 x x**5 x**4 x**3
* y y**2 * y * y**2 y**3 * y * y't2 * y**3 y**4 * y * y**2
Creating a benchmarking data 5 5 5
+ 0.7374896027039002d+02 - 0.1719491617650300d+03 - 0.2337205082972909d+01
* * *
x**2 X
6 6
+
0.4096914036691392d+02 0.1800067385766880d+03
* *
6
+ 0.4921510167702959d+02
*
6 6 6 6 7
+ + + +
* * * * *
x**6 x**5 x**4 x**3 x**2 X
7 7
- 0.4307000132325771d+02 - 0.4599078834361538d+02
7 7 7
+ 0.3292136704825329d+02 + 0.1942510235575868d+02
+ + + + + +
7 7 8 8 8 8 8 8 8 8 8
c c ..... comment
0.1009944820808015d+03 0.5909168187693680d+02 0.5728848597492885d+02 0.1623304513282438d+01 0.2736188709320656d+02
this
out
if no n o i s e
149
* y**3 * y**4 y**5 * * * * *
y y**2 y**3 y**4 y**5 y**6
* * * * * *
y y**2 y**3 y**4 y**5
x**7
* x**6 * X't5 * x**4 * x**3 * x**2 * x
0.3905609419940548d+02 0.2226225744026893d+01 0.5349349866537782d+01 0.4342582899349250d+01 0.2926227853887345d+01 0.8844803734347037d+01 0.1270059588163441d+01 0.7180735683607920d+01 0.1517592187117766d+01 0.5429985121330023d+01 0.2320392475634378d+01 0.1322757693769944d+01
set
*
* * * * * * * * *
y**6 y**7
x**8 x**7
* y
x**6
* y**2
X**5 X**4 X**3
* y**3 * y**4 * y**5
X**2
* y**6
X
* y**7 y**8
is d e s i r e d
on the surface
C
zfunc - zfunc
+ 0.025*sin(x*50.*pi)
+ 0.025*sin(y*50.*pi)
C return end C a n m m m | m | m m m | n m m m l m | l | | m m m m |
real*4
function
|
|
|
m
m
|
|
l
l
m
|
|
|
|
| | | m u m
| | m
m l m | i m m m n m
|
| | |
m | l u | m
ran(ieeed)
C-----C
C C C C C C C C C C C C C . C C
Uniform random number generator returns • random number between 0.0 a n d 1,0 u s i n g t h r e e l i n e a r c o n g r u e n t l a l g e n e r a t o r s . T h r e e g e n e r a t o r s a r e u s e d to so t h a t t h e v a l u e s c a n be s h u f f l e d to a v o i d s e q u e n t i a l c o r r e l a t i o n s . O n e g e n e r a t o r is u s e d to g e n e r a t e t h e m o s t s i g n i f i c a n t p a r t of t h e o u t p u t . The second is for t h e l e a s t s i g n i f i c a n t part, a n d t h e t h i r d is to c o n t r o l the shuffling routine. seed value (initialize
-
ieeed (non-destructive) generator with a seed
value
< zero,
i.e.,
a d a p t e d from: N u m e r i c a l R e c i p e s , T h e A r t of S c i e n t i f i c C o m p u t i n g W. H. Press, et el, 1986, C a m b r i d g e U n i v e r s i t y Press, F U N C T I O N RAN1, p g 196, 197
-i)
NY,
818 pp
C
c . . . . . s e e d n u m b e r for r a n d o m n u m b e r g e n e r a t o r integer*4 iseed c ..... three linear congruentlal random number generators i n t e g e r * 4 ml, Is1, icl, m2, Is2, ic2, m3, is3, ic3 r e a l * 4 rml, rm2 p a r a m e t e r (ffil - 2 5 9 2 0 0 , Is1 - 7141, Icl - 54773, rml - 1.0/ml) parameter (m2 - 1 3 4 4 5 6 , Is2 - 8121, ic2 - 28411, rm2 - 1.0/m2) p a r a m e t e r (m3 - 2 4 3 0 0 0 , is3 = 4561, Ic3 - 51349) c . . . . . d a t a a r r a y of 97 r a n d o m n u m b e r s r e a l * 4 r(97) integer*2 J c . . . . . f I T s t t i m e in f l a g to m a k e s u r e r o u t i n e is i n i t i a l i z e d logical*l first d~ta first /.true./ c. . . . .
initialize i f iseed < 0 o r f i r s t if (iseed.lt.O .or. first) then f i r s t - .false. c ..... seed the first routine Ixl - m o d ( i c l - iseed, ml) ix1 - m o d ( l a 1 * i x l + Icl, ml) c . . . . . u s e it to s e e d t h e s e c o n d ix2 - m o d ( i x l , m2) ix1 - m o d ( i a l * i x l + icl, ml) c ..... now seed the third routine
time in
M. A. HERKOMM~R
150
Ix3 = ~ o d ( I x l , 13) c . . . . . g e n e r a t e t h e I n t l a l s e t of r a n d o m n u m b e r s do i00 j - i, 97 Ix1 - m o d ( i a l t i x l + Icl, ml) ix2 - m o d ( l a 2 * I x 2 + Ic2, m2) r(J) = ( f l o a t ( I x l ) + f l o a t ( i x 2 ) * r m 2 ) * r m l 100 continue endif c c ..... generate a random number ix1 = m o d ( i a l * i x l + icl, ml) ix2 - m o d ( i a 2 * i x 2 + ic2, m2) ix3 = m o d ( i a 3 * i x 3 + ic3, m3) c . . . . . u s e t h e t h i r d s e q u e n c e to g e t an i n t e g e r b e t w e e n 1 a n d 97 J - 1 + (97*ix3)/m3 if ( j . g t . 9 7 .or. J.lt.1) s t o p ' e r r o r in R A N f u n c t i o n -- a b o r t ' c c . . . . . r e t u r n t a b l e e n t r y a n d r e f i l l it r a n = r(j) r(j) = ( f l o a t ( i x 1 ) + f l o a t ( i x 2 ) * r m 2 ) * r m l return end
APPENDIX 2
1 2 3
0 0 0
4
0
5 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9
0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3
POO000Ol P0000002 P0000003 PO000004 PO000005 P00OO006 LOO00OO1 LO000001 LOOOOO01 LO000001 LO000001 LO000001 LO000001 LO000001 LO000002 LO000002 LO000002 L0000002 LO000002 LO000002 L0000002 C0000001 C000OOOI C0000001 C0000001 C0000001 C0000002 C0000002 C0000002 C0000002 C0000002 C0000003 C0000003 C0000003 C0000003 C0000003 G 1 G 2 G 3 G 1 G 2 G 3 G 1 G 2 G 3
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 1 1 2 2 2 3 3 3
0.2544519E-01 0.8582137 1.933947 1.313418 1.549796 0.3276081 1.742353 1.741360 1.740368 1.739375 1.738382 1.737390 1.736397 1.735404 0.7317179 0.7545426 0.7773674 0.8001921 0.8230168 0.8458415 0.8686662 0.7380269 0.7180122 0.7158826 0.8859564 0.6519554 0.7244509 0.7319114 0.7742820 0.6956647 0.8219379 0.6584604 0.5796580 0.6865289 0.6465911 0.7001495 0.0000000E+00 1.000000 2.000000 O.0000000E+O0 1.000000 2.000000 O.O000000E+O0 1.000000 2.000000
1.729282 1.045147 0.9790166 1.718119 1.251619 0.5064837 0.6506258E-01 -7.321241 1.980748 1.935788 0.8993605 4.638555 1.229443 1.437661 1.124652 1 624052 1.019861 1 891277 0.9150692 2 035509 0.8102778 2 068561 0.7054864 2 035987 0.6006950 1 805792 0.4959036 1.574549 0.4166911 -1.294396 0.5195094 -0.3699233 0.6223276 0.3944130 0.7251459 1.099516 0.8279641 1.585529 0.9307824 1.754326 1.033601 1.826368 1.372225 0.2453439 1.451359 0.8855955E 1.486418 0.1062388 1.388239 0.6747203 1.444965 -0.1519847 0.2958233 -1.816824 0.2769058 -2.252512 0.1810041 -4.318176 0.1569328 -2.734675 0.2612939 -4.452825 0.6370473 3.206604 0.5680976 4.300440 0.7255970 3.133302 0.6737955 3.560420 0.7438686 2.928964 0.0000000E+00 -23.96610 0.0000000E+00 -8.737357 O.0000000E+00 9.755122 1.000000 0.1293079 1.000000 1.462967 1.000000 -0.4603546E 2.000000 5.094895 2.000000 2.628688 2.000000 2.868078
0098-3004/88 $3.00 + 0.00 Copyright ~ 1988 Pergamon Prt~ pie
PROCEDURES FOR CREATING A BENCHMARKING DATA SET MARK A. HERKOMM'ER PETROSPEC Computer Geology and Geophysics. 4005 Burning Tree Lane. Garland. TX 75042. U.S.A. (Received 6 May 1987) Abstract--The testing of mapping (and other spatially oriented) software usually requires a test-data set. When benchmarking a program, it is worthwhile to use a data set that simulates the various types of data distributions observed in the real world: random (random wells), cluster (oil fields), grid (mineral surveys). or linear (seismic surveys). A FORTRAN V program, BENCHMARK, is presented that samples a 2-D speatial function with these data distributions. Key Words: Data distribution. Software testing, Spatial analysis.
INTRODUCTION When tcsting ncw software, it may be useful to have a benchmarking data set with known characteristics. The program presented in this paper, BENCHMARK, is designcd to produce data sets of the form Y, Y, Z. Data sets of this type are used for gridding and contouring, 2-D fihcring, and othcr spatially oricntcd applications. BENCHMARK was crcated originally to assist in testing of the accuracy of computer-contouring programs. Computer-mapping programs employ many different algorithms for handling the spatial input data. Many of these are designed for special data distributions. Broadly speaking, there are three basic data distributions: random, clustered, and regular (Davis, 1973). When benchmarking software, it is worthwhile to use data sets that reflect these types of data distribution, Each of the aforementioned data distributions has a real world analog in mapping applications. The random distribution is similar to random well locations, as seen in infantile or mature petroleum basins. The clustered distribution in the oil-well field groups are observed in juvenile to middle-age petroleum provinces. The regular grid sampling method is used in mining operations. A special (fourth) situation of the clustered distribution is data clustered along lines. Seismic survey data are a well-known example of this linear type of data distribution.
THE BENCHMARK PROGRAM B E N C H M A R K is a F O R T R A N V program (Appendix 1) designed to simulate the four types of data distributions sampling from a known 2-D surface. The input to the program is an output filename to receive the data, the range to sample the X and Y data from type of sampling mode (random, line, clustered, or regular), and the number of sample points. The Z values are determined by evaluating a single. continuous 2-D function at each of the selected X 139
and Y locations in the "map" area. For the function defined with B E N C H M A R K , the map area is 2 units by 2 units square. The minimum X a n d Yvalues used are 0.0, making the maximum X and Y values 2.0. For the random well-!ocation simulation, X and Y values were located in the range of 0.0-1.0 using RAN, a pseudorandom number generator (Press and others, 1986). They are scaled to the preselected size of the map area. The X and Y values then are put into the function Z F U N C and a Z value calculated. The seismic survey shotpoint simulation data set is better described as regularly spaced control points gathered along random lines. The initial and terminal points (X, Y pairs) of the line are random numbers. The length of the line is calculated and if it is determined that there would be < 5 or > 500 shot points on the line, it is rejected and new initial and terminal points are selected. After a line is accepted, the initial point of line is written to the data set as shot point No. 1. Subsequent points are offset from the previous point by a user specified distance along the line. Points are added along the line until it is populated fully, at which time a new line is selected. The oil-field grouping simulation data set have the characteristic of being clustered. The clusters are synthesized by first selecting a center for the cluster at random. The center of the cluster is cluster point No. I. Then the balance of the control points are determined at random having a deviation of no more than the search distance in X and Y cluster's center. The general shape of a cluster is square. The mineral-survey grid simulation is a rectangular grid of X and Y locations. The grid is created from dividing the Xand Y ranges into uniform sample intervals as determined by the number of samples along the X and Y axes. Samples are taken from the surface based on the X and Y locations. Output from B E N C H M A R K is a file containing the sampled surface. The record format of each distribution is similar so different distributions can be mixed. The general output record format is:
140
M.A.
HERKOMMER
Z.O,
x 1.9 1.6 x
2.6 X 1.9 1.6X X
1.8
1.0 -0.6 0.0 x x x-0.8
1.6
~.,~ 1.4
1.2
-2-2x
-o.fV
-1.5x
-xl.9
x 0.0 x
-1.4 x-1.3
-2.3 -2.5 x
O.9 X 0.4
~.9
0.Tx
x
4.6 X
x 0.5 2.2 x
1.7 x x
5.9 x 5.9 x
~7 4.0 x
-O.8 x
2.0 x
4.2 x
-3.3 x
0.8
0.0
-2.1x -4.4 Ix 0,2
-2.2 1 0.4
I~ 0,6
2.0
1.3x
x2.0 0.8 x
-2.0 x
x
2.2 x x 1.3
-3.5
x-1,3
3.5 x 4.8 x
x-5.2
x-3-2 x 0.2
xi.4 -O.2 x
-3.1 x
-2.9 x
6.3 x
2.4 x 2.3 x xl.8 2.2 x 2.0 1.9 x # 1.~ x
1.8 x
0.9 x
4.t
6.8 X
-0.1 xX 0.3
2,2 x
0.8
0.4
0.2x
x 1.1
1.8
4.6
1.4 1.Z x
1.4
2.0
0.6
2.1 X
x1.1
0.1 x 1.6~
x
2.2 2.1 X X 1.8 X
0.8 x
1.7 x 1.0
2.0 X
1.6 x
x -5.7
-6.9 x
-6.3 x
4,9 x 3.~x
-Z6 x
-5.8 x
x -7-3 I
1
l
1
0.8
1.0
1.2
x -3"91
1.4
1.6
I 1.8
I 2.0
Figure I. Random distribution--simulated well locations. pointtlsymblltyllclass]lindiv[l--x value--II--)' valuc--lJ--z value value--l, ] where point is the sample point number, symb is the symbol code, ty is the type, class is the group number of the particular type, indiv is the individual number withing the particular group, x value is the X location, ). value is the Y location, and z value is the sampled surface value at the determined location. Appendix Z shows an output file from B E N C H M A R K that contains all four distribution types. The random distribution record format is (2i6, 3x, 'P', i7.7, 7x, 3glS.7e2). The data values it contains are: point number, symbol code 'P' type, point number (zero padded), X, Y,.and Z. Figure ! shows a map with 100 random locations posted. The linear distribution has the format (2i6, 3x, 'L', i7.7, i7, 3g15.7e2). These data values are: point number, symbol code, 'L' type, line number (zero padded), "shot" point number on the line, X, Y, and Z. Figure 2 shows a map with 100 simulated seismic line data with a sample rate of 0.10 units along each line. The cluster distribution record format is (2i6, 3x, 'C', i7.7. i7, 3g15.7e2). The data values on its record are: point number, symbol code, 'C' type, cluster number (zero padded), cluster point number in cluster, X, Y, and Z. Figure 3 shows a map with 100 clustered data points with each cluster having 20 points and a search distance of 0.15 units. The regular grid distribution has the format: (2i6, 3x, 'G', 2i7, 3g15.7e2). These data values are: point number, symbol code. "G' type. X grid node index, Y
grid node index, X, Y, and Z. Figure 4 shows a map with a grid of 10 x 10 grid nodes. T I l E A R B I T R A R Y S U R F A C E U S E D BY BENCIlMARK
The arbitrary surface sampled by BENCHMARK is a single, continuous 2-D function. It is the sum of two discrete functions; an 8th-order polynomial and interference wave forms intersecting at an angle of 90.0- (Appendix i). Users of BENCHMARK can substitute their own function if they wish. Figure 5 shows the surface contoured on a grid of 101 x 101 grid nodes. The 8th-order polynomial part of the function was constructed in several stages. First, a set of 45 random X, Y, and Z coordinates was generated. The Xand Y values of the random numbers were in the range of - 1.0-3.0, and the range of the Z values were in the range of 0.0-2.0. The random numbers then were fitted with an 8th-order polynomial function by the method of least squares, the coefficients of which are shown in the program listing (Appendix I, subroutine ZFUNC). Because the test map area considers X and Y values sampled only in the range of 0.0-2.0, there are only minor "'edge effects" in the surface. A simple interference wave form function was added to the 8th-order polynomial surface to simulate noise and lower the degree of autocorrelation within
Creating a benchmarking
d a t a set
141
2.0
1.8
1.6
1.4
1.2
1.0
1.0 Xl. 0 x -2.1
2.2 x 1.0 x
2.4 x
1.9 x
x2.0
1.6 x 1.3 x -2.3 x 1.2 x 1.6 -2.1 0.9 x x x..1. 5 0.6x 1.9 x x x 2.2 _ 0 . -9-0 . ; ¢ 0.5 x 2.4 x 0.2 x 0.1 ~0.7 x2.5 x 1.1 xO.9 x 1.7 1.5 1.8 2.4x ~1.4 x x 2.3 .62 xl.O x 6~ ~(8 1.8x 2.1 X3. x
4.0 ~
0.8
1.6 ~
4.9 x 5.3 x
0.6
-0.2 X -1.3 x ~Om4 -2.3 x x 1.3 -3.1 x x - 3x. 6
0.2 -0.3 I~ 0.2
3.1 x -3.7 1.5 x x_3.9 ~5.7. -1.8 x 0.0 x xO.O 0.8 x -1.4 x lzr I l 0.6 0.8 0.4
1.4 x 0.5 1.6 x x 1.9 0.9 x x 1.4 2.0 x
2.41"~ _,11"" 2., 2.~2~.
1.1 x 0.4 x
5.1 x 4.3 x
0.4
0.0
2.4 x
Z.2 x
I 1.0
x 2.02.0 1.6 X x xl.5 1.8 0.8 x X 0.6 x 1.6 x - 0 . 6 x x- 0 . 2 -1.9 -2-0 x x -3.5 -2.3 x x -5.0 -3.2 x x -6.4 x -3.8 x ! I [ 1.2 1.4 1.6
F i g u r e 2. L i n e a r d i s t r i b u t i o n - - s i m u l a t e d
I 1.0
I 2.0
seismic survey.
2.0
1.0
1.6
0.4 x -0.8 -1.8d~1.7 X x - 0 . 9
~, x-2"~-2.,
-oi2
-2.0 1.4
1.3 x-O.7 x x-0.2 -0.4 x 0.5
-1.5 x -I.1
1.2
1 . ~ "1 2.3
1.0
1.2 x 2.3 x x 1k4 xl.6 2.1 x 5.4 :3.7 1.3 1.0 x 2.0 2.0 5.5#x5.6 x x x x x x 1.2 Xl.9 x x3.7 - 0 . 4 - 0 . 9 0.3 12~k 1.6 x5-6 3.9 -0.5~F x 0.SXoX91.3 x 6x8 5x'O x -I 2 x -2.3 -27 1.4 2.3 1.1 _.. -2.9 x ;(" x x 1.9 x x , . 5 " ~ X x ..x32 1 4 2"x32.6 x2.9 ; ~ 0.6 -3.4 ' " g. ~¢ X3.4 • x-O.9 X , = -0.5 16~2.4"29 ~z'0"3-"'~ - ~ . 4 "2.3 xx x 3.7 -1.4 -4.3X~4"? x x x 3.4 ~35Xx -4.7 -1.5 1.1 XxX4.2 • x - 4 . 6 -7.4 4~0
0.8
0.6
0.4
O.Z
-4,8
0.0
I 0.2
1 0.4
I 0.6
-4.7. [ 0.8
I 1.0
[ 1.2
F i g u r e 3. C l u s t e r d i s t r i b u t i o n - - s i m u l a t e d
I 1.4
oil fields.
[ 1.6
I 1.8
I 2.0
M. A. I-IERKOMM]LR
142
2.0x5.1
2.9
x
3.1
3,0
2.8
2.6
2.5
1.9
1.2
2.9
1.0 2.5
0.5
1.1 X
1.7 X
2.0 X
2.3 X
2.4 X
1,9 X
1.3 X
2.0
1.6 I-1.9 X
-2.3 X
-1.0 X
0.1 X
0.9 X
1.5 X
2.0 X
1.'/' X
1.1 X
1.4
1,4 ~ 3 . 8 X
-2.6 X
--0.9 X
0.1 X
0.8 X
1.5 X
2.2 X
2.1 X
1.3 X
0.2
1.2 ~ 1 . 8 x
0.1 x
1.4 X
1.5 X
1.4 X
1.7 X
2.4 X
2.5 X
1.7 X
-0.3
1"0!.2
4.4x
4.7 X
3.1 X
1.6 X
1.0 X
1.6 X
2.:) X
2.0 X
-0.4
6.4 X
3.2 X
0.0 X
-1,3 X
-0.6 X
0.7 X
1.9 X
2.1
5.0 x
1,4 x
-2.9 x
-4,8 x
-3.8 x
-1,2 x
1.8 x
4.9
2.9 X
-1.4 X
-6.1 X
-0.1 X
-6.4 X
-2.6 X
2.2 X
7.9
l°i ° o4
r2i ' 0.6
I-7i4 o.0
%31 ,6
~o ,o
98~ 20
1.8
2.0
x
x
v
K
X
X
X
X
X
t
0.8 t . 6 x
7.2 x
0.6 1 . 9
6.5
0,41
x
0.1
1.6
0.2i
0.0
x
x
x
x
x
X
x
x
x
I - 24.0 -5.6 x ix 0.2
j - 9 i, ~ - ~ 8 , ,o ,2 ,4
Figure 4. Grid distribution--simulated mineral survey.
2.0 1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.2
0.4
0.6
0.8
1.0
L2
1.4
1.6
Figure 5. Benchmark surface used to sample distributions.
Creating a benchmarking data set the surface (Harbaugh, Doveton, and Davis, 1977). Adding the wave function makes it more difficult to model the surface. Only those data sets using at least 10,000 control points would be able the "see" this overlain interference surface. To create the "noise", two wave patterns are intersected in phase at 90.0 ° with a wavelength of 0.04 units and an amplitude of 0.05. This is seen as the last line in the equation for the function in Appendix I.
143
F O R T R A N may have difficulty with the number of continuation lines in ZFUNC. This can be overcome by breaking the function into pieces and summing the parts. The program requires approximately 64 K bytes of memory to execute. Although run times may vary, on an IBM PC/AT, 10,000 points can be generated in about 4 min. REFERENCES
COMPUTER REQUIREMENTS FOR USING BENCHMARK B E N C H M A R K was written originally on a VAX 11/780 in ANSI standard F O R T R A N V with minor VAX extensions. The version presented here has been ported to an IBM PC/AT and compiled with R y a n McFarland FORTRAN. Some other versions of
Davis, J. C.. 1973. Statistics and data analysis in geology: John Wiley & Sons, New York. 550p. Harbaugh. J. W., Doveton, J. H., Davis, J. C.. 1977, Probability methods in oil exploration: John Wiley & Sons, New York. 267p. Press, W. H., Flannery, B. P.. Teukolsky. S. A., and Vetteding, W. T., 1986. Numerical recipes: Cambridge Univ. Press, Cambridge, 818p.
APPENDIX 1
Program BENCHMARK c cnnmmnlaiun~nnl|ummnilamammllimmlamnmnmmmaeallmul|mmN|iim|anaugmmmnulnn
c c c c c c c
program to select data points randomly from a given surface data selection is in 4 modes; random points and points along random lines, points in random clusters, and points on a grid. written
by Mark A. Herkommer
c
modified
for pc on 4-28-87
on 6-6-83
Cm.minmnnmmanmmnn.uumnennminmimmnmmminnmm.wmi.mmmmmmmm...mnnnm.iwenlnl~ c
character cyn*l,coutfil*40 integer*2
cbell
parameter
(cbell - 7)
c c common /range/ C
. . . . . . . . . .
"
. . . . . . .
xmin, xmax, xdif, ymln, ymax, ydlf
"
100
wrlte(*,1100) &'Welcome to BENCHMARK -- a Fortran V language program used to', &'generate data sets for benchmarklng purposes. Given a user', &'defined 3-d function, benchmark samples the surface at random', &'locatlons. The program uses 4 sampllng techniques, selectlng', &'random points (mode-0), points along random lines (mode-l),', &'points clustered around a random point (mode-2), or points', &'arranged in the pattern of a regular grid (mode-3). For', &'"polnts-along-random-lines" mode, the sampling distance along', &'llne selected must be specified. Lines will have more than', &'5 and less than 500 points. For "polnts-clustered-around', &'a-random-polnt", the cluster distance and number of points per', &'per cluster must be specified. For "points-on-a-regular-grldO',' &'an x and y increment must be selected. The range for the data', &'range to be selected is input as xmin, ymin, xmax, and ymax.' 1100 format(6x,a) wrlte(*,1100) ' ', &'The output file from benchmark is a sequential data file with', &'each record having the followlng format:'," ', &'point number (16), llne/cluster Id (aS), x (g15.7),'// &'y (g15.7), z (g15.7}'
c c ..... get parameters for data point selection c write(*,1200} 1200 format(//," enter the output filename read(*,'(a40)')coutfll open(unlt - 10,file - coutfil,status - 'unknown')
: ')
144
M . A . H£RKOM.'~R
xmin - aread('enter the x-mlnimum ymln - aread('enter the y-minimum xmax - aread('enter the x-maximum ymax - aread('enter the y-maximum xdif - xmax - xmin ydlf - ymax - ymin i£ ( x d l f . l e . O . O .or. y d l f . l e . O . O ) & stop 'minlmum/maximum r a n g e error"
=') =') =') =')
c c ..... sampling mode c 200 wrlte(*,1210) 1210 f o r m a t ( / , ' e n t e r t h e s a m p l i n g m o d e ; ' , t 3 0 , ' O - random points',/ & t30,'1 - points along random lines',/ & t30,'2 - points clustered around a random point',/ & t 3 0 , ' 3 - p o i n t s on a r e g u l a r g r i d ' , / ) c mode = aread('enter the sampling mode c if(mode.eq.O) then npnt a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eg.1) t h e n d i s t s a m - a r e a d ( ' e n t e r the s a m p l e d i s t a n c e a l o n g l i n e s npnt = a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eq.2) t h e n srad a r e a d ( ' e n t e r the s e a r c h d i s t a n c e for all c l u s t e r s n c l u p n t - a r e a d ( ' e n t e r the n u m b e r of p o i n t s p e r c l u s t e r npnt a r e a d ( ' e n t e r the n u m b e r of p o i n t s to be s a m p l e d c e l s e i f (mode.eq.3) t h e n nx a r e a d ( ' e n t e r n u m b e r of s a m p l e s a l o n g t h e x a x i s ny a r e a d ( ' e n t e r n u m b e r of s a m p l e s a l o n g t h e y a x i s If (nx.le.l .or. ny.le.l) & s t o p ' s e l e c t s a m p l e r a t e s g r e a t e r t h a n I' x i n c - x d i f / ( n x - I) y i n c - y d i f / ( n y - i) npnt - nx*ny c else write(*,'(/,a)') ' p l e a s e s e l e c t O, I, 2, or 3...' g o t o 200 endif isymb a r e a d ( ' e n t e r t h e s y m b o l c o d e for t h e p o i n t s c c ..... write out selection parameters c wrlte(*,1250)coutfll 1250 f o r m a t ( / / / / ' d a t a p o i n t s e l e c t i o n p a r a m e t e r s ' , / & , ............................... ,,/ &' o u t p u t f i l e n a m e : ',a40) call awritef('x-minlmum - ',xmin) call awritef('y-mlnimum - ',ymin) call awrltef('x-maximum = ',x max) call awritef('y-maximum = ',y max) if(mode.eq.O) then write(*,1260) 1260 format(/,' sampling mode - ', & 'random points') e l s e i f (mode.eq.l) t h e n write(*,1270) 1270 format(/, ' sampling mode - ', & ' p o i n t s a l o n g r a n d o m lines') call awriteE('sample distance along lines - ',distsam) ishtpnt - 0 illne - 0 c e l e e l f (mode.eq.2) t h e n write(*,1280) I ', format(/, ' sampling mode 1280 & 'points clustered around a random point') call awrltef('search d i s t a n c e for a l l c l u s t e r s = call awrltei('number of p o i n t s p e r c l u s t e r = iclupnt - 0 icluster - 0
elseif
(mode.eq.3)
then
',srad) ',nclupnt)
=')
=')
=') =')
=') =') =')
-') =')
=')
Creatingabenchmarkingdata
1290
145
~t
write(*,1290) format(/, ' s a m p l i n g m o d e " ', & . ' p o i n t s o n a r e g u l a r grid') call a w r i t e i ( ' n u m b e r of s a m p l e s a l o n g x a x i s call a w r l t e i ( ' n u m b e r of s a m p l e s a l o n g y a x i s ix - 0 iy - 0 endif call s w r l t e i ( ' n u m b e r of p o i n t s to be s a m p l e d call a w r l t e i ( ' s y m b o l c o d e for p o i n t s -
- ",nx) - ',ny)
',npnt) ',isymb)
c c . . . . . c h e c k to m a k e s u r e d a t a is ok c write(*,1300) 1300 f o r m a t ( / , ' are t h e s e v a l u e s ok
&
i,lsymb, icluster,lclupnt,xval,yval,zval
c e l s e i f (mode.eq.3) t h e n c c . . . . . p o i n t s on a r e g u l a r g r i d c c a l l s a m g r d ( x v a l , y v a l , z v a l , x l n c , y l h c , l x , ly,nx,ny) write(10,'(216,3x,''G'',217,3glS.7e2)') & i,lsymb, lx, i y , x v a l , y v a l , z v a l endif 300 continue c wrlte(*,1400) cbell 1400 f o r m a t ( I x , a l , ' s a m p l l n g complete.',//, &' d o y o u w a n t s a m p l e a g a i n w i t h the s a m e p a r a m e t e r s
C function aread(c) c h a r a c t e r c*(*) wrlte(*,1000)c 1000 f o r m a t ( I x , a , ' ') read(*,'(f20.0)')aread return
end
146
M . A . HERKOMMER c C
i
l
l
l
| a a n e m n
N i l
n m e | g a a u m
i i
e | a m m
l
nlnllemm
i
i
l
i
mm
l
i
l
l
i
mnlmmn
i i
| m
n m n n m n
| n | n m l |
c
subroutine awritef(c,val) character c.32
write(*,lOOO)c,val 1000
format(/,Ix,a32,f20.6) return end
c c subroutine a w r i t e i ( c , i v a l ) character c.32 write(*,lOOO)c,ival i000 f o r m a t ( / , i x , a 3 2 , i 2 0 ) return end c c subroutine sampnt(xval,yval,zval) c c ....................................................................... c subroutine samples random points on surface c (in an e f f o r t to s i m u l a t e r a n d o m w e l l l o c a t i o n data) C .......................................................................
C common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
c data
iseed
/1/
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c ..... generate the random point and evaluate c xval - ran(iseed)*xdif + xmin yval - ran(iseed)*ydif + ymin zval - z f u n c ( x v a l , y v a l ) c c . . . . . r e t u r n to m a i n p r o g r a m c return end c
-- . . . . . . . . . . . . . . . . . . . . . . . . . . .
z
C ~ = m m m m m m m m m m m m m m m m n m m m m m m m m m m m m m m m m m m l m m m m n m m m m m m m m m m n m m m m m m m m ~ m m ~ = = a
c subroutine samlin(xval,yval,zval,distsam, iline,ishtpnt) c c ....................................................................... c s u b r o u t i n e s a m p l e s p o i n t s at r e g u l a r i n t e r v a l s a l o n g r a n d o m lines c (in an e f f o r t to s i m u l a t e s e i s m i c l l n e data) C .......................................................................
C common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
c data
iseed
/1/
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c . . . . . t e s t if we n e e d to s t a r t a n e w s a m p l i n g l i n e c if ( i s h t p n t . l e . O .or. i s h t p n t . g e . n s h t p n t ) then c c ..... new line--get initial and terminal points c iline - iline + 1 I00 xval - ran(iseed)*xdif + xmin yval - ran(iseed)*ydif + ymin xend - ran(iseed)*xdif + xmin yend - ran(iseedl*ydif + ymin xlindif - xend - xval ylindiF - yend - yval Ishtpnt - 0 c c . . . . . f i n d o u t h o w m a n y p o i n t s w i l l be a l o n g t h i s l i n e c r e j e c t l i n e if it h a s l e s s t h a n 20 o r m o r e t h a n 500 s h o t p o i n t s c nshtpnt - int(sqrt(xlindif**2 + yllndif**2l/distsam) if ( n s h t p n t . l t . 5 .or. n s h t p n t . g t . 5 0 0 ) g o t o I00 c c . . . . d e t e r m i n e t h e i n c r e m e n t s in x a n d y d i r e c t i o n s
Creating a benchmarking data set
147
c xinc yinc
- xllndlf/nshtpnt - yllndlf/nshtpnt
c endif c c ..... determine c
the
z value
at the
selected
zval - zfunc(xval,yval) c c ..... set x and ¥ values for next c along this line c xval - Kval + xinc yval - yval + yinc ishtpnt - ishtpnt + 1 c c ..... return to main program c return end c
point
point
end
increment
number
of points
c subroutine
samclu(xval,yval,zval,srad,nclupnt,
icluster,lclupnt)
c c c
subroutine samples points (in a n e f f o r t t o s i m u l a t e
C - - - - - - "
. . . . .
" - - - - "
clustered ell field
. . . . . . . . . . . . . . . . . . .
"
. . . . .
"
around a random point distributions of well
. . . . . . . . . . . . . . . . . . . . . .
"
data) . . . . .
"------
c common
/range/
xmln,xmax,xdif,ymin,ymax,ydlf
c data C
iseed
. . . . . . . . . . . . . . .
/I/ "
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
c . . . . . t e s t if w e n e e d t o s t a r t a n e w s a m p l i n g c l u s t e r c if (iclupnt.le.O .or. i c l u p n t . g e . n c l u p n t ) then c c ..... new cluster--get the center of the cluster c xcenter - ran(Iseed)*(xdlf - 2.*srad) + (xmln + erad) ycenter - ran(Iseed)*(ydlf - 2.*srad) + (ymln + srad) xval - xcenter yval - ycenter Iclupnt ~ 0 Icluster - Icluster + 1 else c c ..... still the same old cluster c xval - ran(Iseed)*(2.*srad) + (xcenter - srad) yval - ran(iseed)*(2.*srad) + (ycenter - erad) endif c c ..... determine the z value at the selected point c zval - zfunc(xval,yval) c c ..... increment the number of points in this cluster c iclupnt - iclupnt + 1 c c ..... return to main program c return end c C
subroutine
eamgrd(xval,yval,zval,xinc,yinc, ix,iy,nx,ny)
c C
c c C
. . . . . . . . . . . . . . . . . . . . .
"
. . . . .
"
. . . . .
subroutine samples points (in a n e f f o r t t o s l m u l a t e . . . . . . . . . . . . . . . . . . . . .
"
. . . . .
"
. . . . . . .
"
. . . . . . . . . . . . . .
on a regular grid regular drill hole
" - - - - ' - - - - "
. . . . . . . .
-- . . . . .
------
c common
/range/
xmin, xmax, xdif, ymln, ymax, ydlf
--
. . . . . . . . . . . . . .
sampling . . . . . . . . . . .
methods) "
. . . . . . . .
148
M . A . H£RKOMM~R c data iseed /i/ c ............................................. c ..... test if we need to start a new grid c if ( i x . l e . O .or. ix.ge.nx) then
- ........................ row
c c ..... new grid row--get the center of the cluster c xval - xmin yval = ymin + iy*yinc ix = 0 iy = iy + I else c c ..... still the same old row c xval = xmin + ix*xinc endif c c ..... determine the z value at the selected point c zval = zfunc(xval,yval) c c ..... increment the grid in the x direction c ix = ix + 1 c c ..... return to main program c return end c C function
zfunc(x,y)
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
zfunct
returns
the
z-value
at
the
x,y
c c c c c c c c c
The surface was generated by random numbers o n X, ¥, a n d Z. The X and ¥ values of the random numbers were in the range of -i.0 to 3.0, and the range of the Z values were in the range of 0.0 to 2.0. Since the test map area considers X and ¥ values sampled only in the range of 0.0 to 2.0, there are no significant "edge effects" in the surface.
c c c c c
and
A noise factor is sumulated by an intersecting at 90.0 degrees with an
amplitude
Z where
pi
of
0.05.
0.025*sln(X
* 50
-
3.14...,
C . . . . . . . . . . . . . . . . . . . .
and
The
sin
is
interference a wavelength
surface
* pi)
location
+
* sin(¥
triggnometric
the
surface.
wave pattern of 0.04 units
is d e f i n e d
0.025
the
on
as: * 50 Sine
* pi) function.
-- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c real*8
pi
parameter
(pi
-
3.141592653589793238462643)
c • common
/range/
xmin,xmax,xdif,ymin,ymax,ydif
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C zfunc 1 1 2 2 2 3 3 3 3 4 4 4 4 4 5 5 5
+ + + + + + + + +
-.2396609820441817d+02 0.ii07854299397272d+03 0.7971325018951808d+02 0.1073037791619365d+03 0.2724598506594927d+03 0.1754729851507735d+02 0.1331986343733640d+03 0.2608997594233118d+03 0.2012544672059821d+03 0.9876180255358249d+02 0.2296436705657264d+03 0.7478698842146132d+02 0.2424236688359301d+03 0.8488048414202117d+02 0.6543174966233367d+02 0.6674810997854953d+02 0.2848169230914072d+03 0.8317310483769853d+02
* * * * * * * * * * * * * * * * *
x y x**2 x x**3 x**2 x x**4 x**3 x**2 x x**5 x**4 x**3
* y y**2 * y * y**2 y**3 * y * y't2 * y**3 y**4 * y * y**2
Creating a benchmarking data 5 5 5
+ 0.7374896027039002d+02 - 0.1719491617650300d+03 - 0.2337205082972909d+01
* * *
x**2 X
6 6
+
0.4096914036691392d+02 0.1800067385766880d+03
* *
6
+ 0.4921510167702959d+02
*
6 6 6 6 7
+ + + +
* * * * *
x**6 x**5 x**4 x**3 x**2 X
7 7
- 0.4307000132325771d+02 - 0.4599078834361538d+02
7 7 7
+ 0.3292136704825329d+02 + 0.1942510235575868d+02
+ + + + + +
7 7 8 8 8 8 8 8 8 8 8
c c ..... comment
0.1009944820808015d+03 0.5909168187693680d+02 0.5728848597492885d+02 0.1623304513282438d+01 0.2736188709320656d+02
this
out
if no n o i s e
149
* y**3 * y**4 y**5 * * * * *
y y**2 y**3 y**4 y**5 y**6
* * * * * *
y y**2 y**3 y**4 y**5
x**7
* x**6 * X't5 * x**4 * x**3 * x**2 * x
0.3905609419940548d+02 0.2226225744026893d+01 0.5349349866537782d+01 0.4342582899349250d+01 0.2926227853887345d+01 0.8844803734347037d+01 0.1270059588163441d+01 0.7180735683607920d+01 0.1517592187117766d+01 0.5429985121330023d+01 0.2320392475634378d+01 0.1322757693769944d+01
set
*
* * * * * * * * *
y**6 y**7
x**8 x**7
* y
x**6
* y**2
X**5 X**4 X**3
* y**3 * y**4 * y**5
X**2
* y**6
X
* y**7 y**8
is d e s i r e d
on the surface
C
zfunc - zfunc
+ 0.025*sin(x*50.*pi)
+ 0.025*sin(y*50.*pi)
C return end C a n m m m | m | m m m | n m m m l m | l | | m m m m |
real*4
function
|
|
|
m
m
|
|
l
l
m
|
|
|
|
| | | m u m
| | m
m l m | i m m m n m
|
| | |
m | l u | m
ran(ieeed)
C-----C
C C C C C C C C C C C C C . C C
Uniform random number generator returns • random number between 0.0 a n d 1,0 u s i n g t h r e e l i n e a r c o n g r u e n t l a l g e n e r a t o r s . T h r e e g e n e r a t o r s a r e u s e d to so t h a t t h e v a l u e s c a n be s h u f f l e d to a v o i d s e q u e n t i a l c o r r e l a t i o n s . O n e g e n e r a t o r is u s e d to g e n e r a t e t h e m o s t s i g n i f i c a n t p a r t of t h e o u t p u t . The second is for t h e l e a s t s i g n i f i c a n t part, a n d t h e t h i r d is to c o n t r o l the shuffling routine. seed value (initialize
-
ieeed (non-destructive) generator with a seed
value
< zero,
i.e.,
a d a p t e d from: N u m e r i c a l R e c i p e s , T h e A r t of S c i e n t i f i c C o m p u t i n g W. H. Press, et el, 1986, C a m b r i d g e U n i v e r s i t y Press, F U N C T I O N RAN1, p g 196, 197
-i)
NY,
818 pp
C
c . . . . . s e e d n u m b e r for r a n d o m n u m b e r g e n e r a t o r integer*4 iseed c ..... three linear congruentlal random number generators i n t e g e r * 4 ml, Is1, icl, m2, Is2, ic2, m3, is3, ic3 r e a l * 4 rml, rm2 p a r a m e t e r (ffil - 2 5 9 2 0 0 , Is1 - 7141, Icl - 54773, rml - 1.0/ml) parameter (m2 - 1 3 4 4 5 6 , Is2 - 8121, ic2 - 28411, rm2 - 1.0/m2) p a r a m e t e r (m3 - 2 4 3 0 0 0 , is3 = 4561, Ic3 - 51349) c . . . . . d a t a a r r a y of 97 r a n d o m n u m b e r s r e a l * 4 r(97) integer*2 J c . . . . . f I T s t t i m e in f l a g to m a k e s u r e r o u t i n e is i n i t i a l i z e d logical*l first d~ta first /.true./ c. . . . .
initialize i f iseed < 0 o r f i r s t if (iseed.lt.O .or. first) then f i r s t - .false. c ..... seed the first routine Ixl - m o d ( i c l - iseed, ml) ix1 - m o d ( l a 1 * i x l + Icl, ml) c . . . . . u s e it to s e e d t h e s e c o n d ix2 - m o d ( i x l , m2) ix1 - m o d ( i a l * i x l + icl, ml) c ..... now seed the third routine
time in
M. A. HERKOMM~R
150
Ix3 = ~ o d ( I x l , 13) c . . . . . g e n e r a t e t h e I n t l a l s e t of r a n d o m n u m b e r s do i00 j - i, 97 Ix1 - m o d ( i a l t i x l + Icl, ml) ix2 - m o d ( l a 2 * I x 2 + Ic2, m2) r(J) = ( f l o a t ( I x l ) + f l o a t ( i x 2 ) * r m 2 ) * r m l 100 continue endif c c ..... generate a random number ix1 = m o d ( i a l * i x l + icl, ml) ix2 - m o d ( i a 2 * i x 2 + ic2, m2) ix3 = m o d ( i a 3 * i x 3 + ic3, m3) c . . . . . u s e t h e t h i r d s e q u e n c e to g e t an i n t e g e r b e t w e e n 1 a n d 97 J - 1 + (97*ix3)/m3 if ( j . g t . 9 7 .or. J.lt.1) s t o p ' e r r o r in R A N f u n c t i o n -- a b o r t ' c c . . . . . r e t u r n t a b l e e n t r y a n d r e f i l l it r a n = r(j) r(j) = ( f l o a t ( i x 1 ) + f l o a t ( i x 2 ) * r m 2 ) * r m l return end
APPENDIX 2
1 2 3
0 0 0
4
0
5 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9
0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3
POO000Ol P0000002 P0000003 PO000004 PO000005 P00OO006 LOO00OO1 LO000001 LOOOOO01 LO000001 LO000001 LO000001 LO000001 LO000001 LO000002 LO000002 LO000002 L0000002 LO000002 LO000002 L0000002 C0000001 C000OOOI C0000001 C0000001 C0000001 C0000002 C0000002 C0000002 C0000002 C0000002 C0000003 C0000003 C0000003 C0000003 C0000003 G 1 G 2 G 3 G 1 G 2 G 3 G 1 G 2 G 3
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 1 1 2 2 2 3 3 3
0.2544519E-01 0.8582137 1.933947 1.313418 1.549796 0.3276081 1.742353 1.741360 1.740368 1.739375 1.738382 1.737390 1.736397 1.735404 0.7317179 0.7545426 0.7773674 0.8001921 0.8230168 0.8458415 0.8686662 0.7380269 0.7180122 0.7158826 0.8859564 0.6519554 0.7244509 0.7319114 0.7742820 0.6956647 0.8219379 0.6584604 0.5796580 0.6865289 0.6465911 0.7001495 0.0000000E+00 1.000000 2.000000 O.0000000E+O0 1.000000 2.000000 O.O000000E+O0 1.000000 2.000000
1.729282 1.045147 0.9790166 1.718119 1.251619 0.5064837 0.6506258E-01 -7.321241 1.980748 1.935788 0.8993605 4.638555 1.229443 1.437661 1.124652 1 624052 1.019861 1 891277 0.9150692 2 035509 0.8102778 2 068561 0.7054864 2 035987 0.6006950 1 805792 0.4959036 1.574549 0.4166911 -1.294396 0.5195094 -0.3699233 0.6223276 0.3944130 0.7251459 1.099516 0.8279641 1.585529 0.9307824 1.754326 1.033601 1.826368 1.372225 0.2453439 1.451359 0.8855955E 1.486418 0.1062388 1.388239 0.6747203 1.444965 -0.1519847 0.2958233 -1.816824 0.2769058 -2.252512 0.1810041 -4.318176 0.1569328 -2.734675 0.2612939 -4.452825 0.6370473 3.206604 0.5680976 4.300440 0.7255970 3.133302 0.6737955 3.560420 0.7438686 2.928964 0.0000000E+00 -23.96610 0.0000000E+00 -8.737357 O.0000000E+00 9.755122 1.000000 0.1293079 1.000000 1.462967 1.000000 -0.4603546E 2.000000 5.094895 2.000000 2.628688 2.000000 2.868078