def fit(t,x,d,sd,
        network_model,
        network_params,
        station_model,
        station_params):
  ''' 
  Fit network and station processes to the observations, not
  distinguishing between signal and noise.
  '''
  t = np.asarray(t,dtype=float)
  x = np.asarray(x,dtype=float)
  d = np.array(d,dtype=float)
  sd = np.array(sd,dtype=float)
  diff = np.array([0,0,0])

  net_gp = composite(network_model,network_params,gpnetwork.CONSTRUCTORS)
  sta_gp = composite(station_model,station_params,gpstation.CONSTRUCTORS)

  t_grid,x0_grid = np.meshgrid(t,x[:,0],indexing='ij')
  t_grid,x1_grid = np.meshgrid(t,x[:,1],indexing='ij')
  # flat observation times and positions
  z = np.array([t_grid.ravel(),
                x0_grid.ravel(),
                x1_grid.ravel()]).T

  # mask indicates missing data
  mask = np.isinf(sd)
  z,d,sd = z[~mask.ravel()],d[~mask],sd[~mask]

  # Build covariance and basis vectors for the combined process. Do
  # not evaluated at masked points
  sta_sigma,sta_p = station_sigma_and_p(sta_gp,t,mask)
  net_sigma = net_gp._covariance(z,z,diff,diff)
  net_p = net_gp._basis(z,diff)
  # combine station gp with the network gp
  mu = np.zeros(z.shape[0])
  sigma = as_sparse_or_array(sta_sigma + net_sigma)
  p = np.hstack((sta_p,net_p))
  del sta_sigma,net_sigma,sta_p,net_p
  # best fit combination of signal and noise to the observations
  uf,suf = _fit(d,sd,mu,sigma,p)
  # fold back into 2d arrays
  u = np.full((t.shape[0],x.shape[0]),np.nan)
  u[~mask] = uf
  su = np.full((t.shape[0],x.shape[0]),np.inf)
  su[~mask] = suf
  return u,su def reml(t,x,d,sd,
         network_model,
         network_params,
         network_fix,
         station_model,
         station_params,
         station_fix):
  ''' 
  Returns the Restricted Maximum Likelihood (REML) estimatates of the
  unknown hyperparameters.
  '''
  t = np.asarray(t,dtype=float)
  x = np.asarray(x,dtype=float)
  d = np.array(d,dtype=float)
  sd = np.array(sd,dtype=float)
  diff = np.array([0,0,0])

  t_grid,x0_grid = np.meshgrid(t,x[:,0],indexing='ij')
  t_grid,x1_grid = np.meshgrid(t,x[:,1],indexing='ij')
  # flat observation times and positions
  z = np.array([t_grid.ravel(),
                x0_grid.ravel(),
                x1_grid.ravel()]).T

  # mask indicates missing data
  mask = np.isinf(sd)
  z,d,sd = z[~mask.ravel()],d[~mask],sd[~mask]
  # number of network hyperparameters
  n = len(network_params)
  # combined network and station parameters
  params = np.hstack((network_params,station_params))
  fix = np.hstack((network_fix,station_fix+n))
  free = np.array([i for i in range(len(params)) if i not in fix],dtype=int)
  
  def objective(theta):
    logger.debug('Current hyperparameters : ' + ' '.join('%0.4e' % i for i in theta))
    test_params = np.copy(params)
    test_params[free] = theta 
    test_network_params = test_params[:n]
    test_station_params = test_params[n:]
    net_gp = composite(network_model,test_network_params,gpnetwork.CONSTRUCTORS)
    sta_gp = composite(station_model,test_station_params,gpstation.CONSTRUCTORS)
    # station process
    sta_sigma,sta_p = station_sigma_and_p(sta_gp,t,mask)
    # add data noise to the diagonals of sta_sigma. Both matrices are
    # sparse so this is efficient
    obs_sigma = _as_covariance(sd)
    sta_sigma = as_sparse_or_array(sta_sigma + obs_sigma)
    # network process
    net_sigma = net_gp._covariance(z,z,diff,diff)
    net_p = net_gp._basis(z,diff)
    # combine station gp with the network gp
    mu = np.zeros(z.shape[0])
    sigma = as_sparse_or_array(sta_sigma + net_sigma)
    p = np.hstack((sta_p,net_p))
    del sta_sigma,net_sigma,obs_sigma,sta_p,net_p
    try:
      out = likelihood(d,mu,sigma,p=p)
    except np.linalg.LinAlgError as err:
      logger.warning(
        'An error was raised while computing the log '
        'likelihood:\n\n%s\n' % repr(err))
      logger.warning('Returning -INF for the log likelihood')   
      out = -np.inf
      
    logger.debug('Log likelihood : %.8e' % out)
    return out  

  opt,val = fmax_pos(objective,params[free],disp=False)
  logger.info('Optimal hyperparameters : ' + ' '.join('%.4e' % i for i in opt))
  params[free] = opt
  out_network_params = params[:n]
  out_station_params = params[n:]
  out_likelihood = val
  return out_network_params,out_station_params,out_likelihood def strain(t,x,d,sd,
           network_prior_model,
           network_prior_params,
           network_noise_model,
           network_noise_params,
           station_noise_model,
           station_noise_params,
           out_t,out_x,rate,
           covariance):
  ''' 
  Computes deformation gradients from displacement data.
  '''  
  t = np.asarray(t,dtype=float)
  x = np.asarray(x,dtype=float)
  d = np.array(d,dtype=float)
  sd = np.array(sd,dtype=float)
  diff = np.array([0,0,0])

  t_grid,x0_grid = np.meshgrid(t,x[:,0],indexing='ij')  
  t_grid,x1_grid = np.meshgrid(t,x[:,1],indexing='ij')  
  # flat observation times and positions
  z = np.array([t_grid.ravel(),
                x0_grid.ravel(),
                x1_grid.ravel()]).T

  t_grid,x0_grid = np.meshgrid(out_t,out_x[:,0],indexing='ij')  
  t_grid,x1_grid = np.meshgrid(out_t,out_x[:,1],indexing='ij')  
  # flat observation times and positions
  out_z = np.array([t_grid.ravel(),
                    x0_grid.ravel(),
                    x1_grid.ravel()]).T

  prior_gp = composite(network_prior_model,network_prior_params,gpnetwork.CONSTRUCTORS)
  noise_gp = composite(network_noise_model,network_noise_params,gpnetwork.CONSTRUCTORS)
  sta_gp   = composite(station_noise_model,station_noise_params,gpstation.CONSTRUCTORS)

  # find missing data
  mask = np.isinf(sd)
  # get unmasked data and uncertainties
  z,d,sd = z[~mask.ravel()],d[~mask],sd[~mask]
  # build noise covariance and basis vectors
  sta_sigma,sta_p = station_sigma_and_p(sta_gp,t,mask)
  # add data noise to the station noise
  obs_sigma = _as_covariance(sd)
  sta_sigma = as_sparse_or_array(sta_sigma + obs_sigma)
  # make network noise
  net_sigma = noise_gp._covariance(z,z,diff,diff)
  net_p = noise_gp._basis(z,diff)
  # combine noise processes
  noise_sigma = as_sparse_or_array(sta_sigma + net_sigma)
  noise_p = np.hstack((sta_p,net_p))
  del sta_sigma,net_sigma,obs_sigma,sta_p,net_p
  # condition the prior with the data
  post_gp = prior_gp.condition(z,d,sigma=noise_sigma,p=noise_p)
  if rate:
    dudx_gp = post_gp.differentiate((1,1,0)) # x derivative of velocity
    dudy_gp = post_gp.differentiate((1,0,1)) # y derivative of velocity

  else:  
    dudx_gp = post_gp.differentiate((0,1,0)) # x derivative of displacement
    dudy_gp = post_gp.differentiate((0,0,1)) # y derivative of displacement

  if covariance:
    # Evaluate the mean and covariances of the posterior
    dudx = dudx_gp.mean(out_z)
    cdudx = dudx_gp.covariance(out_z,out_z)
    sdudx = np.sqrt(np.diag(cdudx))

    dudy = dudy_gp.mean(out_z)
    cdudy = dudy_gp.covariance(out_z,out_z)
    sdudy = np.sqrt(np.diag(cdudy))

    dudx  = dudx.reshape((out_t.shape[0],out_x.shape[0]))
    sdudx = sdudx.reshape((out_t.shape[0],out_x.shape[0]))
    cdudx = cdudx.reshape((out_t.shape[0],out_x.shape[0],
                           out_t.shape[0],out_x.shape[0]))
    dudy  = dudy.reshape((out_t.shape[0],out_x.shape[0]))
    sdudy = sdudy.reshape((out_t.shape[0],out_x.shape[0]))
    cdudy = cdudy.reshape((out_t.shape[0],out_x.shape[0],
                           out_t.shape[0],out_x.shape[0]))
    out = (dudx,sdudx,cdudx,dudy,sdudy,cdudy)

  else:
    # Just evaluate the mean and standard deviations of the posterior
    dudx,sdudx = dudx_gp.meansd(out_z,chunk_size=1000)
    dudy,sdudy = dudy_gp.meansd(out_z,chunk_size=1000)

    dudx = dudx.reshape((out_t.shape[0],out_x.shape[0]))
    sdudx = sdudx.reshape((out_t.shape[0],out_x.shape[0]))
    dudy = dudy.reshape((out_t.shape[0],out_x.shape[0]))
    sdudy = sdudy.reshape((out_t.shape[0],out_x.shape[0]))
    out = (dudx,sdudx,dudy,sdudy)
        
  return out def autoclean(t,x,d,sd,
              network_model,
              network_params,
              station_model,
              station_params,
              tol):
  ''' 
  Returns a dataset that has been cleaned of outliers using a data
  editing algorithm.
  '''
  t = np.asarray(t,dtype=float)
  x = np.asarray(x,dtype=float)
  de = np.array(d,dtype=float,copy=True)
  sde = np.array(sd,dtype=float,copy=True)
  diff = np.array([0,0,0])

  net_gp = composite(network_model,network_params,gpnetwork.CONSTRUCTORS)
  sta_gp = composite(station_model,station_params,gpstation.CONSTRUCTORS)

  t_grid,x0_grid = np.meshgrid(t,x[:,0],indexing='ij')
  t_grid,x1_grid = np.meshgrid(t,x[:,1],indexing='ij')
  # flat observation times and positions
  z = np.array([t_grid.ravel(),
                x0_grid.ravel(),
                x1_grid.ravel()]).T

  # mask indicates missing data
  mask = np.isinf(sde)
  zu,du,sdu = z[~mask.ravel()],de[~mask],sde[~mask]
  # Build covariance and basis vectors for the combined process. Do
  # not evaluated at masked points
  sta_sigma,sta_p = station_sigma_and_p(sta_gp,t,mask)
  net_sigma = net_gp._covariance(zu,zu,diff,diff)
  net_p = net_gp._basis(zu,diff)
  # combine station gp with the network gp
  mu = np.zeros(zu.shape[0])  
  sigma = as_sparse_or_array(sta_sigma + net_sigma)
  p = np.hstack((sta_p,net_p))
  del sta_sigma,net_sigma,sta_p,net_p
  # returns the indices of outliers 
  out_idx = outliers(du,sdu,
                     mu=mu,sigma=sigma,p=p,
                     tol=tol)
  # mask the outliers in *de* and *sde*
  r,c = np.nonzero(~mask)
  de[r[out_idx],c[out_idx]] = np.nan
  sde[r[out_idx],c[out_idx]] = np.inf
  return (de,sde) def composite_units(components,constructors):
  ''' 
  returns the units for the composite Gaussian process hyperparameters
  '''
  components = list(components)
  try:
    cs = [constructors[m] for m in components]
  except KeyError as err:
    raise ValueError(
      '"%s" is not a valid Gaussian process. Use Gaussian processes '
      'from the following list:\n%s' % 
      (err.args[0],', '.join(['"%s"' % i for i in constructors.keys()])))

  units = []
  for ci in cs:
    units += ci.units
  
  return units   def se(sigma,cls):
def exp(sigma,cls):
def mat32(sigma,cls):
def mat52(sigma,cls):
def wen32(sigma,cls):
def spwen32(sigma,cls):
def se_se(sigma,cts,cls):
def exp_se(sigma,cts,cls):
def wen11_se(sigma,cts,cls):
def wen12_se(sigma,cts,cls):
def spwen11_se(sigma,cts,cls):
def spwen12_se(sigma,cts,cls):
def fogm_se(sigma,fc,cls):
def bm_se(sigma,t0,cls):
def ibm_se(sigma,t0,cls):
def mat32_se(sigma,cts,cls):
def mat52_se(sigma,cts,cls):
CONSTRUCTORS = {'bm-se':bm_se,
                'ibm-se':ibm_se,
                'fogm-se':fogm_se,
                'se-se':se_se,
                'exp-se':exp_se,
                'mat32-se':mat32_se,
                'mat52-se':mat52_se,
                'wen11-se':wen11_se,
                'wen12-se':wen12_se,
                'spwen11-se':spwen11_se,
                'spwen12-se':spwen12_se} def const():
  def basis(x,diff):
def linear():
  def basis(x,diff):
def per():
  def basis(x):
def step(t0):
  def basis(x,diff):
def mat32(sigma,cts):
def mat52(sigma,cts):
def wen11(sigma,cts):
def wen12(sigma,cts):
def wen30(sigma,cts):
def spwen11(sigma,cts):
def spwen12(sigma,cts):
def spwen30(sigma,cts):
def se(sigma,cts):
def exp(sigma,cts):
def fogm(sigma,w):
def bm(sigma,t0):
  def mean(x):
  def cov(x1,x2):
