tverify_ssa_inv.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
---
tverify_ssa_inv.py (4015B)
---
1 #! /usr/bin/env python3
2 #
3 # Copyright (C) 2012, 2014 David Maxwell
4 #
5 # This file is part of PISM.
6 #
7 # PISM is free software; you can redistribute it and/or modify it under the
8 # terms of the GNU General Public License as published by the Free Software
9 # Foundation; either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 # details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with PISM; if not, write to the Free Software
19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21 help_description = """Verifies that the end result of an SSA inversion using 'pismi.py'
22 terminates after an expected number of iterations and with an expected final misfit."""
23 help_epilog = """E.g.
24
25 verify_ssa_inv.py inversion_result.nc --desired_misfit 100 --iter_max 50 --morozov
26
27 will succeed if the inversion terminated with a misfit less than 100 and terminated using no more
28 than 50 iterations. The --morozov flag indicates that the next-to-last iteration should
29 checked to have a misfit of at least 100. On the other hand
30
31 verify_ssa_inv.py inversion_result.nc --desired_misfit 100 --iter_max 50 --misfit_tolerance 2.5
32
33 will succeed if the inversion terminated with a misfit within 2.5 of 100 and terminated using no more
34 than 50 iterations.
35
36 """
37
38
39 def fail(msg):
40 print("Test Failed: %s" % msg)
41 exit(1)
42
43
44 def success():
45 print("Test Succeeded.")
46 exit(0)
47
48
49 netCDF = None
50 try:
51 import netCDF4 as netCDF
52 except:
53 fail('Unable to import netCDF4.')
54 exit(1)
55
56 import argparse
57
58 parser = argparse.ArgumentParser(description=help_description,
59 epilog=help_epilog, formatter_class=argparse.RawDescriptionHelpFormatter)
60 parser.add_argument('inv_filename', help='NC file with inversion results', type=str)
61 parser.add_argument('-m', '--desired_misfit', required=True, type=float, help='desired final misfit')
62 parser.add_argument('-i', '--iter_max', required=True, type=int, help='maximum number of iterations')
63 parser.add_argument('-z', '--morozov', required=False, action='store_true',
64 help='verify that next-to-last misfit is larger than desired misfit')
65 parser.add_argument('-t', '--misfit_tolerance', required=False, type=float,
66 help='maximum difference between desired and computed misfit')
67
68 args = parser.parse_args()
69
70 inv_filename = args.inv_filename
71 try:
72 data = netCDF.Dataset(inv_filename)
73 except:
74 fail('Unable to open inversion file "%s"\nPerhaps the inversion run failed to converge, or terminated unexpectedly.' % inv_filename)
75
76 # Grab the misfit history from the file.
77 misfit = data.variables.get('inv_ssa_misfit')
78 if misfit is None:
79 fail('Inversion file %s missing misfit history variable "inv_ssa_misfit".\nPerhaps the inversion run failed to converge, or terminated unexpectedly.' % inv_filename)
80
81 desired_misfit = args.desired_misfit
82
83 iter_count = len(misfit)
84 last_misfit = misfit[-1]
85
86 if iter_count > args.iter_max:
87 fail("Inversion took an excessive number of iterations: %d > %d" % (iter_count, args.iter_max))
88
89 if args.misfit_tolerance:
90 d_misfit = abs(last_misfit - desired_misfit)
91 if d_misfit > args.misfit_tolerance:
92 fail("Final misfit too far from desired misfit: |final-desired| = |%0.4g-%0.4g| = %0.4g > %0.4g = desired" %
93 (last_misfit, desired_misfit, d_misfit, args.misfit_tolerance))
94 else:
95 if last_misfit >= desired_misfit:
96 fail("Desired final misfit not met: computed = %0.4g >= %0.4g = desired" % (last_misfit, desired_misfit))
97
98 if args.morozov:
99 penultimate_misfit = misfit[-2]
100 if penultimate_misfit < desired_misfit:
101 fail("Penultimate misfit too small: computed = %0.4g < %0.4g = desired" % (penultimate_misfit, desired_misfit))
102
103
104 success()