tssa_siple.py - 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
---
tssa_siple.py (24210B)
---
1 # Copyright (C) 2012, 2014, 2015, 2016, 2017, 2018, 2019 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 """Inverse SSA solvers using the siple library."""
20
21 import PISM
22 from PISM.logging import logError, logMessage
23 from PISM.util import convert, Bunch
24
25 if not PISM.imported_from_sphinx:
26 from petsc4py import PETSc
27
28 import sys
29 import math
30 import traceback
31
32 import siple
33 import PISM.invert.sipletools
34 from PISM.invert.ssa import InvSSASolver
35 from PISM.invert.sipletools import PISMLocalVector
36 from siple.gradient.forward import NonlinearForwardProblem
37 from siple.gradient.nonlinear import InvertNLCG, InvertIGN
38
39 siple.reporting.clear_loggers()
40 siple.reporting.add_logger(PISM.invert.sipletools.pism_logger)
41 siple.reporting.set_pause_callback(PISM.invert.sipletools.pism_pause)
42
43 WIDE_STENCIL = 2
44
45
46 class InvSSASolver_Gradient(InvSSASolver):
47
48 """Inverse SSA solver based on `siple` iterative gradient methods."""
49
50 def __init__(self, ssarun, method):
51 """
52 :param ssarun: The :class:`PISM.ssa.SSARun` defining the forward problem. See :class:`PISM.inv_ssa.SSAForwardRunFromInputFile`.
53 :param method: String describing the actual ``siple`` algorithm to use. One of ``sd``, ``nlcg`` or ``ign``."""
54 InvSSASolver.__init__(self, ssarun, method)
55
56 monitor_adjoint = PISM.OptionBool("-inv_monitor_adjoint",
57 "Track accuracy of the adjoint during computation")
58 morozov_scale_factor = PISM.OptionReal("-inv_morozov_scale",
59 "Scale factor (>=1) for Morozov discrepancy principle",
60 1.0).value()
61 ls_verbose = PISM.OptionBool("-inv_ls_verbose", "Turn on a verbose linesearch.")
62 ign_theta = PISM.OptionReal("-ign_theta", "theta parameter for IGN algorithm", 0.5)
63 max_it = int(self.config.get_number("inverse.max_iterations"))
64
65 target_misfit = PISM.OptionReal("-inv_target_misfit",
66 "m/year; desired root misfit for inversions", 0.0)
67
68 if target_misfit.is_set():
69 self.target_misfit = target_misfit.value()
70 else:
71 raise RuntimeError("Missing required option -inv_target_misfit")
72
73 velocity_scale = ssarun.grid.ctx().config().get_number("inverse.ssa.velocity_scale", "m/year")
74 self.target_misfit /= velocity_scale
75
76 self.forward_problem = SSAForwardProblem(ssarun)
77
78 # Determine the inversion algorithm, and set up its arguments.
79 if self.method == "ign":
80 Solver = InvertSSAIGN
81 else:
82 Solver = InvertSSANLCG
83
84 self.u_i = None
85 self.zeta_i = None
86
87 params = Solver.defaultParameters()
88 params.ITER_MAX = max_it
89 if self.method == "sd":
90 params.steepest_descent = True
91 params.ITER_MAX = 10000
92 elif self.method == "ign":
93 params.linearsolver.ITER_MAX = 10000
94 params.linearsolver.verbose = True
95 if ls_verbose:
96 params.linesearch.verbose = True
97 params.verbose = True
98 params.thetaMax = ign_theta
99 params.mu = morozov_scale_factor
100
101 params.deriv_eps = 0.
102
103 self.siple_solver = Solver(self.forward_problem, params=params)
104
105 if monitor_adjoint:
106 self.addIterationListener(MonitorAdjoint())
107 if self.method == 'ign':
108 self.addLinearIterationListener(MonitorAdjointLin())
109
110 def solveForward(self, zeta, out=None):
111 r"""Given a parameterized design variable value :math:`\zeta`, solve the SSA.
112 See :cpp:class:`IP_TaucParam` for a discussion of parameterizations.
113
114 :param zeta: :cpp:class:`IceModelVec` containing :math:`\zeta`.
115 :param out: optional :cpp:class:`IceModelVec` for storage of the computation result.
116 :returns: An :cpp:class:`IceModelVec` contianing the computation result.
117 """
118 if out is None:
119 out = self.forward_problem.F(PISMLocalVector(zeta))
120 else:
121 out = self.forward_problem.F(PISMLocalVector(zeta), out=PISMLocalVector(out))
122 return out.core()
123
124 def addIterationListener(self, listener):
125 """Add a listener to be called after each iteration. See FIXME."""
126 self.siple_solver.addIterationListener(SipleIterationListenerAdaptor(self, listener))
127
128 def addDesignUpdateListener(self, listener):
129 self.siple_solver.addXUpdateListener(SipleDesignUpdateListenerAdaptor(self, listener))
130
131 def addLinearIterationListener(self, listener):
132 "Add a linear iteration listener."
133 self.siple_solver.addLinearIterationListener(SipleLinearIterationListenerAdaptor(self, listener))
134
135 def solveInverse(self, zeta0, u_obs, zeta_inv):
136 r"""Executes the inversion algorithm.
137
138 :param zeta0: The best `a-priori` guess for the value of the parameterized design variable value :math:`\zeta`.
139 Ignored by siple gradient algorithms in deference to `zeta_inv`.
140 :param u_obs: :cpp:class:`IceModelVec2V` of observed surface velocities.
141 :param zeta_inv: :cpp:class:`zeta_inv` starting value of :math:`\zeta` for iterative minimization.
142 :returns: A :cpp:class:`TerminationReason`.
143 """
144 try:
145 vecs = self.ssarun.modeldata.vecs
146 if vecs.has('zeta_fixed_mask'):
147 self.ssarun.ssa.set_tauc_fixed_locations(vecs.zeta_fixed_mask)
148
149 (self.zeta_i, self.u_i) = self.siple_solver.solve(zeta_inv, u_obs, self.target_misfit)
150 except Exception:
151 import traceback
152 exc_type, exc_value, exc_traceback = sys.exc_info()
153 description = ""
154 for l in traceback.format_exception(exc_type, exc_value, exc_traceback):
155 description += l
156 # It would be nice to make siple so that if the inverse solve fails
157 # you can still keep the most recent iteration.
158 self.u_i = None
159 self.zeta_i = None
160 return PISM.GenericTerminationReason(-1, description)
161 return PISM.GenericTerminationReason(1, "Morozov Discrepancy Met")
162
163 def inverseSolution(self):
164 """Returns a tuple ``(zeta, u)`` of :cpp:class:`IceModelVec`'s corresponding to the values
165 of the design and state variables at the end of inversion."""
166 return (self.zeta_i, self.u_i)
167
168
169 class SSAForwardProblem(NonlinearForwardProblem):
170
171 """Subclass of a :class:`siple.NonlinearForwardProblem` defining the forward problem for ``siple``.
172 The heavy lifting is forwarded to a :cpp:class:`IP_SSATaucForwardProblem` or
173 :cpp:class:`IP_SSAHardnessForwardProblem` living inside a :class:`PISM.invert.ssa.SSATaucForwardRun`."""
174
175 def __init__(self, ssarun):
176 """
177 :param ssarun:
178
179 A :class:`PISM.invert.ssa.SSAForwardRun`
180 or :class:`PISM.invert.ssa.SSAForwardRunFromInputFile` defining
181 the forward problem."""
182
183 self.ssarun = ssarun
184 self.ssa = self.ssarun.ssa
185 self.grid = ssarun.grid
186
187 self.tmpV = PISM.IceModelVec2V()
188 self.tmpV.create(self.grid, "work vector (2V)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)
189
190 self.tmpS = PISM.IceModelVec2S()
191 self.tmpS.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)
192
193 self.tmpS2 = PISM.IceModelVec2S()
194 self.tmpS2.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL)
195
196 ksp_rtol = 1e-12
197 self.ksp = PETSc.KSP()
198 self.ksp.create(self.grid.com)
199 self.ksp.setTolerances(ksp_rtol, PETSc.DEFAULT, PETSc.DEFAULT, PETSc.DEFAULT)
200 self.ksp.getPC().setType('bjacobi')
201 self.ksp.setFromOptions()
202
203 (self.designFunctional, self.stateFunctional) = PISM.invert.ssa.createGradientFunctionals(ssarun)
204
205 self.designForm = None
206
207 def F(self, x, out=None, guess=None):
208 """
209 Returns the value of the forward problem at the design variable ``x``.
210
211 Nonlinear problems often make use of an initial guess; this can be provided in ``guess``.
212
213 Storage in ``out``, if given, is used for the return value.
214 """
215 if out is None:
216 out = self.rangeVector()
217 reason = self.ssa.linearize_at(x.core())
218 if reason.failed():
219 raise PISM.AlgorithmFailureException(reason.description())
220 out.core().copy_from(self.ssa.solution())
221 return out
222
223 def T(self, d, out=None):
224 """
225 Returns the value of the linearization :math:`T` of the forward problem at the design variable :math:`x`
226 specified previously in :meth:`linearizeAt`, in the direction `d`.
227
228 Storage in `out`, if given, is used for the return value.
229 """
230 if out is None:
231 out = self.rangeVector()
232 self.ssa.apply_linearization(d.core(), out.core())
233 return out
234
235 def TStar(self, r, out=None):
236 """
237 Let :math:`T` be the linearization of the forward problem :math:`F` at the design variable `x` (as specified previously in :meth:`linearizeAt`).
238 Its adjoint is :math:`T^*`. This method returns the value of :math:`T^*` in the direction `r`.
239
240 Storage in `out`, if given, is used for the return value.
241 """
242 if out is None:
243 out = self.domainVector()
244
245 if self.designForm is None:
246 stencil_width = int(self.grid.ctx().config().get_number("grid.max_stencil_width"))
247 da2 = self.grid.get_dm(1, stencil_width)
248
249 if PISM.PETSc.Sys.getVersion() < (3, 5, 0):
250 self.designForm = da2.get().getMatrix("baij")
251 else:
252 da2.get().setMatType("baij")
253 self.designForm = da2.get().getMatrix()
254
255 self.designFunctional.assemble_form(self.designForm)
256
257 # First step
258 self.stateFunctional.gradientAt(r.core(), self.tmpV)
259 self.tmpV.scale(0.5)
260 self.ssa.apply_linearization_transpose(self.tmpV, self.tmpS)
261
262 # Second step
263 if PISM.PETSc.Sys.getVersion() < (3, 5, 0):
264 self.ksp.setOperators(self.designForm, self.designForm,
265 PETSc.Mat.Structure.SAME_NONZERO_PATTERN)
266 else:
267 self.ksp.setOperators(self.designForm, self.designForm)
268
269 self.ksp.solve(self.tmpS.vec(), self.tmpS2.vec())
270
271 reason = self.ksp.getConvergedReason()
272 if reason < 0:
273 raise RuntimeError('TStarB linear solve failed to converge (KSP reason %s)\n\n' % reason)
274 else:
275 PISM.logging.logPrattle("TStarB converged (KSP reason %s)\n" % reason)
276
277 out.core().copy_from(self.tmpS2)
278 return out
279
280 def linearizeAt(self, x, guess=None):
281 """
282 Instructs the class that subsequent calls to T and TStar will be conducted for the given value of x.
283
284 Nonlinear problems often make use of an initial guess; this can be provided in 'guess'.
285 """
286 reason = self.ssa.linearize_at(x.core())
287 if reason.failed():
288 raise Exception(reason.description())
289
290 def evalFandLinearize(self, x, out=None, guess=None):
291 """
292 Computes the value of F(x) and locks in a linearization. Sometimes there are efficiencies that
293 can be acheived this way.
294
295 Default implementation simply calls F, then linearizeAt.
296 """
297 if out is None:
298 out = self.rangeVector()
299 self.linearizeAt(x)
300 out.core().copy_from(self.ssa.solution())
301 return out
302
303 def rangeIP(self, a, b):
304 """
305 Computes the inner product of two vectors in the range (i.e. state) space.
306 """
307 return self.stateFunctional.dot(a.core(), b.core())
308
309 def domainIP(self, a, b):
310 """
311 Computes the inner product of two vectors in the domain (i.e. design) space.
312 """
313 return self.designFunctional.dot(a.core(), b.core())
314
315 def rangeVector(self):
316 """Constructs a brand new vector from the range (i.e. state) vector space"""
317 v = PISM.IceModelVec2V()
318 v.create(self.grid, "", True, WIDE_STENCIL)
319
320 # Add appropriate meta data.
321 intent = "?inverse?" # FIXME
322 desc = "SSA velocity computed by inversion"
323 v.set_attrs(intent, "%s%s" % ("X-component of the ", desc), "m s-1", "m year-1", "", 0)
324 v.set_attrs(intent, "%s%s" % ("Y-component of the ", desc), "m s-1", "m year-1", "", 1)
325
326 huge_vel = convert(1e6, "m/year", "m/second")
327 attrs = [("valid_min", -huge_vel), ("valid_max", huge_vel), ("_FillValue", 2 * huge_vel)]
328 for a in attrs:
329 for component in range(2):
330 v.metadata(component).set_number(a[0], a[1])
331
332 return PISMLocalVector(v)
333
334 def domainVector(self):
335 """Constructs a brand new vector from the domain (i.e. design) vector space"""
336 v = PISM.IceModelVec2S()
337 v.create(self.grid, "", True, WIDE_STENCIL)
338 return PISMLocalVector(v)
339
340
341 class InvertSSANLCG(InvertNLCG):
342
343 r"""Subclass of :class:`siple.gradient.nonlinear.InvertNLCG` for inversion of SSA velocities
344 from :math:`\tau_c` or hardness."""
345
346 @staticmethod
347 def defaultParameters():
348 params = InvertNLCG.defaultParameters()
349 return params
350
351 def __init__(self, forward_problem, params=None):
352 """:param forward_problem: A :class:`SSATaucForwardProblem` that defines the forward problem.
353 :param params: A :class:`siple.params.Parameters` containing run-time parameters.
354 """
355 InvertNLCG.__init__(self, params)
356 self.forward_problem = forward_problem
357 self.misfit_goal = 0.0
358
359 def forwardProblem(self):
360 """:returns: the associated :class:`SSATaucForwardProblem` provided at construction"""
361 return self.forward_problem
362
363 def stopConditionMet(self, count, x, Fx, y, r):
364 """
365 Determines if minimization should be halted (based, e.g. on a Morozov discrepancy principle)
366
367 :param count: current iteration count
368 :param x: point in domain of potential minimizer.
369 :param Fx: value of nonlinear function at `x`
370 :param r: current residual, i.e. :math:`y-F(x)`
371 :returns: boolean, `True` if termination condition is met
372 """
373
374 misfit = math.sqrt(abs(self.forward_problem.rangeIP(r, r)))
375
376 if (misfit < self.misfit_goal):
377 siple.reporting.msg('Stop condition met')
378 return True
379 return False
380
381 def initialize(self, x, y, deltaLInf):
382 """
383 Hook called at the start of :meth:`solve`. This gives the class a chance to massage the input.
384
385 The remaining arguments are passed directly from solve, and can be used for determining the
386 final stopping criterion.
387
388 Returns vectors corresponding to the initial value of `x` and the desired value of `y` where :math:`y=F(x)`.
389 """
390 xv = PISMLocalVector(x)
391 yv = PISMLocalVector(y)
392
393 self.misfit_goal = self.params.mu * deltaLInf
394
395 return (xv, yv)
396
397 def finalize(self, x, y):
398 """
399 Hook called at the end of :meth:`solve`. Gives the chance to massage the return values.
400 """
401 zeta = x.core()
402 u = y.core()
403 return (zeta, u)
404
405
406 class InvertSSAIGN(InvertIGN):
407
408 r"""Subclass of :class:`siple.gradient.nonlinear.InvertIGN` for inversion of SSA velocities
409 from :math:`\tau_c` or hardness."""
410
411 @staticmethod
412 def defaultParameters():
413 params = InvertIGN.defaultParameters()
414 return params
415
416 def __init__(self, forward_problem, params=None):
417 """:param forward_problem: A :class:`SSATaucForwardProblem` that defines the forward problem.
418 :param params: A :class:`siple.params.Parameters` containing run-time parameters.
419 """
420 InvertIGN.__init__(self, params)
421 self.forward_problem = forward_problem
422
423 def forwardProblem(self):
424 """:returns: the associated :class:`SSATaucForwardProblem` provided at construction"""
425 return self.forward_problem
426
427 def temper_d(self, x, d, y, r):
428 """Method called during iterations in :meth:`solve` to ensure we don't make to large a step.
429 :param x: current design parameter
430 :param d: step in design space
431 :param y: desired value of state paramter
432 :param r: residual of desired and current state parameter values (:math:`y-F(x)`)
433
434 Changes `d` as a side-effect to temper step size.
435 """
436 dnorm = d.norm('linf')
437 xnorm = x.norm('linf')
438 if dnorm > 2 * xnorm:
439 siple.reporting.msg('wild change predicted by linear step. scaling')
440 d.scale(2 * xnorm / dnorm)
441
442 def initialize(self, x, y, target_misfit):
443 """
444 Hook called at the start of :meth:`solve`. This gives the class a chance to massage the input.
445
446 The remaining arguments are passed directly from solve, and can be used for determining the
447 final stopping criterion.
448
449 Returns vectors corresponding to the initial value of `x` and the desired value of `y` where :math:`y=F(x)`.
450 """
451 xv = PISMLocalVector(x)
452 yv = PISMLocalVector(y)
453
454 return (xv, yv, target_misfit)
455
456 def finalize(self, x, y):
457 """
458 Hook called at the end of :meth:`solve`. Gives the chance to massage the return values.
459 """
460
461 zeta = x.core()
462 u = y.core()
463
464 return (zeta, u)
465
466
467 class SipleIterationListenerAdaptor(object):
468
469 """Adaptor for passing listening events from `siple`-based solvers to a python object."""
470
471 def __init__(self, owner, listener):
472 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us
473 :param listener: The python-based listener.
474 """
475 self.owner = owner
476 self.listener = listener
477
478 def __call__(self, siplesolver, it, x, Fx, y, d, r, *args):
479 """Callback from `siple`. Gathers together the long list of arguments
480 into a dictionary and passes it along in a standard form to the python listener."""
481 data = Bunch(zeta=x.core(), u=Fx.core(), zeta_step=d.core(), residual=r.core(), target_misfit=self.owner.target_misfit)
482
483 if self.owner.method == 'ign':
484 data.update(T_zeta_step=args[0].core())
485 else:
486 data.update(TStar_residual=args[0].core())
487 try:
488 self.listener(self.owner, it, data)
489 except Exception:
490 logError("\nERROR: Exception occured during an inverse solver listener callback:\n\n")
491 traceback.print_exc(file=sys.stdout)
492 raise
493
494 class SipleLinearIterationListenerAdaptor(object):
495
496 """Adaptor for passing listening events the linear steps of `siple`-based `ign` solvers to a python object."""
497
498 def __init__(self, owner, listener):
499 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us
500 :param listener: The python-based listener.
501 """
502 self.owner = owner
503 self.listener = listener
504
505 def __call__(self, siplesolver, it, x, y, d, r, Td, TStarR):
506 """Callback from `siple`. Gathers together the long list of arguments
507 into a dictionary and passes it along in a standard form to the python listener."""
508
509 data = Bunch(x=x.core(), y=y.core(), r=r.core(), d=d.core(), Td=Td.core(), TSTarR=TStarR.core())
510 try:
511 self.listener(self.owner, it, data)
512 except Exception:
513 logError("\nERROR: Exception occured during an inverse solver listener callback:\n\n")
514 traceback.print_exc(file=sys.stdout)
515 raise
516
517 class SipleDesignUpdateListenerAdaptor(object):
518
519 """Adaptor for design variable update events of `siple`-based solvers to a python object."""
520
521 def __init__(self, owner, listener):
522 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us
523 :param listener: The python-based listener.
524 """
525 self.owner = owner
526 self.listener = listener
527
528 def __call__(self, siplesolver, it, zeta, u, u_obs, r):
529 """Callback from `siple`. Gathers together the long list of arguments
530 into a dictionary and passes it along in a standard form to the python listener."""
531 data = Bunch(zeta=zeta.core(), u=u.core(), r=r.core(), u_obs=u_obs.core())
532 try:
533 self.listener(self.owner, it, data)
534 except Exception as e:
535 logError("\nWARNING: Exception occured during an inverse solver DesignUpdate listener callback:\n%s\n\n" % str(e))
536
537
538 class MonitorAdjoint(object):
539
540 r"""Iteration listener that can be used to verify the correctness of the implementation of an adjoint.
541 For adjoint-based interative inverse methods, a residual ``r`` is known in state space and a step direction ``d`` is
542 known in design space. A linearized forward problem :math:`T` maps from design space to state space, and its adjoint :math:`T^*`
543 goes in the opposite direction. The inner products :math:`\left<Td,r\right>_{\rm State}`
544 and :math:`\left<d,T^*r\right>_{\rm Design}` should always be the same; this is a good diagnostic to determine
545 of an adjoint has been coded correctly. The listener prints a comparison of the values of the two inner products
546 at each iteration.
547 """
548
549 def __init__(self):
550 self.Td = None
551 self.TStarR = None
552 self.didWarning = False
553
554 def __call__(self, inverse_solver, count, data):
555 """
556 :param inverse_sovler: the solver (e.g. :class:`~InvSolver_Tikhonov`) we are listening to.
557 :param count: the iteration number.
558 :param data: dictionary of data related to the iteration.
559 """
560 method = inverse_solver.method
561 if method != 'sd' and method != 'nlcg' and method != 'ign':
562 if not self.didWarning:
563 PISM.verbPrintf(1, PISM.Context().com, '\nWarning: unable to monitor adjoint for inverse method: %s\nOption -inv_monitor_adjoint ignored\n' % method)
564 self.didWarning = True
565 return
566 fp = inverse_solver.forward_problem
567 d = PISM.invert.sipletools.PISMLocalVector(data.d)
568 r = PISM.invert.sipletools.PISMLocalVector(data.r)
569 self.Td = fp.T(d, self.Td)
570 self.TStarR = fp.TStar(r, out=self.TStarR)
571 ip1 = fp.domainIP(d, self.TStarR)
572 ip2 = fp.rangeIP(self.Td, r)
573 logMessage("adjoint test: <Td,r>=%g <d,T^*r>=%g (percent error %g)", ip1, ip2, (abs(ip1 - ip2)) / max(abs(ip1), abs(ip2)))
574
575
576 class MonitorAdjointLin(object):
577
578 def __init__(self):
579 self.Td = None
580 self.TStarR = None
581
582 def __call__(self, inverse_solver, count, data):
583 """
584 :param inverse_sovler: the solver (e.g. :class:`~InvSSASolver_Tikhonov`) we are listening to.
585 :param count: the iteration number.
586 :param data: dictionary of data related to the iteration.
587 """
588 fp = inverse_solver.forward_problem
589 r = PISM.invert.sipletools.PISMLocalVector(data.r)
590 d = PISM.invert.sipletools.PISMLocalVector(data.d)
591 self.Td = fp.T(d, self.Td)
592 self.TStarR = fp.TStar(r, out=self.TStarR)
593 ip1 = fp.domainIP(d, self.TStarR)
594 ip2 = fp.rangeIP(self.Td, r)
595 logMessage("adjoint test: <Td,r>=%g <d,T^*r>=%g (percent error %g)", ip1, ip2, (abs(ip1 - ip2)) / max(abs(ip1), abs(ip2)))