tIPLogRelativeFunctional.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
---
tIPLogRelativeFunctional.cc (3203B)
---
1 // Copyright (C) 2012, 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 "IPLogRelativeFunctional.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}\f$ everywhere.
32 */
33 void IPLogRelativeFunctional::normalize(double scale) {
34
35 // The local value of the weights
36 double value = 0;
37
38 double scale_sq = scale*scale;
39
40 double w = 1.;
41
42 IceModelVec::AccessList list(m_u_observed);
43
44 if (m_weights) {
45 list.add(*m_weights);
46 }
47
48 for (Points p(*m_grid); p; p.next()) {
49 const int i = p.i(), j = p.j();
50
51 Vector2 &u_obs_ij = m_u_observed(i, j);
52 if (m_weights) {
53 w = (*m_weights)(i, j);
54 }
55 double obsMagSq = (u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v) + m_eps*m_eps;
56 value += log(1 + w*scale_sq/obsMagSq);
57 }
58
59 m_normalization = GlobalSum(m_grid->com, value);
60 }
61
62 void IPLogRelativeFunctional::valueAt(IceModelVec2V &x, double *OUTPUT) {
63
64 // The value of the objective
65 double value = 0;
66
67 double w = 1;
68
69 IceModelVec::AccessList list{&x, &m_u_observed};
70 if (m_weights) {
71 list.add(*m_weights);
72 }
73
74 for (Points p(*m_grid); p; p.next()) {
75 const int i = p.i(), j = p.j();
76
77 Vector2 &x_ij = x(i, j);
78 Vector2 &u_obs_ij = m_u_observed(i, j);
79 if (m_weights) {
80 w = (*m_weights)(i, j);
81 }
82 double obsMagSq = (u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v) + m_eps*m_eps;
83 value += log(1 + w*(x_ij.u*x_ij.u + x_ij.v*x_ij.v)/obsMagSq);
84 }
85
86 value /= m_normalization;
87
88 GlobalSum(m_grid->com, &value, OUTPUT, 1);
89 }
90
91 void IPLogRelativeFunctional::gradientAt(IceModelVec2V &x, IceModelVec2V &gradient) {
92 gradient.set(0);
93
94 double w = 1;
95
96 IceModelVec::AccessList list{&x, &gradient, &m_u_observed};
97 if (m_weights) {
98 list.add(*m_weights);
99 }
100
101 for (Points p(*m_grid); p; p.next()) {
102 const int i = p.i(), j = p.j();
103
104 Vector2 &x_ij = x(i, j);
105 Vector2 &u_obs_ij = m_u_observed(i, j);
106 if (m_weights) {
107 w = (*m_weights)(i, j);
108 }
109 double obsMagSq = u_obs_ij.u*u_obs_ij.u + u_obs_ij.v*u_obs_ij.v + m_eps*m_eps;
110 double dJdxsq = w/(obsMagSq + w*(x_ij.u*x_ij.u + x_ij.v*x_ij.v));
111
112 gradient(i, j).u = dJdxsq*2*x_ij.u/m_normalization;
113 gradient(i, j).v = dJdxsq*2*x_ij.v/m_normalization;
114 }
115 }
116
117 } // end of namespace inverse
118 } // end of namespace pism