forked from simpeg/simpeg
/
DataMisfit.py
169 lines (127 loc) · 4.47 KB
/
DataMisfit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
from __future__ import print_function
import numpy as np
import properties
from . import Utils
from . import Survey
from . import ObjectiveFunction
class BaseDataMisfit(ObjectiveFunction.L2ObjectiveFunction):
"""
BaseDataMisfit
.. note::
You should inherit from this class to create your own data misfit
term.
"""
debug = False #: Print debugging information
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
_hasFields = True #: Data Misfits take fields, handy to store them
def __init__(self, survey, **kwargs):
assert survey.ispaired, 'The survey must be paired to a problem.'
if isinstance(survey, Survey.BaseSurvey):
self.survey = survey
self.prob = survey.prob
super(BaseDataMisfit, self).__init__(**kwargs)
@property
def nP(self):
if self._mapping is not None:
return self.mapping.nP
elif self.prob.model is not None:
return len(self.prob.model)
else:
return '*'
@property
def Wd(self):
raise AttributeError(
'The `Wd` property been depreciated, please use: `W` instead'
)
class l2_DataMisfit(BaseDataMisfit):
"""
The data misfit with an l_2 norm:
.. math::
\mu_\\text{data} = {1\over 2}\left|
\mathbf{W}_d (\mathbf{d}_\\text{pred} -
\mathbf{d}_\\text{obs}) \\right|_2^2
"""
std = 0.05 #: default standard deviation if not provided by survey
eps = None #: default floor
eps_factor = 1e-5 #: factor to multiply by the norm of the data to create floor
def __init__(self, survey, **kwargs):
BaseDataMisfit.__init__(self, survey, **kwargs)
if self.std is None:
if getattr(self.survey, 'std', None) is not None:
print(
'SimPEG.DataMisfit.l2_DataMisfit assigning default std '
'of 5%'
)
else:
self.std = self.survey.std
if self.eps is None:
if getattr(self.survey, 'eps', None) is None:
print(
'SimPEG.DataMisfit.l2_DataMisfit assigning default eps '
'of 1e-5 * ||dobs||'
)
self.eps = (
np.linalg.norm(Utils.mkvc(survey.dobs), 2)*self.eps_factor
) # default
else:
self.eps = self.survey.eps
@property
def W(self):
"""W
The data weighting matrix.
The default is based on the norm of the data plus a noise floor.
:rtype: scipy.sparse.csr_matrix
:return: W
"""
if getattr(self, '_W', None) is None:
survey = self.survey
self._W = Utils.sdiag(1/(abs(survey.dobs)*self.std+self.eps))
return self._W
@W.setter
def W(self, value):
if len(value.shape) < 2:
value = Utils.sdiag(value)
assert value.shape == (self.survey.nD, self.survey.nD), (
'W must have shape ({nD},{nD}), not ({val0}, val{1})'.format(
nD=self.survey.nD, val0=value.shape[0], val1=value.shape[1]
)
)
self._W = value
@Utils.timeIt
def __call__(self, m, f=None):
"__call__(m, f=None)"
if f is None:
f = self.prob.fields(m)
R = self.W * self.survey.residual(m, f)
return 0.5*np.vdot(R, R)
@Utils.timeIt
def deriv(self, m, f=None):
"""
deriv(m, f=None)
Derivative of the data misfit
.. math::
\mathbf{J}^{\top} \mathbf{W}^{\top} \mathbf{W}
(\mathbf{d} - \mathbf{d}^{obs})
:param numpy.ndarray m: model
:param SimPEG.Fields.Fields f: fields object
"""
if f is None:
f = self.prob.fields(m)
return self.prob.Jtvec(
m, self.W.T * (self.W * self.survey.residual(m, f=f)), f=f
)
@Utils.timeIt
def deriv2(self, m, v, f=None):
"""
deriv2(m, v, f=None)
.. math::
\mathbf{J}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{J}
:param numpy.ndarray m: model
:param numpy.ndarray v: vector
:param SimPEG.Fields.Fields f: fields object
"""
if f is None:
f = self.prob.fields(m)
return self.prob.Jtvec_approx(
m, self.W * (self.W * self.prob.Jvec_approx(m, v, f=f)), f=f
)