def ibm(sigma,t0):
  def mean(x,diff):
  def cov(x1,x2,diff1,diff2):
CONSTRUCTORS = {'const':const,
                'linear':linear,
                'per':per, 
                'step':step,
                'bm':bm,
                'ibm':ibm,
                'fogm':fogm,
                'mat32':mat32,
                'mat52':mat52,
                'wen11':wen11,
                'wen12':wen12,
                'wen30':wen30,
                'spwen11':spwen11,
                'spwen12':spwen12,
                'spwen30':spwen30,
                'se':se,
                'exp':exp} @_memoize
def mjd_inv(m,fmt):
  ''' 
  Converts Modified Julian Date (MJD) to a date string 
  
  Parameters
  ----------
  m : int
    Modified Julian Date

  fmt : string
    format string indicating how to form *out*
      
  Returns
  -------
  out : str
    Date string
    
  '''
  m = int(m)
  d = _REFERENCE_DATETIME + timedelta(m)
  out = d.strftime(fmt)
  return out @_memoize
def mjd(s,fmt):
  ''' 
  Converts a date string into Modified Julian Date (MJD)
  
  Parameters
  ----------
  s : string
    Date string

  fmt : string
    Format string indicating how to parse *s*
      
  Returns
  -------
  out : int
    Modified Julian Date

  '''
  d = datetime.strptime(s,fmt)
  out = (d - _REFERENCE_DATETIME).days
  return out def make_basemap(lon,lat,resolution=None):
  ''' 
  Creates a transverse mercator projection which is centered about the 
  given positions. 

  Modifying this function should be sufficient to change the way that 
  all map projections are generated in PyGeoNS
  '''
  lon = np.asarray(lon)
  lat = np.asarray(lat)
  n = lon.shape[0]
  # the buffer should be the average shortest distance between 
  # points... but this is good enough
  lon_buff = max(0.1,lon.ptp()/np.sqrt(n))
  lat_buff = max(0.1,lat.ptp()/np.sqrt(n))
  llcrnrlon = min(lon) - lon_buff
  llcrnrlat = min(lat) - lat_buff
  urcrnrlon = max(lon) + lon_buff
  urcrnrlat = max(lat) + lat_buff
  lon_0 = (llcrnrlon + urcrnrlon)/2.0
  lat_0 = (llcrnrlat + urcrnrlat)/2.0
  return Basemap(projection='tmerc',
                 resolution=resolution,
                 lon_0 = lon_0,
                 lat_0 = lat_0,
                 llcrnrlon = llcrnrlon,
                 llcrnrlat = llcrnrlat,
                 urcrnrlon = urcrnrlon,
                 urcrnrlat = urcrnrlat) def dict_from_hdf5(infile):
  ''' 
  Loads a data dictionary from an hdf5 file.
  
  Parameters
  ----------
  infile : str
    Name of hdf5 file
    
  Returns
  -------
  out : dict
    Data dictionary

  '''
  out = {}
  fin = h5py.File(infile,'r')
  for k in fin.keys():
    out[k] = fin[k][...]

  fin.close()
  check_data(out)
  return out   def hdf5_from_dict(outfile,data):
  ''' 
  Writes an hdf5 file from the data dictionary.
  
  Parameters
  ----------
  outfile : str
    Name of the output file
  
  data : dict
    Data dictionary      
  
  '''
  check_data(data)
  fout = h5py.File(outfile,'w') 
  for k in data.keys():
    fout[k] = data[k]
    
  fout.close()
  return
