Python SimPEG.DataMisfit类代码示例

    def setUp(self):

        cs = 25.
        hx = [(cs, 0, -1.3), (cs, 21), (cs, 0, 1.3)]
        hy = [(cs, 0, -1.3), (cs, 21), (cs, 0, 1.3)]
        hz = [(cs, 0, -1.3), (cs, 20)]
        mesh = Mesh.TensorMesh([hx, hy, hz], x0="CCN")
        blkind0 = Utils.ModelBuilder.getIndicesSphere(
            np.r_[-100., -100., -200.], 75., mesh.gridCC
        )
        blkind1 = Utils.ModelBuilder.getIndicesSphere(
            np.r_[100., 100., -200.], 75., mesh.gridCC
        )
        sigma = np.ones(mesh.nC)*1e-2
        eta = np.zeros(mesh.nC)
        tau = np.ones_like(sigma)*1.
        eta[blkind0] = 0.1
        eta[blkind1] = 0.1
        tau[blkind0] = 0.1
        tau[blkind1] = 0.01

        x = mesh.vectorCCx[(mesh.vectorCCx > -155.) & (mesh.vectorCCx < 155.)]
        y = mesh.vectorCCx[(mesh.vectorCCy > -155.) & (mesh.vectorCCy < 155.)]
        Aloc = np.r_[-200., 0., 0.]
        Bloc = np.r_[200., 0., 0.]
        M = Utils.ndgrid(x-25., y, np.r_[0.])
        N = Utils.ndgrid(x+25., y, np.r_[0.])

        times = np.arange(10)*1e-3 + 1e-3
        rx = SIP.Rx.Dipole(M, N, times)
        src = SIP.Src.Dipole([rx], Aloc, Bloc)
        survey = SIP.Survey([src])
        wires = Maps.Wires(('eta', mesh.nC), ('taui', mesh.nC))
        problem = SIP.Problem3D_N(
            mesh,
            sigma=sigma,
            etaMap=wires.eta,
            tauiMap=wires.taui
        )
        problem.Solver = Solver
        problem.pair(survey)
        mSynth = np.r_[eta, 1./tau]
        survey.makeSyntheticData(mSynth)
        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

    def test_basic_inversion(self):

        """
        Test to see if inversion recovers model
        """

        h = [(2, 30)]
        meshObj = Mesh.TensorMesh((h, h, [(2, 10)]), x0='CCN')

        mod = 0.00025*np.ones(meshObj.nC)
        mod[(meshObj.gridCC[:, 0] > -4.) &
            (meshObj.gridCC[:, 1] > -4.) &
            (meshObj.gridCC[:, 0] < 4.) &
            (meshObj.gridCC[:, 1] < 4.)] = 0.001

        times = np.logspace(-4, -2, 5)
        waveObj = VRM.WaveformVRM.SquarePulse(delt=0.02)

        x, y = np.meshgrid(np.linspace(-17, 17, 16), np.linspace(-17, 17, 16))
        x, y, z = mkvc(x), mkvc(y), 0.5*np.ones(np.size(x))
        rxList = [VRM.Rx.Point(np.c_[x, y, z], times=times, fieldType='dbdt', fieldComp='z')]

        txNodes = np.array([[-20, -20, 0.001],
                            [20, -20, 0.001],
                            [20, 20, 0.001],
                            [-20, 20, 0.01],
                            [-20, -20, 0.001]])
        txList = [VRM.Src.LineCurrent(rxList, txNodes, 1., waveObj)]

        Survey = VRM.Survey(txList)
        Survey.t_active = np.zeros(Survey.nD, dtype=bool)
        Survey.set_active_interval(-1e6, 1e6)
        Problem = VRM.Problem_Linear(meshObj, ref_factor=2)
        Problem.pair(Survey)
        Survey.makeSyntheticData(mod)
        Survey.eps = 1e-11

        dmis = DataMisfit.l2_DataMisfit(Survey)
        W = mkvc((np.sum(np.array(Problem.A)**2, axis=0)))**0.25
        reg = Regularization.Simple(
            meshObj, alpha_s=0.01, alpha_x=1., alpha_y=1., alpha_z=1., cell_weights=W
            )
        opt = Optimization.ProjectedGNCG(
            maxIter=20, lower=0., upper=1e-2, maxIterLS=20, tolCG=1e-4
            )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
        directives = [
            Directives.BetaSchedule(coolingFactor=2, coolingRate=1),
            Directives.TargetMisfit()
        ]
        inv = Inversion.BaseInversion(invProb, directiveList=directives)

        m0 = 1e-6*np.ones(len(mod))
        mrec = inv.run(m0)

        dmis_final = np.sum((dmis.W.diagonal()*(Survey.dobs - Problem.fields(mrec)))**2)
        mod_err_2 = np.sqrt(np.sum((mrec-mod)**2))/np.size(mod)
        mod_err_inf = np.max(np.abs(mrec-mod))

        self.assertTrue(dmis_final < Survey.nD and mod_err_2 < 5e-6 and mod_err_inf < np.max(mod))

    def setUp(self):
        mesh = Mesh.TensorMesh([30, 30], x0=[-0.5, -1.])
        sigma = np.ones(mesh.nC)
        model = np.log(sigma)

        prob = DC.Problem3D_CC(mesh, rhoMap=Maps.ExpMap(mesh))

        rx = DC.Rx.Pole(
            Utils.ndgrid([mesh.vectorCCx, np.r_[mesh.vectorCCy.max()]])
        )
        src = DC.Src.Dipole(
            [rx], np.r_[-0.25, mesh.vectorCCy.max()],
            np.r_[0.25, mesh.vectorCCy.max()]
        )
        survey = DC.Survey([src])

        prob.pair(survey)

        dobs = survey.makeSyntheticData(model)

        dmis = DataMisfit.l2_DataMisfit(survey)

        self.model = model
        self.mesh = mesh
        self.survey = survey
        self.prob = prob
        self.dobs = dobs
        self.dmis = dmis

