tIPMeanSquareFunctional.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tIPMeanSquareFunctional.cc (6258B)
---
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 "IPMeanSquareFunctional.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/pism_utilities.hh"
22
23 namespace pism {
24 namespace inverse {
25
26 //! Implicitly set the normalization constant for the functional.
27 /*! The normalization constant is selected so that if an input
28 IceModelVec2V has component vectors all of length \a scale, then the funtional value will be 1. I.e.
29 \f[
30 c_N^{-1} = \sum_{i} w_i {\tt scale}^2.
31 \f]*/
32 void IPMeanSquareFunctional2V::normalize(double scale) {
33
34 // The local value of the weights
35 double value = 0;
36
37 if (m_weights) {
38 IceModelVec::AccessList list{m_weights};
39
40 for (Points p(*m_grid); p; p.next()) {
41 const int i = p.i(), j = p.j();
42
43 value += (*m_weights)(i, j);
44 }
45 } else {
46 for (Points p(*m_grid); p; p.next()) {
47 value += 1;
48 }
49 }
50
51 m_normalization = GlobalSum(m_grid->com, value);
52 m_normalization *= (scale*scale);
53 }
54
55 void IPMeanSquareFunctional2V::valueAt(IceModelVec2V &x, double *OUTPUT) {
56
57 // The value of the objective
58 double value = 0;
59
60 IceModelVec::AccessList list{&x};
61
62 if (m_weights) {
63 list.add(*m_weights);
64 for (Points p(*m_grid); p; p.next()) {
65 const int i = p.i(), j = p.j();
66
67 Vector2 &x_ij = x(i, j);
68 value += (x_ij.u*x_ij.u + x_ij.v*x_ij.v)*(*m_weights)(i, j);
69 }
70 } else {
71 for (Points p(*m_grid); p; p.next()) {
72 const int i = p.i(), j = p.j();
73
74 Vector2 &x_ij = x(i, j);
75 value += (x_ij.u*x_ij.u + x_ij.v*x_ij.v);
76 }
77 }
78 value /= m_normalization;
79
80 GlobalSum( m_grid->com, &value, OUTPUT, 1);
81 }
82
83 void IPMeanSquareFunctional2V::dot(IceModelVec2V &a, IceModelVec2V &b, double *OUTPUT) {
84
85 // The value of the objective
86 double value = 0;
87
88 IceModelVec::AccessList list{&a, &b};
89
90 if (m_weights) {
91 list.add(*m_weights);
92 for (Points p(*m_grid); p; p.next()) {
93 const int i = p.i(), j = p.j();
94
95 Vector2 &a_ij = a(i, j);
96 Vector2 &b_ij = b(i, j);
97 value += (a_ij.u*b_ij.u + a_ij.v*b_ij.v)*(*m_weights)(i, j);
98 }
99 } else {
100 for (Points p(*m_grid); p; p.next()) {
101 const int i = p.i(), j = p.j();
102
103 Vector2 &a_ij = a(i, j);
104 Vector2 &b_ij = b(i, j);
105 value += (a_ij.u*b_ij.u + a_ij.v*b_ij.v);
106 }
107 }
108 value /= m_normalization;
109
110 GlobalSum( m_grid->com, &value, OUTPUT, 1);
111 }
112
113 void IPMeanSquareFunctional2V::gradientAt(IceModelVec2V &x, IceModelVec2V &gradient) {
114 gradient.set(0);
115
116 IceModelVec::AccessList list{&x, &gradient};
117
118 if (m_weights) {
119 list.add(*m_weights);
120 for (Points p(*m_grid); p; p.next()) {
121 const int i = p.i(), j = p.j();
122
123 gradient(i, j).u = 2*x(i, j).u*(*m_weights)(i, j) / m_normalization;
124 gradient(i, j).v = 2*x(i, j).v*(*m_weights)(i, j) / m_normalization;
125 }
126 } else {
127 for (Points p(*m_grid); p; p.next()) {
128 const int i = p.i(), j = p.j();
129
130 gradient(i, j).u = 2*x(i, j).u / m_normalization;
131 gradient(i, j).v = 2*x(i, j).v / m_normalization;
132 }
133 }
134 }
135
136 //! Implicitly set the normalization constant for the functional.
137 /*! The normalization constant is selected so that if an input
138 IceModelVec2S has entries all equal to \a scale, then the funtional value will be 1. I.e.
139 \f[
140 c_N^{-1} = \sum_{i} w_i {\tt scale}^2.
141 \f]*/
142 void IPMeanSquareFunctional2S::normalize(double scale) {
143
144 // The local value of the weights
145 double value = 0;
146
147 if (m_weights) {
148 IceModelVec::AccessList list(*m_weights);
149 for (Points p(*m_grid); p; p.next()) {
150 const int i = p.i(), j = p.j();
151
152 value += (*m_weights)(i, j);
153 }
154 } else {
155 for (Points p(*m_grid); p; p.next()) {
156 value += 1;
157 }
158 }
159
160 m_normalization = GlobalSum(m_grid->com, value);
161 m_normalization *= (scale*scale);
162 }
163
164 void IPMeanSquareFunctional2S::valueAt(IceModelVec2S &x, double *OUTPUT) {
165
166 // The value of the objective
167 double value = 0;
168
169 IceModelVec::AccessList list(x);
170
171 if (m_weights) {
172 list.add(*m_weights);
173 for (Points p(*m_grid); p; p.next()) {
174 const int i = p.i(), j = p.j();
175
176 double &x_ij = x(i, j);
177 value += x_ij*x_ij*(*m_weights)(i, j);
178 }
179 } else {
180 for (Points p(*m_grid); p; p.next()) {
181 const int i = p.i(), j = p.j();
182
183 double &x_ij = x(i, j);
184 value += x_ij*x_ij;
185 }
186 }
187 value /= m_normalization;
188
189 GlobalSum(m_grid->com, &value, OUTPUT, 1);
190 }
191
192 void IPMeanSquareFunctional2S::dot(IceModelVec2S &a, IceModelVec2S &b, double *OUTPUT) {
193
194 // The value of the objective
195 double value = 0;
196
197 IceModelVec::AccessList list{&a, &b};
198
199 if (m_weights) {
200 list.add(*m_weights);
201 for (Points p(*m_grid); p; p.next()) {
202 const int i = p.i(), j = p.j();
203
204 value += (a(i, j)*b(i, j))*(*m_weights)(i, j);
205 }
206 } else {
207 for (Points p(*m_grid); p; p.next()) {
208 const int i = p.i(), j = p.j();
209
210 value += (a(i, j)*b(i, j));
211 }
212 }
213 value /= m_normalization;
214
215 GlobalSum(m_grid->com, &value, OUTPUT, 1);
216 }
217
218
219 void IPMeanSquareFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &gradient) {
220 gradient.set(0);
221
222 IceModelVec::AccessList list{&x, &gradient};
223
224 if (m_weights) {
225 list.add(*m_weights);
226 for (Points p(*m_grid); p; p.next()) {
227 const int i = p.i(), j = p.j();
228
229 gradient(i, j) = 2*x(i, j)*(*m_weights)(i, j) / m_normalization;
230 }
231 } else {
232 for (Points p(*m_grid); p; p.next()) {
233 const int i = p.i(), j = p.j();
234
235 gradient(i, j) = 2*x(i, j) / m_normalization;
236 }
237 }
238 }
239
240 } // end of namespace inverse
241 } // end of namespace pism