<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Demo 2 - Remote parallel computation [distributed]
Demo for site visit | Brendan Smithyman | April 8, 2015
Choice of IPython / jupyter cluster profile
Step1:  Importing libraries
numpy is the de facto standard Python numerical computing library
zephyr.Dispatcher is zephyr's primary parallel remote problem interface
IPython.parallel provides parallel task control (nominally, this is to be handled inside the Dispatcher object)
Step2:  Plotting configuration
These lines import matplotlib, which is a standard Python plotting library, and configure the output formats for figures.
Step3:  These lines define some plotting functions, which are used later.
Step4:  System / modelling configuration
This code sets up the seismic problem; see the comments inline. In a live inversion problem this would most likely be read from a configuration file (but could be defined interactively for development purposes).
Properties of the grid and forward modelling
Step5:  Properties of the model
Step6:  Array geometry
Step7:  Numerical / parallel parameters
Step8:  Computed properties
Step9:  Problem geometry
Step10:  Configuration dictionary
(assembled from previous sections)
Step11:  Parallel computations
This section runs each of the parallel computations on the remote worker nodes.
Set up problem
Create the Dispatcher object using the systemConfig dictionary as input
Spawn survey and problem interfaces, which implement the SimPEG standard properties
Generate a set of "transmitter" objects, each of which knows about its respective "receivers" (in seismic parlance, these would be "sources" and "receivers"; the term "transmitter" is more common in EM and potential fields geophysics)
Tell the dispatcher object about the transmitters
Step12:  Forward modelling and backpropagation
Example (commented out) showing how to generate synthetic data using the SimPEG-style survey and problem interfaces. In this implementation, both are essentially expressions of the Dispatcher. The Dispatcher API has yet to be merged into SimPEG.
Step13:  This code runs the forward modelling on the [remote] workers. It returns asynchronously, so the code can run in the background.
Step14:  However, it will block if we ask for the data or wavefields
Step15:  Results
We show the resulting data and wavefield properties.
Data selection
Step16:  Geometry, data and forward wavefield
<ASSISTANT_TASK:>
Python Code:
# profile = 'phobos'   # remote workstation
# profile = 'pantheon' # remote cluster
profile = 'mpi' # local machine
Explanation: Demo 2 - Remote parallel computation [distributed]
Demo for site visit | Brendan Smithyman | April 8, 2015
Choice of IPython / jupyter cluster profile
End of explanation
import numpy as np
from zephyr.Dispatcher import SeisFDFDDispatcher
from IPython.parallel import Reference
Explanation: Importing libraries
numpy is the de facto standard Python numerical computing library
zephyr.Dispatcher is zephyr's primary parallel remote problem interface
IPython.parallel provides parallel task control (nominally, this is to be handled inside the Dispatcher object)
End of explanation
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
%matplotlib inline
import mpld3
mpld3.enable_notebook()
Explanation: Plotting configuration
These lines import matplotlib, which is a standard Python plotting library, and configure the output formats for figures.
End of explanation
lclip = 2000
hclip = 3000
clipscale = 0.1
sms = 0.5
rms = 0.5
def plotField(u):
    clip = clipscale*abs(u).max()
    plt.imshow(u.real, cmap=cm.bwr, vmin=-clip, vmax=clip)
def plotModel(v):
    plt.imshow(v.real, cmap=cm.jet, vmin=lclip, vmax=hclip)
def plotGeometry(geom):
    
    srcpos = geom['src'][:,::2]
    recpos = geom['rec'][:,::2]
    
    axistemp = plt.axis()
    plt.plot(srcpos[:,0], srcpos[:,1], 'kx', markersize=sms)
    plt.plot(recpos[:,0], recpos[:,1], 'kv', markersize=rms)
    plt.axis(axistemp)
