DISPATCH
_averages.py
1 # -*- coding: utf-8 -*-
2 """
3 Created on Tue Aug 07 23:05:36 2018
4 
5 @author: Aake
6 """
7 # Pythn 2/3 compatibility
8 from __future__ import print_function
9 
10 import numpy as np
11 import dispatch
12 #import types
13 
14 def aver(s,iv=0):
15  f=0.0
16  v=0.0
17  p=s.patches[0]
18  data=p.var(iv)
19  for p in s.patches:
20  dv=np.product(p.ds)
21  v=v+dv
22  f=f+dv*np.sum(data)
23  return f/v
24 
25 def rms(s,iv=0):
26  f=0.0
27  v=0.0
28  p=s.patches[0]
29  for p in s.patches:
30  if iv in p.all_keys:
31  data=p.var(iv)
32  dv=np.product(p.ds)
33  v=v+dv
34  f=f+dv*np.sum(data)
35  fa=f/v
36  f=0.0
37  for p in s.patches:
38  if iv in p.all_keys:
39  data=p.var(iv)
40  dv=np.product(p.ds)
41  v=v+dv
42  f=f+dv*np.sum((data-fa)**2)
43  return np.sqrt(f/v)
44 
45 def minloc(a):
46  ''' Return the 3-D index with the minimum of array a'''
47  i = a.argmin()
48  return np.unravel_index(i,np.shape(a))
49 
50 def maxloc(a):
51  ''' Return the 3-D index with the maximum of array a'''
52  i = a.argmax()
53  return np.unravel_index(i,np.shape(a))
54 
55 def _patches(s):
56  """ Allow a snapshot, a patch, or a list of patches """
57  if hasattr(s,'dict'):
58  pp=s.patches
59  elif isinstance(s,dispatch._dispatch._patch):
60  pp=[s]
61  else:
62  pp=s
63  return pp
64 
65 def stat(s,iv=0,i4=0):
66  f=0.0
67  v=0.0
68  pp=_patches(s)
69  p=pp[0]
70  data=p.var(iv,i4=i4)
71  d=data
72  m=maxloc(d)
73  dmin=float(d.min())
74  dmax=float(d.max())
75  rmax=[p.x[m[0]],p.y[m[1]],p.z[m[2]]]
76  m=minloc(d)
77  rmin=[p.x[m[0]],p.y[m[1]],p.z[m[2]]]
78  for p in pp:
79  if iv in p.all_keys:
80  data=p.var(iv)
81  dv=np.product(p.ds)
82  v=v+dv
83  d=data
84  dmx=float(d.max())
85  dmn=float(d.min())
86  if dmx>dmax:
87  dmax=dmx
88  m=maxloc(d)
89  rmax=[p.x[m[0]],p.y[m[1]],p.z[m[2]]]
90  if dmn < dmin:
91  dmin=dmn
92  m=minloc(d)
93  rmin=[p.x[m[0]],p.y[m[1]],p.z[m[2]]]
94  f=f+dv*np.sum(d)
95  av=float(f/v)
96  f=0.0
97  for p in pp:
98  if iv in p.all_keys:
99  data=p.var(iv)
100  dv=np.product(p.ds)
101  v=v+dv
102  f=f+dv*np.sum((data-av)**2)
103  rm=np.sqrt(float(f/v))
104  print('average: {:12.4e}'.format(av))
105  print(' rms: {:12.4e}'.format(rm))
106  print(' max: {:12.4e}, at {}'.format(dmax,rmax))
107  print(' min: {:12.4e}, at {}'.format(dmin,rmin))
108 
109 def haver1(s,iv=0,point=0.0,dir=2,x=None,y=None,z=None,all=False,i4=0,verbose=0):
110  """ Horizontal average in one plane """
111  pp=_patches(s)
112  n_patch=0
113  f=0.0
114  v=0.0
115  if x:
116  dir=0
117  point=x
118  elif y:
119  dir=1
120  point=y
121  elif z:
122  dir=2
123  point=z
124  for p in pp:
125  ds=p.ds[dir]
126  lc=p.position[dir]-p.size[dir]/2.0
127  uc=p.position[dir]+p.size[dir]/2.0
128  if all and p.guard_zones:
129  lc=lc-p.ng[dir]*p.ds[dir]
130  uc=uc+p.ng[dir]*p.ds[dir]
131  n=p.gn[dir]
132  else:
133  n=p.n[dir]
134  if uc >= point and lc <= point:
135  n_patch+=1
136  if p.no_mans_land:
137  lc=lc+p.ds[dir]/2
138  fi=(point-lc)/ds
139  i=max(0,min(n-2,np.floor(fi)))
140  fi=fi-i
141  if fi < 0:
142  w0=1.0+fi
143  w1=None
144  elif fi > 1.0:
145  w0=2.0-fi
146  w1=None
147  else:
148  w0=1.0-fi
149  w1=fi
150  if verbose>1:
151  print ('id:',p.id,'i:',i,'w:',w0,w1)
152  n+=1
153  dv=np.product(p.ds)
154  data=p.var(iv,i4=i4,all=all)
155  v=v+dv*w0
156  if w1:
157  v=v+dv*w1
158  if dir==0:
159  f=f+dv*np.sum(data[i,:,:])*w0
160  if w1:
161  f=f+dv*np.sum(data[i+1,:,:])*w1
162  elif dir==1:
163  f=f+dv*np.sum(data[:,i,:])*w0
164  if w1:
165  f=f+dv*np.sum(data[:,i+1,:])*w1
166  else:
167  f=f+dv*np.sum(data[:,:,i])*w0
168  if w1:
169  f=f+dv*np.sum(data[:,:,i+1])*w1
170  if verbose>0:
171  print('using',n_patch,'patches')
172  return f/v
173 
174 def hminmax(s,iv=0,dir=2,i4=0,all=False,verbose=0):
175  """ Find the smallest and largest values in planes perpendicular to dir """
176  pp=_patches(s)
177  p=pp[0]
178  xx=[]
179  hmin=[]
180  hmax=[]
181  for p in pp:
182  if all:
183  rr=p.xyz[dir]
184  else:
185  rr=p.xyzi[dir]
186  if iv in p.all_keys:
187  data=p.var(iv,all=all,i4=i4)
188  for i in range(data.shape[dir]):
189  if dir==0:
190  f=data[i,:,:]
191  elif dir==1:
192  f=data[:,i,:]
193  else:
194  f=data[:,:,i]
195  xx.append(rr[i])
196  hmin.append(f.min())
197  hmax.append(f.max())
198  # sort the results
199  xx=np.array(xx)
200  hmin=np.array(hmin)
201  hmax=np.array(hmax)
202  w=xx.argsort()
203  xx=xx[w]
204  hmin=hmin[w]
205  hmax=hmax[w]
206  # harvest
207  ss=[]
208  ymin=[]
209  ymax=[]
210  i=0
211  while i<len(xx):
212  x0=xx[i]
213  h1=hmin[i]
214  h2=hmax[i]
215  x0=xx[i]
216  # if there are nore at the same coordinate
217  while i<len(xx) and xx[i]==x0:
218  h1=min(h1,hmin[i])
219  h2=max(h2,hmax[i])
220  i+=1
221  ss.append(x0)
222  ymin.append(h1)
223  ymax.append(h2)
224  xx=np.array(ss)
225  hmin=np.array(ymin)
226  hmax=np.array(ymax)
227  return xx,hmin,hmax
228 
229 def haver(s,iv=0,dir=2,i4=0,all=False,verbose=0):
230  """ Horizontal average in all planes """
231  pp=_patches(s)
232  p=pp[0]
233  xx=[]
234  hav=[]
235  jv=dispatch.map_var(p,iv)
236  for p in pp:
237  es=0.5 if p.no_mans_land else 0.0
238  es=es+p.idx.h[dir,jv]
239  n=p.gn if all else p.n
240  if all:
241  rr=p.xyz[dir]
242  else:
243  rr=p.xyzi[dir]
244  if iv in p.all_keys:
245  data=p.var(iv,i4=i4,all=all)
246  for i in range(data.shape[dir]):
247  if dir==0:
248  f=data[i,:,:]
249  elif dir==1:
250  f=data[:,i,:]
251  else:
252  f=data[:,:,i]
253  xx.append(rr[i])
254  hav.append(f.mean())
255  # sort the results
256  xx=np.array(xx)
257  hav=np.array(hav)
258  w=xx.argsort()
259  xx=xx[w]
260  hav=hav[w]
261  # harvest
262  ss=[]
263  hh=[]
264  i=0
265  while i<len(xx):
266  x0=xx[i]
267  n=0
268  ha=0.0
269  while i<len(xx) and xx[i]==x0:
270  ha+=hav[i]
271  i+=1
272  n+=1
273  ss.append(x0)
274  hh.append(ha/n)
275  xx=np.array(ss)
276  hh=np.array(hh)
277  return xx,hh