tIPTaoTikhonovProblem.hh - 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
---
tIPTaoTikhonovProblem.hh (16956B)
---
1 // Copyright (C) 2012,2013,2014,2015,2016,2017 David Maxwell and Constantine Khroulev
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 #ifndef IPTAOTIKHONOVPROBLEM_HH_4NMM724B
20 #define IPTAOTIKHONOVPROBLEM_HH_4NMM724B
21
22 #include <memory>
23
24 #include "TaoUtil.hh"
25 #include "functional/IPFunctional.hh"
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/IceGrid.hh"
28 #include "pism/util/Logger.hh"
29
30 namespace pism {
31 namespace inverse {
32
33 template<class ForwardProblem> class IPTaoTikhonovProblem;
34
35 //! Iteration callback class for IPTaoTikhonovProblem
36 /** A class for objects receiving iteration callbacks from a
37 * IPTaoTikhonovProblem. These callbacks can be used to monitor the
38 * solution, plot iterations, print diagnostic messages, etc.
39 * IPTaoTikhonovProblemListeners are ususally used via a reference
40 * counted pointer IPTaoTikhonovProblemListener::Ptr to allow for good
41 * memory management when Listeners are created as subclasses of
42 * python classes. It would have been better to nest this inside of
43 * IPTaoTikhonovProblem, but SWIG has a hard time with nested classes,
44 * so it's outer instead.
45 */
46 template<class ForwardProblem> class IPTaoTikhonovProblemListener {
47 public:
48 typedef std::shared_ptr<IPTaoTikhonovProblemListener> Ptr;
49
50
51 typedef typename ForwardProblem::DesignVec::Ptr DesignVecPtr;
52 typedef typename ForwardProblem::StateVec::Ptr StateVecPtr;
53
54 IPTaoTikhonovProblemListener() {}
55 virtual ~IPTaoTikhonovProblemListener() {}
56
57 //! The method called after each minimization iteration.
58 virtual void iteration(IPTaoTikhonovProblem<ForwardProblem> &problem,
59 double eta, int iter,
60 double objectiveValue, double designValue,
61 const DesignVecPtr d,
62 const DesignVecPtr diff_d,
63 const DesignVecPtr grad_d,
64 const StateVecPtr u,
65 const StateVecPtr diff_u,
66 const DesignVecPtr grad_u,
67 const DesignVecPtr gradient) = 0;
68 };
69
70
71 //! \brief Defines a Tikhonov minimization problem to be solved with a TaoBasicSolver.
72 /*! Suppose \f$F\f$ is a map from a space \f$D\f$ of design variables to a space \f$S\f$ of
73 state variables and we wish to solve a possibly ill-posed problem of the form
74 \f[ F(d) = u \f]
75 where \f$u\f$ is know and \f$d\f$ is unknown. Approximate solutions can be obtained by
76 finding minimizers of an associated Tikhonov functional
77 \f[
78 J(d) = J_{S}(F(d)-u) + \frac{1}{\eta}J_{D}(d-d_0)
79 \f]
80 where \$J_{D}\$ and \$J_{S}\$ are functionals on the spaces \f$D\f$ and \f$S\f$ respectively,
81 \f$\eta\f$ is a penalty parameter, and \f$d_0\f$ is a best a-priori guess for the the solution.
82 The IPTaoTikhonovProblem class encapuslates all of the data required to formulate the minimization
83 problem as a Problem tha can be solved using a TaoBasicSolver. It is templated on the
84 the class ForwardProblem which defines the class of the forward map \f$F\f$ as well as the
85 spaces \f$D\f$ and \f$S\f$. An instance of ForwardProblem, along with
86 specific functionals \f$J_D\f$ and \f$J_S\f$, the parameter \f$\eta\f$, and the data
87 \f$y\f$ and \f$x_0\f$ are provided on constructing a IPTaoTikhonovProblem.
88
89 For example, if the SSATaucForwardProblem class defines the map taking yield stresses \f$\tau_c\f$
90 to the corresponding surface velocity field solving the SSA, a schematic setup of solving
91 the associated Tikhonov problem goes as follows.
92
93 \code
94 SSATaucForwardProblem forwardProblem(grid);
95 L2NormFunctional2S designFunctional(grid); //J_X
96 L2NormFunctional2V stateFunctional(grid); //J_Y
97 IceModelVec2V u_obs; // Set this to the surface velocity observations.
98 IceModelVec2S tauc_0; // Set this to the initial guess for tauc.
99 double eta; // Set this to the desired penalty parameter.
100
101 typedef InvSSATauc IPTaoTikhonovProblem<SSATaucForwardProblem>;
102 InvSSATauc tikhonovProblem(forwardProblem,tauc_0,u_obs,eta,designFunctional,stateFunctional);
103
104 TaoBasicSolver<InvSSATauc> solver(com, "cg", tikhonovProblem);
105
106 TerminationReason::Ptr reason = solver.solve();
107
108 if (reason->succeeded()) {
109 printf("Success: %s\n",reason->description().c_str());
110 } else {
111 printf("Failure: %s\n",reason->description().c_str());
112 }
113 \endcode
114
115 The class ForwardProblem that defines the forward problem
116 must have the following characteristics:
117
118 <ol>
119 <li> Contains typedefs for DesignVec and StateVec that effectively
120 define the function spaces \f$D\f$ and \f$S\f$. E.g.
121
122 \code
123 typedef IceModelVec2S DesignVec;
124 typedef IceModelVec2V StateVec;
125 \endcode
126 would be appropriate for a map from basal yeild stress to surface velocities.
127
128 <li> A method
129 \code
130 TerminationReason::Ptr linearize_at(DesignVec &d);
131 \endcode
132 that instructs the class to compute the value of F and
133 anything needed to compute its linearization at \a d. This is the first method
134 called when working with a new iterate of \a d.
135
136 <li> A method
137 \code
138 StateVec &solution()
139 \endcode
140 that returns the most recently computed value of \f$F(d)\f$
141 as computed by a call to linearize_at.
142
143 <li> A method
144 \code
145 void apply_linearization_transpose(StateVec &du, DesignVec &dzeta);
146 \endcode
147 that computes the action of \f$(F')^t\f$,
148 where \f$F'\f$ is the linearization of \f$F\f$ at the current iterate, and the transpose
149 is computed in the standard sense (i.e. thinking of \f$F'\f$ as a matrix with respect
150 to the bases implied by the DesignVec and StateVec spaces). The need for a transpose arises
151 because
152 \f[
153 \frac{d}{dt} J_{S}(F(d+t\delta d)-u) = [DJ_S]_{k}\; F'_{kj} \; \delta d
154 \f]
155 and hence the gradient of the term \f$J_{S}(F(d)-u)\f$ with respect to \f$d\f$ is given
156 by
157 \f[
158 (F')^t (\nabla J_S)^t.
159 \f]
160 </ol>
161 */
162 template<class ForwardProblem> class IPTaoTikhonovProblem
163 {
164 public:
165
166 typedef typename ForwardProblem::DesignVec DesignVec;
167 typedef typename ForwardProblem::StateVec StateVec;
168
169 typedef typename ForwardProblem::DesignVec::Ptr DesignVecPtr;
170 typedef typename ForwardProblem::StateVec::Ptr StateVecPtr;
171
172 /*! Constructs a Tikhonov problem:
173
174 Minimize \f$J(d) = J_S(F(d)-u_obs) + \frac{1}{\eta} J_D(d-d0) \f$
175
176 that can be solved with a TaoBasicSolver.
177
178 @param forward Class defining the map F. See class-level documentation for requirements of F.
179 @param d0 Best a-priori guess for the design parameter.
180 @param u_obs State parameter to match (i.e. approximately solve F(d)=u_obs)
181 @param eta Penalty parameter/Lagrange multiplier. Take eta to zero to impose more regularization to an ill posed problem.
182 @param designFunctional The functional \f$J_D\f$
183 @param stateFunctional The functional \f$J_S\f$
184 */
185
186 IPTaoTikhonovProblem(ForwardProblem &forward, DesignVec &d0, StateVec &u_obs, double eta,
187 IPFunctional<DesignVec>&designFunctional, IPFunctional<StateVec>&stateFunctional);
188
189 virtual ~IPTaoTikhonovProblem();
190
191
192 //! Sets the initial guess for minimization iterations. If this isn't set explicitly,
193 // the parameter \f$d0\f$ appearing the in the Tikhonov functional will be used.
194 virtual void setInitialGuess(DesignVec &d) {
195 m_dGlobal.copy_from(d);
196 }
197
198 //! Callback provided to TAO for objective evaluation.
199 virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient);
200
201 //! Add an object to the list of objects to be called after each iteration.
202 virtual void addListener(typename IPTaoTikhonovProblemListener<ForwardProblem>::Ptr listener) {
203 m_listeners.push_back(listener);
204 }
205
206 //! Final value of \f$F(d)\f$, where \f$d\f$ is the solution of the minimization.
207 virtual StateVecPtr stateSolution() {
208 return m_forward.solution();
209 }
210
211 //! Value of \f$d\f$, the solution of the minimization problem.
212 virtual DesignVecPtr designSolution() {
213 return m_d;
214 }
215
216 //! Callback from TaoBasicSolver, used to wire the connections between a Tao and
217 // the current class.
218 virtual void connect(Tao tao);
219
220 //! Callback from TAO after each iteration. The call is forwarded to each element of our list of listeners.
221 virtual void monitorTao(Tao tao);
222
223 //! Callback from TAO to detect convergence. Allows us to implement a custom convergence check.
224 virtual void convergenceTest(Tao tao);
225
226 //! Callback from TaoBasicSolver to form the starting iterate for the minimization. See also
227 // setInitialGuess.
228 virtual TerminationReason::Ptr formInitialGuess(Vec *v) {
229 *v = m_dGlobal.vec();
230 return GenericTerminationReason::success();
231 }
232
233 protected:
234
235 IceGrid::ConstPtr m_grid;
236
237 ForwardProblem &m_forward;
238
239 /// Current iterate of design parameter
240 DesignVecPtr m_d;
241 /// Initial iterate of design parameter, stored without ghosts for the benefit of TAO.
242 DesignVec m_dGlobal;
243 /// A-priori estimate of design parameter
244 DesignVec &m_d0;
245 /// Storage for (m_d-m_d0)
246 DesignVecPtr m_d_diff;
247
248 /// State parameter to match via F(d)=u_obs
249 StateVec &m_u_obs;
250 /// Storage for F(d)-u_obs
251 StateVecPtr m_u_diff;
252
253 /// Temporary storage used in gradient computation.
254 StateVec m_adjointRHS;
255
256 /// Gradient of \f$J_D\f$ at the current iterate.
257 DesignVecPtr m_grad_design;
258 /// Gradient of \f$J_S\f$ at the current iterate.
259 DesignVecPtr m_grad_state;
260 /** Weighted sum of the design and state gradients corresponding to
261 * the gradient of the Tikhonov functional \f$J\f$.
262 */
263 DesignVecPtr m_grad;
264
265 /// Penalty parameter/Lagrange multiplier.
266 double m_eta;
267
268 /// Value of \f$J_D\f$ at the current iterate.
269 double m_val_design;
270 /// Value of \f$J_S\f$ at the current iterate.
271 double m_val_state;
272
273 /// Implementation of \f$J_D\f$.
274 IPFunctional<IceModelVec2S> &m_designFunctional;
275 /// Implementation of \f$J_S\f$.
276 IPFunctional<IceModelVec2V> &m_stateFunctional;
277
278 /// List of iteration callbacks.
279 std::vector<typename IPTaoTikhonovProblemListener<ForwardProblem>::Ptr> m_listeners;
280
281 /// Convergence parameter: convergence stops when \f$||J_D||_2 <\f$ m_tikhonov_rtol.
282 double m_tikhonov_atol;
283
284 /** Convergence parameter: convergence stops when \f$||J_D||_2 \f$
285 * is less than m_tikhonov_rtol times the maximum of the gradient of
286 * \f$J_S\f$ and \f$(1/\eta)J_D\f$. This occurs when the two terms
287 * forming the sum of the gradient of \f$J\f$ point in roughly
288 * opposite directions with the same magnitude.
289 */
290 double m_tikhonov_rtol;
291
292 };
293
294 template<class ForwardProblem>
295 IPTaoTikhonovProblem<ForwardProblem>::IPTaoTikhonovProblem(ForwardProblem &forward,
296 DesignVec &d0, StateVec &u_obs,
297 double eta,
298 IPFunctional<DesignVec> &designFunctional,
299 IPFunctional<StateVec> &stateFunctional)
300 : m_forward(forward), m_d0(d0), m_u_obs(u_obs), m_eta(eta),
301 m_designFunctional(designFunctional), m_stateFunctional(stateFunctional) {
302
303 m_grid = m_d0.grid();
304
305 m_tikhonov_atol = m_grid->ctx()->config()->get_number("inverse.tikhonov.atol");
306 m_tikhonov_rtol = m_grid->ctx()->config()->get_number("inverse.tikhonov.rtol");
307
308 int design_stencil_width = m_d0.stencil_width();
309 int state_stencil_width = m_u_obs.stencil_width();
310
311 m_d.reset(new DesignVec);
312 m_d->create(m_grid, "design variable", WITH_GHOSTS, design_stencil_width);
313
314 m_dGlobal.create(m_grid, "design variable (global)", WITHOUT_GHOSTS, design_stencil_width);
315 m_dGlobal.copy_from(m_d0);
316
317 m_u_diff.reset(new StateVec);
318 m_u_diff->create(m_grid, "state residual", WITH_GHOSTS, state_stencil_width);
319
320 m_d_diff.reset(new DesignVec);
321 m_d_diff->create(m_grid, "design residual", WITH_GHOSTS, design_stencil_width);
322
323 m_grad_state.reset(new DesignVec);
324 m_grad_state->create(m_grid, "state gradient", WITHOUT_GHOSTS, design_stencil_width);
325
326 m_grad_design.reset(new DesignVec);
327 m_grad_design->create(m_grid, "design gradient", WITHOUT_GHOSTS, design_stencil_width);
328
329 m_grad.reset(new DesignVec);
330 m_grad->create(m_grid, "gradient", WITHOUT_GHOSTS, design_stencil_width);
331
332 m_adjointRHS.create(m_grid,"work vector", WITHOUT_GHOSTS, design_stencil_width);
333 }
334
335 template<class ForwardProblem>
336 IPTaoTikhonovProblem<ForwardProblem>::~IPTaoTikhonovProblem() {
337 // empty
338 }
339
340 template<class ForwardProblem>
341 void IPTaoTikhonovProblem<ForwardProblem>::connect(Tao tao) {
342 typedef taoutil::TaoObjGradCallback<IPTaoTikhonovProblem<ForwardProblem>,
343 &IPTaoTikhonovProblem<ForwardProblem>::evaluateObjectiveAndGradient> ObjGradCallback;
344
345 ObjGradCallback::connect(tao,*this);
346
347 taoutil::TaoMonitorCallback< IPTaoTikhonovProblem<ForwardProblem> >::connect(tao,*this);
348
349 taoutil::TaoConvergenceCallback< IPTaoTikhonovProblem<ForwardProblem> >::connect(tao,*this);
350
351 double
352 gatol = PETSC_DEFAULT,
353 grtol = PETSC_DEFAULT,
354 gttol = PETSC_DEFAULT;
355
356 #if PETSC_VERSION_LT(3,7,0)
357 double fatol = 1e-10, frtol = 1e-20;
358 PetscErrorCode ierr = TaoSetTolerances(tao, fatol, frtol, gatol, grtol, gttol);
359 PISM_CHK(ierr, "TaoSetTolerances");
360 #else
361 PetscErrorCode ierr = TaoSetTolerances(tao, gatol, grtol, gttol);
362 PISM_CHK(ierr, "TaoSetTolerances");
363 #endif
364 }
365
366 template<class ForwardProblem>
367 void IPTaoTikhonovProblem<ForwardProblem>::monitorTao(Tao tao) {
368 PetscInt its;
369 TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, NULL);
370
371 int nListeners = m_listeners.size();
372 for (int k=0; k<nListeners; k++) {
373 m_listeners[k]->iteration(*this, m_eta,
374 its, m_val_design, m_val_state,
375 m_d, m_d_diff, m_grad_design,
376 m_forward.solution(), m_u_diff, m_grad_state,
377 m_grad);
378 }
379 }
380
381 template<class ForwardProblem> void IPTaoTikhonovProblem<ForwardProblem>::convergenceTest(Tao tao) {
382 double designNorm, stateNorm, sumNorm;
383 double dWeight, sWeight;
384 dWeight = 1/m_eta;
385 sWeight = 1;
386
387 designNorm = m_grad_design->norm(NORM_2);
388 stateNorm = m_grad_state->norm(NORM_2);
389 sumNorm = m_grad->norm(NORM_2);
390
391 designNorm *= dWeight;
392 stateNorm *= sWeight;
393
394 if (sumNorm < m_tikhonov_atol) {
395 TaoSetConvergedReason(tao, TAO_CONVERGED_GATOL);
396 } else if (sumNorm < m_tikhonov_rtol*std::max(designNorm,stateNorm)) {
397 TaoSetConvergedReason(tao, TAO_CONVERGED_USER);
398 } else {
399 TaoDefaultConvergenceTest(tao, NULL);
400 }
401 }
402
403 template<class ForwardProblem>
404 void IPTaoTikhonovProblem<ForwardProblem>::evaluateObjectiveAndGradient(Tao tao, Vec x,
405 double *value, Vec gradient) {
406 PetscErrorCode ierr;
407 // Variable 'x' has no ghosts. We need ghosts for computation with the design variable.
408 m_d->copy_from_vec(x);
409
410 TerminationReason::Ptr reason = m_forward.linearize_at(*m_d);
411 if (reason->failed()) {
412 Logger::ConstPtr log = m_grid->ctx()->log();
413 log->message(2,
414 "IPTaoTikhonovProblem::evaluateObjectiveAndGradient"
415 " failure in forward solve\n%s\n", reason->description().c_str());
416 ierr = TaoSetConvergedReason(tao, TAO_DIVERGED_USER);
417 PISM_CHK(ierr, "TaoSetConvergedReason");
418 return;
419 }
420
421 m_d_diff->copy_from(*m_d);
422 m_d_diff->add(-1, m_d0);
423 m_designFunctional.gradientAt(*m_d_diff, *m_grad_design);
424
425 m_u_diff->copy_from(*m_forward.solution());
426 m_u_diff->add(-1, m_u_obs);
427
428 // The following computes the reduced gradient.
429 m_stateFunctional.gradientAt(*m_u_diff, m_adjointRHS);
430 m_forward.apply_linearization_transpose(m_adjointRHS, *m_grad_state);
431
432 m_grad->copy_from(*m_grad_design);
433 m_grad->scale(1.0 / m_eta);
434 m_grad->add(1, *m_grad_state);
435
436 ierr = VecCopy(m_grad->vec(), gradient);
437 PISM_CHK(ierr, "VecCopy");
438
439 double valDesign, valState;
440 m_designFunctional.valueAt(*m_d_diff, &valDesign);
441 m_stateFunctional.valueAt(*m_u_diff, &valState);
442
443 m_val_design = valDesign;
444 m_val_state = valState;
445
446 *value = valDesign / m_eta + valState;
447 }
448
449 } // end of namespace inverse
450 } // end of namespace pism
451
452 #endif /* end of include guard: IPTAOTIKHONOVPROBLEM_HH_4NMM724B */