1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
""" Tests for BuildGrid.c module.
To view the grids I used QuikGrid from
http://www.PerspectiveEdge.com
"""
import sys,math,random,time
random.seed(7)
global ncol,nrow,step,xmin,ymin,simple
xmin,ymin,step = 0.0,0.0,1.0
# Surface model:
simple = True # simple: (a1+b1*x)*(a2+b2*y)
# This is almost a plain.
simple = False # !simple: cos(a*x*x)*cos(b*y*y)
# In this case only left bottom part of surface
# should be believable.
def gridsizes(nodes,ratio=1.0):
# ratio = ncol/nrow
nodes = max(100,nodes)
ratio = max(0.1,min(ratio,10.0))
nrow = int(math.sqrt(nodes/ratio))
ncol = nodes/nrow+1
return ncol,nrow
ncol,nrow = gridsizes(10**5) # I used 10**3 - 10**7 nodes
ax,ay = math.pi/(4*ncol+1),math.pi/(4*nrow+1)
percent = 10.2 # how many data points w.r.t. grid nodes:
datapoints = int(0.01 * percent * ncol * nrow) or 5
print "Testing",("complex","simple")[simple],
print "surface. Grid size %d x %d (%d nodes)" % (ncol,nrow,ncol*nrow)
print "About %d datapoints (%.2f%% of nodes)" % (datapoints,percent)
# Test flag
build_only = True # no output to files
build_only = False # output to files
# Trimming distance
trimdist = 5*step # make trimming
trimdist = 0.0 # do not make trimming
if trimdist < step:
print "No trimming"
else:
print "With trimming",trimdist
def surface(x,y):
if simple:
return 100+100*(1+float(x)/ncol)*(1+float(y)/nrow)
return 100+50*(math.cos(ax*x*x)*math.cos(ay*y*y)+1)
def surfaceGen():
while True:
if trimdist < step:
x = random.uniform(-1,ncol) # some points may be
y = random.uniform(-1,nrow) # out of grid
else:
trdist = int(trimdist/step+0.5)
x = random.uniform(-1+trdist,ncol-trdist) # all the points
y = random.uniform(-1+trdist,nrow-trdist) # are inside the grid
z = surface(x,y)
yield x,y,z
def makeInputXYZ(outfile,nvals):
surgen = surfaceGen()
x,y,z = [],[],[]
fo = open(outfile,"wt") # for fromfile()
while(nvals > 0):
xv,yv,zv = surgen.next()
fo.write("%.3f %.3f %.3f\n" % (xv,yv,zv))
x.append(xv) # for fromxyz()
y.append(yv)
z.append(zv)
nvals -= 1
fo.close()
return x,y,z
# create input data file with
xyzfile = "inpdata.xyz" # input for fromfile()
xp,yp,zp = makeInputXYZ(xyzfile, datapoints)
import BuildGrid
# get file statistics
nvals,xmin,xmax,ymin,ymax,zmin,zmax,zavrg,zstnd = \
BuildGrid.filestatistics(xyzfile)
print xyzfile,'statistics. Number of values',nvals
print 'X min %9.1f max %9.1f' % (xmin,xmax)
print 'Y min %9.1f max %9.1f' % (ymin,ymax)
print 'S min %9.1f max %9.1f' % (zmin,zmax)
print 'S avrg %7.2f stnd %7.3f' % (zavrg,zstnd)
# build grid 'fromfile'
t0 = time.clock()
grid = BuildGrid.fromfile(xyzfile=xyzfile,
nx=ncol, ny=nrow, step=step, xmin=xmin, ymin=ymin,
method='Good', # or 'Best' - not implemented
trimming=trimdist, # if no trimming - full grid
unvalue=123.321, # will be used if trimming > 0
abserror=0.001, # may be == 0
relerror=0.0001) # ignores if abserror > 0
print 'fromfile():',xyzfile,len(grid)-ncol*nrow,"(%.2f sec, %.0f nodes/sec)" %\
(time.clock()-t0,ncol*nrow/(time.clock()-t0))
# put grid data to xyz file
if not build_only:
outfile = "outdata1.xyz"
t0 = time.clock()
rv = BuildGrid.tofile(filename=outfile, griddata=grid, gridtype='xyz',
nx=ncol, ny=nrow, step=step, xmin=xmin, ymin=ymin, unvalue=1234.4321)
print 'tofile():',outfile,rv,"(%.2f sec)" % (time.clock()-t0)
# build grid from xyz lists
t0 = time.clock()
grid2 = BuildGrid.fromxyz(xdata=xp,ydata=yp,zdata=zp,
nx=ncol, ny=nrow, step=step, xmin=xmin, ymin=ymin,
method='Good', # or 'Best' (not implemented)
trimming=trimdist, # if no trimming - full grid
unvalue=123.321, # will be used if trimming > 0
abserror=0.001, # may be == 0
relerror=0.0001) # ignores if abserror > 0
print 'fromxyz():',len(grid)-ncol*nrow,"(%.2f sec, %.0f nodes/sec)" %\
(time.clock()-t0,ncol*nrow/(time.clock()-t0))
# put grid to file
if not build_only:
outfile = "outdata2.xyz"
t0 = time.clock()
rv = BuildGrid.tofile(filename=outfile, griddata=grid2, gridtype='xyz',
nx=ncol, ny=nrow, step=step, xmin=xmin, ymin=ymin, unvalue=1234.4321)
print 'tofile():',outfile,rv,"(%.2f sec)" % (time.clock()-t0)
# write full exact grid
if not build_only:
exactfile = "exactdata.xyz"
t0 = time.clock()
fo = open(exactfile,"wt")
for nr in xrange(nrow):
for nc in xrange(ncol):
fo.write("%.2f %.2f %.3f\n" % (nc,nr,surface(nc,nr)))
fo.close()
print 'exact file:',exactfile,"(%.2f sec)" % (time.clock()-t0)
if build_only:
sys.exit()
# Compare Input & Output:
layer = max(5,min(ncol,nrow)/10)
def avst(left,rite,low,top):
av = st = 0.0
nv = 0
for nr in xrange(low,top):
for nc in xrange(left,rite):
d = grid[nr*ncol+nc]-surface(nc,nr)
av += d
st += d*d
nv += 1
av /= nv
st = st/nv-av*av
if st > 0.0: st = math.sqrt(st)
return av,st
print "Comparing %d x %d squares" % (layer,layer),"of grid lines"
av,st = avst(0,layer,0,layer)
print "Bottom Left: av = %8.3f st = %8.3f" % (av,st)
av,st = avst(ncol-layer,ncol,0,layer)
print "Bottom Right: av = %8.3f st = %8.3f" % (av,st)
av,st = avst(0,layer,nrow-layer,nrow)
print "Top Left: av = %8.3f st = %8.3f" % (av,st)
av,st = avst(ncol-layer,ncol,nrow-layer,nrow)
print "Top Right: av = %8.3f st = %8.3f" % (av,st)
av,st = avst(ncol/2-layer/2,ncol/2+layer/2,nrow/2-layer/2,nrow/2+layer/2)
print "Middle: av = %8.3f st = %8.3f" % (av,st)
|
