tIPLogRatioFunctional.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tIPLogRatioFunctional.cc (3535B)
---
1 // Copyright (C) 2013, 2014, 2015, 2016, 2017 David Maxwell
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include "IPLogRatioFunctional.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/pism_utilities.hh"
22
23 namespace pism {
24 namespace inverse {
25
26 //! Determine the normalization constant for the functional.
27 /*! Sets the normalization constant \f$c_N\f$ so that
28 \f[
29 J(x)=1
30 \f]
31 if \f$|x| = \mathtt{scale}|u_{\rm obs}| \f$ everywhere.
32 */
33 void IPLogRatioFunctional::normalize(double scale) {
34
35 double value = 0;
36
37 double w = 1.0;
38
39 IceModelVec::AccessList list(m_u_observed);
40
41 if (m_weights) {
42 list.add(*m_weights);
43 }
44
45 for (Points p(*m_grid); p; p.next()) {
46 const int i = p.i(), j = p.j();
47
48 if (m_weights) {
49 w = (*m_weights)(i, j);
50 }
51
52 Vector2 &u_obs_ij = m_u_observed(i, j);
53 double obsMagSq = u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v + m_eps*m_eps;
54
55 double modelMagSq = scale*scale*(u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v) + m_eps*m_eps;
56
57 double v = log(modelMagSq/obsMagSq);
58 value += w*v*v;
59 }
60
61 m_normalization = GlobalSum(m_grid->com, value);
62 }
63
64 void IPLogRatioFunctional::valueAt(IceModelVec2V &x, double *OUTPUT) {
65
66 // The value of the objective
67 double value = 0;
68
69 double w = 1.;
70
71 IceModelVec::AccessList list{&x, &m_u_observed};
72
73 if (m_weights) {
74 list.add(*m_weights);
75 }
76
77 for (Points p(*m_grid); p; p.next()) {
78 const int i = p.i(), j = p.j();
79
80 if (m_weights) {
81 w = (*m_weights)(i, j);
82 }
83 Vector2 &x_ij = x(i, j);
84 Vector2 &u_obs_ij = m_u_observed(i, j);
85 Vector2 u_model_ij = x_ij+u_obs_ij;
86 double obsMagSq = u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v + m_eps*m_eps;
87
88 double modelMagSq = (u_model_ij.u*u_model_ij.u + u_model_ij.v*u_model_ij.v)+m_eps*m_eps;
89 double v = log(modelMagSq/obsMagSq);
90 value += w*v*v;
91 }
92
93 value /= m_normalization;
94
95 GlobalSum(m_grid->com, &value, OUTPUT, 1);
96 }
97
98 void IPLogRatioFunctional::gradientAt(IceModelVec2V &x, IceModelVec2V &gradient) {
99 gradient.set(0);
100
101 double w = 1.;
102
103 IceModelVec::AccessList list{&x, &gradient, &m_u_observed};
104
105 if (m_weights) {
106 list.add(*m_weights);
107 }
108
109 for (Points p(*m_grid); p; p.next()) {
110 const int i = p.i(), j = p.j();
111
112 if (m_weights) {
113 w = (*m_weights)(i, j);
114 }
115 Vector2 &x_ij = x(i, j);
116 Vector2 &u_obs_ij = m_u_observed(i, j);
117 Vector2 u_model_ij = x_ij+u_obs_ij;
118
119 double obsMagSq = u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v + m_eps*m_eps;
120 double modelMagSq = (u_model_ij.u*u_model_ij.u + u_model_ij.v*u_model_ij.v)+m_eps*m_eps;
121 double v = log(modelMagSq/obsMagSq);
122 double dJdw = 2*w*v/modelMagSq;
123
124 gradient(i, j).u = dJdw*2*u_model_ij.u/m_normalization;
125 gradient(i, j).v = dJdw*2*u_model_ij.v/m_normalization;
126 }
127 }
128
129 } // end of namespace inverse
130 } // end of namespace pism