seis_tomo_joint_with_volume.py

"""
Straight Ray with Volume Data Misfit Term
=========================================

Based on the SEG abstract Heagy, Cockett and Oldenburg, 2014.

Heagy, L. J., Cockett, A. R., & Oldenburg, D. W. (2014, August 5).
Parametrized Inversion Framework for Proppant Volume in a Hydraulically
Fractured Reservoir. SEG Technical Program Expanded Abstracts 2014.
Society of Exploration Geophysicists. doi:10.1190/segam2014-1639.1

This example is a simple joint inversion that consists of a

    - data misfit for the tomography problem
    - data misfit for the volume of the inclusions
      (uses the effective medium theory mapping)
    - model regularization

"""

import numpy as np
import scipy.sparse as sp
import properties
import matplotlib.pyplot as plt

from SimPEG.SEIS import StraightRay
from SimPEG import (
    Mesh, Maps, Utils, Regularization, Optimization,
    InvProblem, Inversion, DataMisfit, Directives, ObjectiveFunction
)


class Volume(ObjectiveFunction.BaseObjectiveFunction):

    """
    A regularization on the volume integral of the model

    .. math::

        \phi_v = \frac{1}{2}|| \int_V m dV - \text{knownVolume} ||^2
    """

    knownVolume = properties.Float("known volume", default=0., min=0.)

    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super(Volume, self).__init__(**kwargs)

    def __call__(self, m):
        return 0.5*(self.estVol(m) - self.knownVolume)**2

    def estVol(self, m):
        return np.inner(self.mesh.vol, m)

    def deriv(self, m):
        # return (self.mesh.vol * np.inner(self.mesh.vol, m))
        return self.mesh.vol * (self.knownVolume - np.inner(self.mesh.vol, m))

    def deriv2(self, m, v=None):
        if v is not None:
            return Utils.mkvc(self.mesh.vol * np.inner(self.mesh.vol, v))
        else:
            # TODO: this is inefficent. It is a fully dense matrix
            return sp.csc_matrix(np.outer(self.mesh.vol, self.mesh.vol))


def run(plotIt=True):

    nC = 40
    de = 1.
    h = np.ones(nC)*de/nC
    M = Mesh.TensorMesh([h, h])

    y = np.linspace(M.vectorCCy[0], M.vectorCCx[-1], int(np.floor(nC/4)))
    rlocs = np.c_[0*y+M.vectorCCx[-1], y]
    rx = StraightRay.Rx(rlocs, None)

    srcList = [
        StraightRay.Src(loc=np.r_[M.vectorCCx[0], yi], rxList=[rx]) for yi in y
    ]

    # phi model
    phi0 = 0
    phi1 = 0.65
    phitrue = Utils.ModelBuilder.defineBlock(
        M.gridCC, [0.4, 0.6], [0.6, 0.4], [phi1, phi0]
    )

    knownVolume = np.sum(phitrue*M.vol)
    print('True Volume: {}'.format(knownVolume))

    # Set up true conductivity model and plot the model transform
    sigma0 = np.exp(1)
    sigma1 = 1e4

    if plotIt:
        fig, ax = plt.subplots(1, 1)
        sigmaMapTest = Maps.SelfConsistentEffectiveMedium(
            nP=1000, sigma0=sigma0, sigma1=sigma1, rel_tol=1e-1, maxIter=150
        )
        testphis = np.linspace(0., 1., 1000)

        sigetest = sigmaMapTest * testphis
        ax.semilogy(testphis, sigetest)
        ax.set_title('Model Transform')
        ax.set_xlabel('$\\varphi$')
        ax.set_ylabel('$\sigma$')

    sigmaMap = Maps.SelfConsistentEffectiveMedium(
        M, sigma0=sigma0, sigma1=sigma1
    )

    # scale the slowness so it is on a ~linear scale
    slownessMap = Maps.LogMap(M) * sigmaMap

    # set up the true sig model and log model dobs
    sigtrue = sigmaMap * phitrue

    # modt = Model.BaseModel(M);
    slownesstrue = slownessMap * phitrue # true model (m = log(sigma))

    # set up the problem and survey
    survey = StraightRay.Survey(srcList)
    problem = StraightRay.Problem(M, slownessMap=slownessMap)
    problem.pair(survey)

    if plotIt:
        fig, ax = plt.subplots(1, 1)
        cb = plt.colorbar(M.plotImage(phitrue, ax=ax)[0], ax=ax)
        survey.plot(ax=ax)
        cb.set_label('$\\varphi$')

    # get observed data
    dobs = survey.makeSyntheticData(phitrue, std=0.03, force=True)
    dpred = survey.dpred(np.zeros(M.nC))

    # objective function pieces
    reg = Regularization.Tikhonov(M)
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmisVol = Volume(mesh=M, knownVolume=knownVolume)
    beta = 0.25
    maxIter = 15

    # without the volume regularization
    opt = Optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
    opt.remember('xc')
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=beta)
    inv = Inversion.BaseInversion(invProb)

    mopt1 = inv.run(np.zeros(M.nC)+1e-16)
    print(
        '\nTotal recovered volume (no vol misfit term in inversion): '
        '{}'.format(
            dmisVol(mopt1)
        )
    )

    # with the volume regularization
    vol_multiplier = 9e4
    reg2 = reg
    dmis2 = dmis + vol_multiplier * dmisVol
    opt2 = Optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
    opt2.remember('xc')
    invProb2 = InvProblem.BaseInvProblem(dmis2, reg2, opt2, beta=beta)
    inv2 = Inversion.BaseInversion(invProb2)

    mopt2 = inv2.run(np.zeros(M.nC)+1e-16)
    print(
        '\nTotal volume (vol misfit term in inversion): {}'.format(
            dmisVol(mopt2)
        )
    )

    # plot results

    if plotIt:

        fig, ax = plt.subplots(1, 1)
        ax.plot(dobs)
        ax.plot(dpred)
        ax.plot(survey.dpred(mopt1), 'o')
        ax.plot(survey.dpred(mopt2), 's')
        ax.legend(['dobs', 'dpred0', 'dpred w/o Vol', 'dpred with Vol'])

        fig, ax = plt.subplots(1, 3, figsize=(16, 4))
        cb0 = plt.colorbar(M.plotImage(phitrue, ax=ax[0])[0], ax=ax[0])
        cb1 = plt.colorbar(M.plotImage(mopt1, ax=ax[1])[0], ax=ax[1])
        cb2 = plt.colorbar(M.plotImage(mopt2, ax=ax[2])[0], ax=ax[2])

        for cb in [cb0, cb1, cb2]:
            cb.set_clim([0., phi1])

        ax[0].set_title('true, vol: {:1.3e}'.format(knownVolume))
        ax[1].set_title(
            'recovered(no Volume term), vol: {:1.3e} '.format(dmisVol(mopt1))
        )
        ax[2].set_title(
            'recovered(with Volume term), vol: {:1.3e} '.format(dmisVol(mopt2))
        )

        plt.tight_layout()

if __name__ == '__main__':
    run()
    plt.show()