In [1]:

    

import sys
sys.path.append('../')
sys.path.append('/usr/local/bin')

import numpy as np
import zephyr.backend
import zephyr.middleware
import SimPEG
import scipy.sparse.linalg


    

In [2]:

    

nx = 100
nz = 200
cBackground = 2500.
cAnomaly = -500.

cInitial = cBackground * np.ones((nz,nx))

cUpdate = np.zeros((nz,nx))
cUpdate[(nz/2)-10:(nz/2)+10,(nx/2)-10:(nx/2)+10] = cAnomaly

cTrue = cInitial + cUpdate

sz = np.arange(25, nz-24, 1)
sx = 25. * np.ones((len(sz),))
rz = np.arange(25, nz-24, 1)
rx = (nx - 25.) * np.ones((len(rz),))

geom = {
    'src':      np.vstack([sx, sz]).T,
    'rec':      np.vstack([rx, rz]).T,
}

systemConfig = {
    'dx':       1.,                             # m
    'dz':       1.,                             # m
    'c':        cTrue.ravel(),                  # m/s
    'rho':      1.,                             # kg/m^3
    'nx':       nx,                             # count
    'nz':       nz,                             # count
    'freqs':    np.arange(50, 450, 50),         # Hz
    'Disc':     zephyr.backend.MiniZephyrHD,    # discretization
    'geom':     geom,                           # dictionary
    'nWorkers': 8,
    'targetGPW':8,
    'cMin':     2000.,
    'Solver':   scipy.sparse.linalg.splu,
}


    

In [3]:

    

problem = zephyr.middleware.Helm2DViscoProblem(systemConfig)
survey  = zephyr.middleware.Helm2DSurvey(systemConfig)
problem.pair(survey)


    

In [4]:

    

%%time
dObs = survey.dpred()


    

    




CPU times: user 919 ms, sys: 688 ms, total: 1.61 s
Wall time: 4.21 s

In [5]:

    

survey.dobs = dObs
survey.std = 1.

dmisfit = SimPEG.DataMisfit.l2_DataMisfit(survey)
# dmisfit.Wd = 1.

opt = SimPEG.Optimization.ProjectedGradient(maxIter=5)
mapping = SimPEG.Maps.IdentityMap(nP=nz*nx)
reg = zephyr.middleware.HelmBaseRegularization(problem.mesh, mapping=mapping)

invProb = SimPEG.InvProblem.BaseInvProblem(dmisfit, reg, opt)
inv = SimPEG.Inversion.BaseInversion(invProb)


    

In [6]:

    

mOpt = inv.run(cInitial.ravel())


    

    




SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same Solver and solverOpts as the problem***
========================================= Projected Gradient =========================================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS   itType   aSet    bSet    Comment   
-----------------------------------------------------------------------------------------------------
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
   0  1.00e+00  2.39e+05  0.00e+00  2.39e+05    2.66e+03      0     SD      0       0                






    




/usr/local/lib/python2.7/dist-packages/SimPEG/Utils/codeutils.py:84: ComplexWarning: Casting complex values to real discards the imaginary part
  values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))






    




------------------------------------------------------------------
0 :    ft     = 2.3966e+05 <= alp*descent     = 2.3947e+05
1 : maxIterLS =      10    <= iterLS          =     10
0 : probSize  =    20000   <= bindingSet      =      0
------------------------- End Linesearch -------------------------
The linesearch got broken. Boo.






    




/usr/local/lib/python2.7/dist-packages/SimPEG/Utils/codeutils.py:109: ComplexWarning: Casting complex values to real discards the imaginary part
  print pad + stopper['str'] % (l<=r,l,r)

In [7]:

    

# import scipy.optimize


    

In [8]:

    

# misfit = lambda x: invProb.evalFunction(x, True, False)
# res = scipy.optimize.minimize(misfit, cInitial.ravel(), jac=True, method='Newton-CG')
# jac = lambda x: -dmisfit.evalDeriv(x)
# res = scipy.optimize.minimize(dmisfit.eval, cInitial.ravel(), jac=jac, method='Newton-CG')


    

In [9]:

    

# import matplotlib.pyplot as plt
# %pylab inline

# plotOpts = {
#     'vmin': min(cTrue.min(), res.x.min()),
#     'vmax': max(cTrue.max(), res.x.max()),
# }

# fig = plt.figure()

# fig.add_subplot(1,2,1)
# plt.imshow(cTrue, **plotOpts)
# plt.colorbar()

# fig.add_subplot(1,2,2)
# plt.imshow(res.x.reshape((nz,nx)), **plotOpts)
# plt.colorbar()


    

In [ ]: