# Path: pygeo/fastio.py
# def getParams (filename = 'for.header'):
#   with open(filename, 'r') as fp:
#     res = fp.read()
# 
#   res = res.strip().split()
#   paramdict = {}
#   paramdict['x'] = [float(item) for item in res[0:2]]
#   paramdict['y'] = [float(item) for item in res[2:4]]
#   paramdict['z'] = [float(item) for item in res[4:6]]
#   paramdict['dx'] = float(res[6])
#   paramdict['dims'] = [int(item) for item in res[7:10]]
# 
#   return paramdict
# 
# def readfast (filename, dims=None):
#   '''
#   Read FAST model file to Python ndarray.
#   '''
# 
#   if (dims is None):
#     dims = getParams()['dims']
# 
#   (nx, ny, nz) = dims
# 
#   fp = open(filename, 'rb')
#   map = mmap.mmap(fp.fileno(), 0, flags=mmap.MAP_PRIVATE)
# 
#   vol = np.zeros(dims, dtype=np.int16)
#   for k in range(nz):
#     for j in range(ny):
#       bytestart = k*ny*(2*rlpad+2*nx) + j*(2*rlpad+2*nx)
#       entries = struct.unpack('%dH'%(nx,), map[bytestart+rlpad:bytestart+rlpad+nx*2])
#       vol[:,j,k] = np.array(entries)
# 
#   map.close()
#   fp.close()
# 
#   return vol
# 
# def writefast (filename, model):
#   '''
#   Write FAST model file from Python ndarray.
#   '''
# 
#   (nx, ny, nz) = model.shape
#   rmark = struct.pack(rlpak, np.int64(nx*2))
# 
#   fp = open(filename, 'wb')
# 
#   for k in range(nz):
#     for j in range(ny):
#       outstring = struct.pack('%dH'%(nx,), *list(model[:,j,k]))
#       fp.write(rmark)
#       fp.write(outstring)
#       fp.write(rmark)
# 
#   fp.close()


import os
import warnings
import numpy as np
from optparse import OptionParser
from pygeo.fastio import getParams, readfast, writefast
from scipy.interpolate import griddata
# Take the weighted average of the cells in each x,z location along the
# y-axis, normalized by the total number of rays in that x,z location
avg2d = np.nan_to_num(((1. * hitcount * model3d) / hitcount3dsum).sum(axis=1))

# Create mask (value of 1 where there are no rays)
selector = hitcount2dsum > 0
baseblank = 1 - selector

if (options.extrap):
  nx, nz = avg2d.shape
  X, Z = np.mgrid[0:nx, 0:nz]
  P = np.array([X.ravel(),Z.ravel()]).T

  x0 = X[selector]
  z0 = Z[selector]
  pos0 = np.array([x0.ravel(),z0.ravel()]).T
  val = avg2d[selector]
  
  model2d = griddata(pos0, val, P, method='nearest').reshape((nx,nz))

else:
  # Create an empty (2D) model in which to store the result
  # Copy the 2D average model into the masked region
  model2d = avg2d.copy() + model2davg.copy() * baseblank

if (options.direct2d):
  writefast(outfile, model2d.reshape(pseudodims))
else:
  newmodel3d = model3d.copy()
  newmodel3d[:] = model2d.reshape(pseudodims)

# Path: pygeo/fastio.py
# def getParams (filename = 'for.header'):
#   with open(filename, 'r') as fp:
#     res = fp.read()
# 
#   res = res.strip().split()
#   paramdict = {}
#   paramdict['x'] = [float(item) for item in res[0:2]]
#   paramdict['y'] = [float(item) for item in res[2:4]]
#   paramdict['z'] = [float(item) for item in res[4:6]]
#   paramdict['dx'] = float(res[6])
#   paramdict['dims'] = [int(item) for item in res[7:10]]
# 
#   return paramdict
# 
# def readfast (filename, dims=None):
#   '''
#   Read FAST model file to Python ndarray.
#   '''
# 
#   if (dims is None):
#     dims = getParams()['dims']
# 
#   (nx, ny, nz) = dims
# 
#   fp = open(filename, 'rb')
#   map = mmap.mmap(fp.fileno(), 0, flags=mmap.MAP_PRIVATE)
# 
#   vol = np.zeros(dims, dtype=np.int16)
#   for k in range(nz):
#     for j in range(ny):
#       bytestart = k*ny*(2*rlpad+2*nx) + j*(2*rlpad+2*nx)
#       entries = struct.unpack('%dH'%(nx,), map[bytestart+rlpad:bytestart+rlpad+nx*2])
#       vol[:,j,k] = np.array(entries)
# 
#   map.close()
#   fp.close()
# 
#   return vol
# 
# def writefast (filename, model):
#   '''
#   Write FAST model file from Python ndarray.
#   '''
# 
#   (nx, ny, nz) = model.shape
#   rmark = struct.pack(rlpak, np.int64(nx*2))
# 
#   fp = open(filename, 'wb')
# 
#   for k in range(nz):
#     for j in range(ny):
#       outstring = struct.pack('%dH'%(nx,), *list(model[:,j,k]))
#       fp.write(rmark)
#       fp.write(outstring)
#       fp.write(rmark)
# 
#   fp.close()


# Path: util/slicefromhitcount.py
import os
import warnings
import numpy as np
from optparse import OptionParser
from pygeo.fastio import getParams, readfast, writefast
from scipy.interpolate import griddata

# Take the weighted average of the cells in each x,z location along the
# y-axis, normalized by the total number of rays in that x,z location
avg2d = np.nan_to_num(((1. * hitcount * model3d) / hitcount3dsum).sum(axis=1))

# Create mask (value of 1 where there are no rays)
selector = hitcount2dsum > 0
baseblank = 1 - selector

if (options.extrap):
  nx, nz = avg2d.shape
  X, Z = np.mgrid[0:nx, 0:nz]
  P = np.array([X.ravel(),Z.ravel()]).T

  x0 = X[selector]
  z0 = Z[selector]
  pos0 = np.array([x0.ravel(),z0.ravel()]).T
  val = avg2d[selector]
  
  model2d = griddata(pos0, val, P, method='nearest').reshape((nx,nz))

else:
  # Create an empty (2D) model in which to store the result
  # Copy the 2D average model into the masked region
  model2d = avg2d.copy() + model2davg.copy() * baseblank

if (options.direct2d):
  writefast(outfile, model2d.reshape(pseudodims))
else:
  newmodel3d = model3d.copy()
  newmodel3d[:] = model2d.reshape(pseudodims)

  writefast(outfile, newmodel3d)