def run(N=100, plotIt=True):

    np.random.seed(1)

    mesh = Mesh.TensorMesh([N])

    nk = 20
    jk = np.linspace(1., 60., nk)
    p = -0.25
    q = 0.25

    def g(k):
        return (
            np.exp(p*jk[k]*mesh.vectorCCx) *
            np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
        )

    G = np.empty((nk, mesh.nC))

    for i in range(nk):
        G[i, :] = g(i)

    mtrue = np.zeros(mesh.nC)
    mtrue[mesh.vectorCCx > 0.3] = 1.
    mtrue[mesh.vectorCCx > 0.45] = -0.5
    mtrue[mesh.vectorCCx > 0.6] = 0

    prob = Problem.LinearProblem(mesh, G=G)
    survey = Survey.LinearSurvey()
    survey.pair(prob)
    survey.makeSyntheticData(mtrue, std=0.01)

    M = prob.mesh

    reg = Regularization.Tikhonov(mesh, alpha_s=1., alpha_x=1.)
    dmis = DataMisfit.l2_DataMisfit(survey)
    opt = Optimization.InexactGaussNewton(maxIter=60)
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
    directives = [
        Directives.BetaEstimate_ByEig(beta0_ratio=1e-2),
        Directives.TargetMisfit()
    ]
    inv = Inversion.BaseInversion(invProb, directiveList=directives)
    m0 = np.zeros_like(survey.mtrue)

    mrec = inv.run(m0)

    if plotIt:
        fig, axes = plt.subplots(1, 2, figsize=(12*1.2, 4*1.2))
        for i in range(prob.G.shape[0]):
            axes[0].plot(prob.G[i, :])
        axes[0].set_title('Columns of matrix G')

        axes[1].plot(M.vectorCCx, survey.mtrue, 'b-')
        axes[1].plot(M.vectorCCx, mrec, 'r-')
        axes[1].legend(('True Model', 'Recovered Model'))
        axes[1].set_ylim([-2, 2])

    return prob, survey, mesh, mrec

    def setUp(self, parallel=True):
        frequency = np.array([900, 7200, 56000], dtype=float)
        hz = get_vertical_discretization_frequency(
            frequency, sigma_background=1./10.
        )
        n_sounding = 10
        dx = 20.
        hx = np.ones(n_sounding) * dx
        mesh = Mesh.TensorMesh([hx, hz], x0='00')
        inds = mesh.gridCC[:, 1] < 25
        inds_1 = mesh.gridCC[:, 1] < 50
        sigma = np.ones(mesh.nC) * 1./100.
        sigma[inds_1] = 1./10.
        sigma[inds] = 1./50.
        sigma_em1d = sigma.reshape(mesh.vnC, order='F').flatten()
        mSynth = np.log(sigma_em1d)

        x = mesh.vectorCCx
        y = np.zeros_like(x)
        z = np.ones_like(x) * 30.
        rx_locations = np.c_[x, y, z]
        src_locations = np.c_[x, y, z]
        topo = np.c_[x, y, z-30.].astype(float)
        mapping = Maps.ExpMap(mesh)
        survey = GlobalEM1DSurveyFD(
            rx_locations=rx_locations,
            src_locations=src_locations,
            frequency=frequency,
            offset=np.ones_like(frequency) * 8.,
            src_type="VMD",
            rx_type="Hz",
            field_type='secondary',
            topo=topo
        )

        problem = GlobalEM1DProblemFD(
            [], sigmaMap=mapping, hz=hz,
            parallel=parallel, n_cpu=2
        )
        problem.pair(survey)
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )

        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=0.)
        inv = Inversion.BaseInversion(invProb)
        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

    def setUp(self):
        mesh = Mesh.TensorMesh([30, 30], x0=[-0.5, -1.])
        sigma = np.random.rand(mesh.nC)
        model = np.log(sigma)

        prob = DC.Problem3D_CC(mesh, rhoMap=Maps.ExpMap(mesh))
        prob1 = DC.Problem3D_CC(mesh, rhoMap=Maps.ExpMap(mesh))

        rx = DC.Rx.Pole(
            Utils.ndgrid([mesh.vectorCCx, np.r_[mesh.vectorCCy.max()]])
        )
        rx1 = DC.Rx.Pole(
            Utils.ndgrid([mesh.vectorCCx, np.r_[mesh.vectorCCy.min()]])
        )
        src = DC.Src.Dipole(
            [rx], np.r_[-0.25, mesh.vectorCCy.max()],
            np.r_[0.25, mesh.vectorCCy.max()]
        )
        src1 = DC.Src.Dipole(
            [rx1], np.r_[-0.25, mesh.vectorCCy.max()],
            np.r_[0.25, mesh.vectorCCy.max()]
        )
        survey = DC.Survey([src])
        prob.pair(survey)

        survey1 = DC.Survey([src1])
        prob1.pair(survey1)

        dobs0 = survey.makeSyntheticData(model)
        dobs1 = survey1.makeSyntheticData(model)

        self.mesh = mesh
        self.model = model

        self.survey0 = survey
        self.prob0 = prob

        self.survey1 = survey1
        self.prob1 = prob1

        self.dmis0 = DataMisfit.l2_DataMisfit(self.survey0)
        self.dmis1 = DataMisfit.l2_DataMisfit(self.survey1)

        self.dmiscobmo = self.dmis0 + self.dmis1

    def setUp(self):

        cs = 12.5
        hx = [(cs, 7, -1.3), (cs, 61), (cs, 7, 1.3)]
        hy = [(cs, 7, -1.3), (cs, 20)]
        mesh = Mesh.TensorMesh([hx, hy], x0="CN")

        # x = np.linspace(-200, 200., 20)
        x = np.linspace(-200, 200., 2)
        M = Utils.ndgrid(x-12.5, np.r_[0.])
        N = Utils.ndgrid(x+12.5, np.r_[0.])

        A0loc = np.r_[-150, 0.]
        A1loc = np.r_[-130, 0.]
        B0loc = np.r_[-130, 0.]
        B1loc = np.r_[-110, 0.]

        rx = DC.Rx.Dipole_ky(M, N)
        src0 = DC.Src.Dipole([rx], A0loc, B0loc)
        src1 = DC.Src.Dipole([rx], A1loc, B1loc)
        survey = IP.Survey([src0, src1])

        sigma = np.ones(mesh.nC) * 1.
        problem = IP.Problem2D_CC(
            mesh, sigma=sigma, etaMap=Maps.IdentityMap(mesh),
            verbose=False
        )
        problem.pair(survey)

        mSynth = np.ones(mesh.nC)*0.1
        survey.makeSyntheticData(mSynth)
        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

    def setUp(self):
        mesh = Mesh.TensorMesh([4, 4, 4])

        # Magnetic inducing field parameter (A,I,D)
        B = [50000, 90, 0]

        # Create a MAGsurvey
        rx = PF.BaseMag.RxObs(
            np.vstack([[0.25, 0.25, 0.25], [-0.25, -0.25, 0.25]])
        )
        srcField = PF.BaseMag.SrcField([rx], param=(B[0], B[1], B[2]))
        survey = PF.BaseMag.LinearSurvey(srcField)

        # Create the forward model operator
        prob = PF.Magnetics.MagneticIntegral(
            mesh, chiMap=Maps.IdentityMap(mesh)
        )

        # Pair the survey and problem
        survey.pair(prob)

        # Compute forward model some data
        m = np.random.rand(mesh.nC)
        survey.makeSyntheticData(m)

        reg = Regularization.Sparse(mesh)
        reg.mref = np.zeros(mesh.nC)

        wr = np.sum(prob.G**2., axis=0)**0.5
        reg.cell_weights = wr
        reg.norms = [0, 1, 1, 1]
        reg.eps_p, reg.eps_q = 1e-3, 1e-3

        # Data misfit function
        dmis = DataMisfit.l2_DataMisfit(survey)
        dmis.W = 1./survey.std

        # Add directives to the inversion
        opt = Optimization.ProjectedGNCG(
            maxIter=2, lower=-10., upper=10.,
            maxIterCG=2
        )

        invProb = InvProblem.BaseInvProblem(dmis, reg, opt)

        self.mesh = mesh
        self.invProb = invProb

    def setUp(self):

        mesh = Mesh.TensorMesh([20, 20, 20], "CCN")
        sigma = np.ones(mesh.nC)*1./100.
        actind = mesh.gridCC[:, 2] < -0.2
        # actMap = Maps.InjectActiveCells(mesh, actind, 0.)

        xyzM = Utils.ndgrid(np.ones_like(mesh.vectorCCx[:-1])*-0.4, np.ones_like(mesh.vectorCCy)*-0.4, np.r_[-0.3])
        xyzN = Utils.ndgrid(mesh.vectorCCx[1:], mesh.vectorCCy, np.r_[-0.3])

        problem = SP.Problem_CC(mesh, sigma=sigma, qMap=Maps.IdentityMap(mesh), Solver=PardisoSolver)
        rx = SP.Rx.Dipole(xyzN, xyzM)
        src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC), mesh=mesh,
                                       modelType="CurrentSource")
        survey = SP.Survey([src])
        survey.pair(problem)

        q = np.zeros(mesh.nC)
        inda = Utils.closestPoints(mesh, np.r_[-0.5, 0., -0.8])
        indb = Utils.closestPoints(mesh, np.r_[0.5, 0., -0.8])
        q[inda] = 1.
        q[indb] = -1.

        mSynth = q.copy()
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Simple(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e-2)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

    def setUp(self):

        cs = 12.5
        hx = [(cs, 2, -1.3), (cs, 61), (cs, 2, 1.3)]
        hy = [(cs, 2, -1.3), (cs, 20)]
        mesh = Mesh.TensorMesh([hx, hy], x0="CN")
        x = np.linspace(-135, 250., 20)
        M = Utils.ndgrid(x-12.5, np.r_[0.])
        N = Utils.ndgrid(x+12.5, np.r_[0.])
        A0loc = np.r_[-150, 0.]
        A1loc = np.r_[-130, 0.]
        # rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
        rx = DC.Rx.Dipole_ky(M, N)
        src0 = DC.Src.Pole([rx], A0loc)
        src1 = DC.Src.Pole([rx], A1loc)
        survey = DC.Survey_ky([src0, src1])
        problem = DC.Problem2D_N(
            mesh, rhoMap=Maps.IdentityMap(mesh),
            Solver=Solver
        )
        problem.pair(survey)

        mSynth = np.ones(mesh.nC)*1.
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

    def setUp(self):

        aSpacing = 2.5
        nElecs = 10

        surveySize = nElecs*aSpacing - aSpacing
        cs = surveySize / nElecs / 4

        mesh = Mesh.TensorMesh([
            [(cs, 10, -1.3), (cs, surveySize / cs), (cs, 10, 1.3)],
            [(cs, 3, -1.3), (cs, 3, 1.3)],
            # [(cs, 5, -1.3), (cs, 10)]
        ], 'CN')

        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
        survey = DC.Survey(srcList)
        problem = DC.Problem3D_N(
            mesh, rhoMap=Maps.IdentityMap(mesh), storeJ=True
            )
        problem.pair(survey)

        mSynth = np.ones(mesh.nC)
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis

cb = plt.colorbar(dat[0], ax=ax[0])
ax[0].set_title("Vertical section")
cb.set_label("Conductivity (S/m)")
ax[0].set_xlabel('Easting (m)')
ax[0].set_ylabel('Depth (m)')
ax[0].set_xlim(-1000., 1000.)
ax[0].set_ylim(-500., 0.)

###############################################################################
# Step 6
# ------
#
# Run inversion

regmesh = Mesh.TensorMesh([31])
dmis = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Tikhonov(regmesh)
opt = Optimization.InexactGaussNewton(maxIter=7, tolX=1e-15)
opt.remember('xc')
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaEstimate_ByEig(beta0_ratio=1e1)
betaSched = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaSched])

# Choose an initial starting model of the background conductivity
m0 = np.log(np.ones(mapping.nP)*sighalf)
mopt = inv.run(m0)


###############################################################################
# Step 7

def resolve_1Dinversions(
    mesh, dobs, src_height, freqs, m0, mref, mapping,
    std=0.08, floor=1e-14, rxOffset=7.86
):
    """
    Perform a single 1D inversion for a RESOLVE sounding for Horizontal
    Coplanar Coil data (both real and imaginary).

    :param discretize.CylMesh mesh: mesh used for the forward simulation
    :param numpy.array dobs: observed data
    :param float src_height: height of the source above the ground
    :param numpy.array freqs: frequencies
    :param numpy.array m0: starting model
    :param numpy.array mref: reference model
    :param Maps.IdentityMap mapping: mapping that maps the model to electrical conductivity
    :param float std: percent error used to construct the data misfit term
    :param float floor: noise floor used to construct the data misfit term
    :param float rxOffset: offset between source and receiver.
    """

    # ------------------- Forward Simulation ------------------- #
    # set up the receivers
    bzr = EM.FDEM.Rx.Point_bSecondary(
        np.array([[rxOffset, 0., src_height]]),
        orientation='z',
        component='real'
    )

    bzi = EM.FDEM.Rx.Point_b(
        np.array([[rxOffset, 0., src_height]]),
        orientation='z',
        component='imag'
    )

    # source location
    srcLoc = np.array([0., 0., src_height])
    srcList = [
        EM.FDEM.Src.MagDipole([bzr, bzi], freq, srcLoc, orientation='Z')
        for freq in freqs
    ]

    # construct a forward simulation
    survey = EM.FDEM.Survey(srcList)
    prb = EM.FDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=PardisoSolver)
    prb.pair(survey)

    # ------------------- Inversion ------------------- #
    # data misfit term
    survey.dobs = dobs
    dmisfit = DataMisfit.l2_DataMisfit(survey)
    uncert = abs(dobs) * std + floor
    dmisfit.W = 1./uncert

    # regularization
    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
    reg = Regularization.Simple(regMesh)
    reg.mref = mref

    # optimization
    opt = Optimization.InexactGaussNewton(maxIter=10)

    # statement of the inverse problem
    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)

    # Inversion directives and parameters
    target = Directives.TargetMisfit()
    inv = Inversion.BaseInversion(invProb, directiveList=[target])

    invProb.beta = 2.   # Fix beta in the nonlinear iterations
    reg.alpha_s = 1e-3
    reg.alpha_x = 1.
    prb.counter = opt.counter = Utils.Counter()
    opt.LSshorten = 0.5
    opt.remember('xc')

    # run the inversion
    mopt = inv.run(m0)
    return mopt, invProb.dpred, survey.dobs

    def setUp(self):

        ndv = -100
        # Create a self.mesh
        dx = 5.

        hxind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
        hyind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
        hzind = [(dx, 5, -1.3), (dx, 6)]

        self.mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')

        # Get index of the center
        midx = int(self.mesh.nCx/2)
        midy = int(self.mesh.nCy/2)

        # Lets create a simple Gaussian topo and set the active cells
        [xx, yy] = np.meshgrid(self.mesh.vectorNx, self.mesh.vectorNy)
        zz = -np.exp((xx**2 + yy**2) / 75**2) + self.mesh.vectorNz[-1]

        # Go from topo to actv cells
        topo = np.c_[Utils.mkvc(xx), Utils.mkvc(yy), Utils.mkvc(zz)]
        actv = Utils.surface2ind_topo(self.mesh, topo, 'N')
        actv = np.where(actv)[0]

        # Create active map to go from reduce space to full
        self.actvMap = Maps.InjectActiveCells(self.mesh, actv, -100)
        nC = len(actv)

        # Create and array of observation points
        xr = np.linspace(-20., 20., 20)
        yr = np.linspace(-20., 20., 20)
        X, Y = np.meshgrid(xr, yr)

        # Move the observation points 5m above the topo
        Z = -np.exp((X**2 + Y**2) / 75**2) + self.mesh.vectorNz[-1] + 5.

        # Create a MAGsurvey
        locXYZ = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
        rxLoc = PF.BaseGrav.RxObs(locXYZ)
        srcField = PF.BaseGrav.SrcField([rxLoc])
        survey = PF.BaseGrav.LinearSurvey(srcField)

        # We can now create a density model and generate data
        # Here a simple block in half-space
        model = np.zeros((self.mesh.nCx, self.mesh.nCy, self.mesh.nCz))
        model[(midx-2):(midx+2), (midy-2):(midy+2), -6:-2] = 0.5
        model = Utils.mkvc(model)
        self.model = model[actv]

        # Create active map to go from reduce set to full
        actvMap = Maps.InjectActiveCells(self.mesh, actv, ndv)

        # Create reduced identity map
        idenMap = Maps.IdentityMap(nP=nC)

        # Create the forward model operator
        prob = PF.Gravity.GravityIntegral(
            self.mesh,
            rhoMap=idenMap,
            actInd=actv
        )

        # Pair the survey and problem
        survey.pair(prob)

        # Compute linear forward operator and compute some data
        d = prob.fields(self.model)

        # Add noise and uncertainties (1nT)
        data = d + np.random.randn(len(d))*0.001
        wd = np.ones(len(data))*.001

        survey.dobs = data
        survey.std = wd

        # PF.Gravity.plot_obs_2D(survey.srcField.rxList[0].locs, d=data)

        # Create sensitivity weights from our linear forward operator
        wr = PF.Magnetics.get_dist_wgt(self.mesh, locXYZ, actv, 2., 2.)
        wr = wr**2.

        # Create a regularization
        reg = Regularization.Sparse(self.mesh, indActive=actv, mapping=idenMap)
        reg.cell_weights = wr
        reg.norms = np.c_[0, 0, 0, 0]
        reg.gradientType = 'component'
        # reg.eps_p, reg.eps_q = 5e-2, 1e-2

        # Data misfit function
        dmis = DataMisfit.l2_DataMisfit(survey)
        dmis.W = 1/wd

        # Add directives to the inversion
        opt = Optimization.ProjectedGNCG(maxIter=100, lower=-1., upper=1.,
                                         maxIterLS=20, maxIterCG=10,
                                         tolCG=1e-3)
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e+8)

        # Here is where the norms are applied