import numpy as np
import logging
import subprocess as sp
from pygeons.main.fit import fit
from pygeons.main.reml import reml
from pygeons.main.strain import strain
from pygeons.main.autoclean import autoclean
from pygeons.main.gptools import composite_units
from pygeons.main import gpnetwork
from pygeons.main import gpstation
from pygeons.mjd import mjd_inv,mjd
from pygeons.basemap import make_basemap
from pygeons.io.convert import dict_from_hdf5,hdf5_from_dict
  
  # convert geodetic positions to cartesian
  bm = make_basemap(data['longitude'],data['latitude'])
  x,y = bm(data['longitude'],data['latitude'])
  xy = np.array([x,y]).T

  _log_fit(input_file,
           network_model,network_params,
           station_model,station_params,
           output_file)
  
  for dir in ['east','north','vertical']:
    u,su = fit(t=data['time'][:,None],
               x=xy,      
               d=data[dir],
               sd=data[dir+'_std_dev'],
               network_model=network_model,
               network_params=network_params[dir],
               station_model=station_model,
               station_params=station_params[dir])
    out[dir] = u
    out[dir+'_std_dev'] = su

  hdf5_from_dict(output_file,out)
  logger.info('Posterior fit written to %s' % output_file)
  return


def pygeons_autoclean(input_file,
                      network_model=('spwen12-se',),
                      network_params=(1.0,0.1,100.0),
                      station_model=('linear',),
                      station_params=(),
                      output_stem=None,
                      outlier_tol=4.0):
  ''' 
  Remove outliers with a data editing algorithm
  '''
  logger.info('Running pygeons autoclean ...')
  data = dict_from_hdf5(input_file)
  if data['time_exponent'] != 0:
    raise ValueError('input dataset must have units of displacement')

  if data['space_exponent'] != 1:
    raise ValueError('input dataset must have units of displacement')

  # dictionary which will contain the edited data
  out = dict((k,np.copy(v)) for k,v in data.iteritems())

  # convert params to a dictionary of hyperparameters for each direction
  network_params = _params_dict(network_params)
  station_params = _params_dict(station_params)

  # make output file name
  if output_stem is None:
    output_stem = _remove_extension(input_file) + '.autoclean'

  output_file = output_stem + '.h5'
  
  # convert geodetic positions to cartesian
  bm = make_basemap(data['longitude'],data['latitude'])
  x,y = bm(data['longitude'],data['latitude'])
  xy = np.array([x,y]).T

  _log_autoclean(input_file,
                 network_model,network_params,
                 station_model,station_params,
                 outlier_tol,
                 output_file)
  
  for dir in ['east','north','vertical']:
    de,sde = autoclean(t=data['time'][:,None],
                       x=xy, 
                       d=data[dir],
                       sd=data[dir+'_std_dev'],
                       network_model=network_model,
                       network_params=network_params[dir],
                       station_model=station_model,
                       station_params=station_params[dir],
                       tol=outlier_tol)
    out[dir] = de
    out[dir+'_std_dev'] = sde

  hdf5_from_dict(output_file,out)
  logger.info('Edited data written to %s' % output_file)
  return


def pygeons_reml(input_file,
                 network_model=('spwen12-se',),
                 network_params=(1.0,0.1,100.0),
                 network_fix=(),
                 station_model=('linear',),
                 station_params=(),
                 station_fix=(),
                 output_stem=None):
  ''' 
  Restricted maximum likelihood estimation
  '''
  logger.info('Running pygeons reml ...')
  data = dict_from_hdf5(input_file)
  if data['time_exponent'] != 0:
    raise ValueError('input dataset must have units of displacement')

  if data['space_exponent'] != 1:
    raise ValueError('input dataset must have units of displacement')

  # convert params to a dictionary of hyperparameters for each direction
  network_params = _params_dict(network_params)
  network_fix = np.asarray(network_fix,dtype=int)
  station_params = _params_dict(station_params)
  station_fix = np.asarray(station_fix,dtype=int)

  # make output file name
  if output_stem is None:
    output_stem = _remove_extension(input_file) + '.reml'

  output_file = output_stem + '.txt'
  
  # convert geodetic positions to cartesian

  bm = make_basemap(data['longitude'],data['latitude'])