Explanation: These lines define some plotting functions, which are used later.
End of explanation
cellSize    = 1             # m
nx          = 164           # count
nz          = 264           # count
freqs       = [1e2]         # Hz
freeSurf    = [False, False, False, False] # t r b l
nPML        = 32            # number of PML points
nky         = 80            # number of y-directional plane-wave components
Explanation: System / modelling configuration
This code sets up the seismic problem; see the comments inline. In a live inversion problem this would most likely be read from a configuration file (but could be defined interactively for development purposes).
Properties of the grid and forward modelling
End of explanation
velocity    = 2500          # m/s
vanom       = 500           # m/s
density     = 2700          # units of density
Q           = 500           # can be inf
Explanation: Properties of the model
End of explanation
srcs        = np.array([np.ones(101)*32, np.zeros(101), np.linspace(32, 232, 101)]).T
recs        = np.array([np.ones(101)*132, np.zeros(101), np.linspace(32, 232, 101)]).T
nsrc        = len(srcs)
nrec        = len(recs)
recmode     = 'fixed'
geom        = {
    'src':  srcs,
    'rec':  recs,
    'mode': 'fixed',
}
Explanation: Array geometry
End of explanation
cache       = False         # whether to cache computed wavefields for a given source
cacheDir    = '.'
parFac = 2
chunksPerWorker = 0.5       # NB: parFac * chunksPerWorker = number of source array subsets
ensembleClear = False
Explanation: Numerical / parallel parameters
End of explanation
dims        = (nx,nz)       # tuple
rho         = np.fliplr(np.ones(dims) * density)
nfreq       = len(freqs)    # number of frequencies
nsp         = nfreq * nky   # total number of 2D subproblems
cPert       = np.zeros(dims)
cPert[(nx/2)-20:(nx/2)+20,(nz/2)-20:(nz/2)+20] = vanom
c           = np.fliplr(np.ones(dims) * velocity)
cFlat       = c
c          += np.fliplr(cPert)
cTrue       = c
Explanation: Computed properties
End of explanation
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
plotModel(c.T)
plotGeometry(geom)
ax1.set_title('Velocity Model')
ax1.set_xlabel('X')
ax1.set_ylabel('Z')
fig.tight_layout()
Explanation: Problem geometry
End of explanation
# Base configuration for all subproblems
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'c':        c.T,        # m/s
    'rho':      rho.T,      # density
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
    'geom':     geom,
    'cache':    cache,
    'cacheDir': cacheDir,
    'freqs':    freqs,
    'nky':      nky,
    'parFac':   parFac,
    'chunksPerWorker':  chunksPerWorker,
    'profile':  profile,
    'ensembleClear':    ensembleClear,
#     'MPI': False,
#    'Solver':   Reference('SimPEG.SolverWrapD(scipy.sparse.linalg.splu)'),#Solver,
}
Explanation: Configuration dictionary
(assembled from previous sections)
End of explanation
%%time
sp = SeisFDFDDispatcher(systemConfig)
survey, problem = sp.spawnInterfaces()
sxs = survey.genSrc()
sp.srcs = sxs
Explanation: Parallel computations
This section runs each of the parallel computations on the remote worker nodes.
Set up problem
Create the Dispatcher object using the systemConfig dictionary as input
Spawn survey and problem interfaces, which implement the SimPEG standard properties
Generate a set of "transmitter" objects, each of which knows about its respective "receivers" (in seismic parlance, these would be "sources" and "receivers"; the term "transmitter" is more common in EM and potential fields geophysics)
Tell the dispatcher object about the transmitters
End of explanation
# d = survey.projectFields()
# uF = problem.fields()
Explanation: Forward modelling and backpropagation
Example (commented out) showing how to generate synthetic data using the SimPEG-style survey and problem interfaces. In this implementation, both are essentially expressions of the Dispatcher. The Dispatcher API has yet to be merged into SimPEG.
End of explanation
%%time
sp.forward()
sp.forwardGraph
Explanation: This code runs the forward modelling on the [remote] workers. It returns asynchronously, so the code can run in the background.
End of explanation
%%time
d = sp.dPred
uF = sp.uF
d[0].shape
Explanation: However, it will block if we ask for the data or wavefields:
End of explanation
freqNum = 0
srcNum = 0
frt = uF[freqNum]
drt = d[freqNum]
clipScaleF = 1e-1 * abs(frt[srcNum]).max()
Explanation: Results
We show the resulting data and wavefield properties.
Data selection:
End of explanation
fig = plt.figure()
ax1 = fig.add_subplot(1,3,1)
plotModel(c.T)
plotGeometry(geom)
ax1.set_title('Velocity Model')
ax1.set_xlabel('X')
ax1.set_ylabel('Z')
ax2 = fig.add_subplot(1,3,2)
plt.imshow(drt.real, cmap=cm.bwr)
ax2.set_title('Real part of d: $\omega = %0.1f$'%(freqs[freqNum],))
ax2.set_xlabel('Receiver #')
ax2.set_ylabel('Source #')
ax3 = fig.add_subplot(1,3,3)
plt.imshow(frt[srcNum].real, vmin=-clipScaleF, vmax=clipScaleF, cmap=cm.bwr)
plt.title('uF: $\omega = %0.1f$, src. %d'%(freqs[freqNum], srcNum))
fig.tight_layout()
Explanation: Geometry, data and forward wavefield:
End of explanation
<END_TASK>