#.........这里部分代码省略.........

    def setUp(self):

        cs = 25.
        hx = [(cs, 0, -1.3), (cs, 21), (cs, 0, 1.3)]
        hy = [(cs, 0, -1.3), (cs, 21), (cs, 0, 1.3)]
        hz = [(cs, 0, -1.3), (cs, 20), (cs, 0, 1.3)]
        mesh = Mesh.TensorMesh([hx, hy, hz], x0="CCC")
        blkind0 = Utils.ModelBuilder.getIndicesSphere(
            np.r_[-100., -100., -200.], 75., mesh.gridCC
        )
        blkind1 = Utils.ModelBuilder.getIndicesSphere(
            np.r_[100., 100., -200.], 75., mesh.gridCC
        )
        sigma = np.ones(mesh.nC)*1e-2
        airind = mesh.gridCC[:, 2] > 0.
        sigma[airind] = 1e-8
        eta = np.zeros(mesh.nC)
        tau = np.ones_like(sigma) * 1.
        c = np.ones_like(sigma) * 0.5

        eta[blkind0] = 0.1
        eta[blkind1] = 0.1
        tau[blkind0] = 0.1
        tau[blkind1] = 0.01

        actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.)
        actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.)
        actmapc = Maps.InjectActiveCells(mesh, ~airind, 1.)

        x = mesh.vectorCCx[(mesh.vectorCCx > -155.) & (mesh.vectorCCx < 155.)]
        y = mesh.vectorCCy[(mesh.vectorCCy > -155.) & (mesh.vectorCCy < 155.)]
        Aloc = np.r_[-200., 0., 0.]
        Bloc = np.r_[200., 0., 0.]
        M = Utils.ndgrid(x-25., y, np.r_[0.])
        N = Utils.ndgrid(x+25., y, np.r_[0.])

        times = np.arange(10)*1e-3 + 1e-3
        rx = SIP.Rx.Dipole(M, N, times)
        src = SIP.Src.Dipole([rx], Aloc, Bloc)
        survey = SIP.Survey([src])

        wires = Maps.Wires(('eta', actmapeta.nP), ('taui', actmaptau.nP), ('c', actmapc.nP))
        problem = SIP.Problem3D_N(
            mesh,
            sigma=sigma,
            etaMap=actmapeta*wires.eta,
            tauiMap=actmaptau*wires.taui,
            cMap=actmapc*wires.c,
            actinds=~airind,
            storeJ = True,
            verbose=False
        )

        problem.Solver = Solver
        problem.pair(survey)
        mSynth = np.r_[eta[~airind], 1./tau[~airind], c[~airind]]
        survey.makeSyntheticData(mSynth)
        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg_eta = Regularization.Simple(mesh, mapping=wires.eta, indActive=~airind)
        reg_taui = Regularization.Simple(mesh, mapping=wires.taui, indActive=~airind)
        reg_c = Regularization.Simple(mesh, mapping=wires.c, indActive=~airind)
        reg = reg_eta + reg_taui + reg_c
        opt = Optimization.InexactGaussNewton(
            maxIterLS=20, maxIter=10, tolF=1e-6,
            tolX=1e-6, tolG=1e-6, maxIterCG=6
        )
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis